1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5from __future__ import absolute_import, print_function 

6 

7import math 

8import random 

9import logging 

10 

11try: 

12 from StringIO import StringIO as BytesIO 

13except ImportError: 

14 from io import BytesIO 

15 

16import numpy as num 

17 

18from pyrocko.guts import (Object, Float, Bool, Int, Tuple, String, List, 

19 Unicode, Dict) 

20from pyrocko.guts_array import Array 

21from pyrocko.dataset import topo 

22from pyrocko import orthodrome as od 

23from . import gmtpy 

24from . import nice_value 

25 

26try: 

27 newstr = unicode 

28except NameError: 

29 newstr = str 

30 

31points_in_region = od.points_in_region 

32 

33logger = logging.getLogger('pyrocko.plot.automap') 

34 

35earthradius = 6371000.0 

36r2d = 180./math.pi 

37d2r = 1./r2d 

38km = 1000. 

39d2m = d2r*earthradius 

40m2d = 1./d2m 

41cm = gmtpy.cm 

42 

43 

44def darken(c, f=0.7): 

45 return (c[0]*f, c[1]*f, c[2]*f) 

46 

47 

48def corners(lon, lat, w, h): 

49 ll_lat, ll_lon = od.ne_to_latlon(lat, lon, -0.5*h, -0.5*w) 

50 ur_lat, ur_lon = od.ne_to_latlon(lat, lon, 0.5*h, 0.5*w) 

51 return ll_lon, ll_lat, ur_lon, ur_lat 

52 

53 

54def extent(lon, lat, w, h, n): 

55 x = num.linspace(-0.5*w, 0.5*w, n) 

56 y = num.linspace(-0.5*h, 0.5*h, n) 

57 slats, slons = od.ne_to_latlon(lat, lon, y[0], x) 

58 nlats, nlons = od.ne_to_latlon(lat, lon, y[-1], x) 

59 south = slats.min() 

60 north = nlats.max() 

61 

62 wlats, wlons = od.ne_to_latlon(lat, lon, y, x[0]) 

63 elats, elons = od.ne_to_latlon(lat, lon, y, x[-1]) 

64 elons = num.where(elons < wlons, elons + 360., elons) 

65 

66 if elons.max() - elons.min() > 180 or wlons.max() - wlons.min() > 180.: 

67 west = -180. 

68 east = 180. 

69 else: 

70 west = wlons.min() 

71 east = elons.max() 

72 

73 return topo.positive_region((west, east, south, north)) 

74 

75 

76class NoTopo(Exception): 

77 pass 

78 

79 

80class OutOfBounds(Exception): 

81 pass 

82 

83 

84class FloatTile(Object): 

85 xmin = Float.T() 

86 ymin = Float.T() 

87 dx = Float.T() 

88 dy = Float.T() 

89 data = Array.T(shape=(None, None), dtype=float, serialize_as='table') 

90 

91 def __init__(self, xmin, ymin, dx, dy, data): 

92 Object.__init__(self, init_props=False) 

93 self.xmin = float(xmin) 

94 self.ymin = float(ymin) 

95 self.dx = float(dx) 

96 self.dy = float(dy) 

97 self.data = data 

98 self._set_maxes() 

99 

100 def _set_maxes(self): 

101 self.ny, self.nx = self.data.shape 

102 self.xmax = self.xmin + (self.nx-1) * self.dx 

103 self.ymax = self.ymin + (self.ny-1) * self.dy 

104 

105 def x(self): 

106 return self.xmin + num.arange(self.nx) * self.dx 

107 

108 def y(self): 

109 return self.ymin + num.arange(self.ny) * self.dy 

110 

111 def get(self, x, y): 

112 ix = int(round((x - self.xmin) / self.dx)) 

113 iy = int(round((y - self.ymin) / self.dy)) 

114 if 0 <= ix < self.nx and 0 <= iy < self.ny: 

115 return self.data[iy, ix] 

116 else: 

117 raise OutOfBounds() 

118 

119 

120class City(Object): 

121 def __init__(self, name, lat, lon, population=None, asciiname=None): 

122 name = newstr(name) 

123 lat = float(lat) 

124 lon = float(lon) 

125 if asciiname is None: 

126 asciiname = name.encode('ascii', errors='replace') 

127 

128 if population is None: 

129 population = 0 

130 else: 

131 population = int(population) 

132 

133 Object.__init__(self, name=name, lat=lat, lon=lon, 

134 population=population, asciiname=asciiname) 

135 

136 name = Unicode.T() 

137 lat = Float.T() 

138 lon = Float.T() 

139 population = Int.T() 

140 asciiname = String.T() 

141 

142 

143class Map(Object): 

144 lat = Float.T(optional=True) 

145 lon = Float.T(optional=True) 

146 radius = Float.T(optional=True) 

147 width = Float.T(default=20.) 

148 height = Float.T(default=14.) 

149 margins = List.T(Float.T()) 

150 illuminate = Bool.T(default=True) 

151 skip_feature_factor = Float.T(default=0.02) 

152 show_grid = Bool.T(default=False) 

153 show_topo = Bool.T(default=True) 

154 show_scale = Bool.T(default=False) 

155 show_topo_scale = Bool.T(default=False) 

156 show_center_mark = Bool.T(default=False) 

157 show_rivers = Bool.T(default=True) 

158 show_plates = Bool.T(default=False) 

159 show_plate_velocities = Bool.T(default=False) 

160 show_plate_names = Bool.T(default=False) 

161 show_boundaries = Bool.T(default=False) 

162 illuminate_factor_land = Float.T(default=0.5) 

163 illuminate_factor_ocean = Float.T(default=0.25) 

164 color_wet = Tuple.T(3, Int.T(), default=(216, 242, 254)) 

165 color_dry = Tuple.T(3, Int.T(), default=(172, 208, 165)) 

166 color_boundaries = Tuple.T(3, Int.T(), default=(1, 1, 1)) 

167 topo_resolution_min = Float.T( 

168 default=40., 

169 help='minimum resolution of topography [dpi]') 

170 topo_resolution_max = Float.T( 

171 default=200., 

172 help='maximum resolution of topography [dpi]') 

173 replace_topo_color_only = FloatTile.T( 

174 optional=True, 

175 help='replace topo color while keeping topographic shading') 

176 topo_cpt_wet = String.T(default='light_sea') 

177 topo_cpt_dry = String.T(default='light_land') 

178 axes_layout = String.T(optional=True) 

179 custom_cities = List.T(City.T()) 

180 gmt_config = Dict.T(String.T(), String.T()) 

181 comment = String.T(optional=True) 

182 approx_ticks = Int.T(default=4) 

183 

184 def __init__(self, gmtversion='newest', **kwargs): 

185 Object.__init__(self, **kwargs) 

186 self._gmt = None 

187 self._scaler = None 

188 self._widget = None 

189 self._corners = None 

190 self._wesn = None 

191 self._minarea = None 

192 self._coastline_resolution = None 

193 self._rivers = None 

194 self._dems = None 

195 self._have_topo_land = None 

196 self._have_topo_ocean = None 

197 self._jxyr = None 

198 self._prep_topo_have = None 

199 self._labels = [] 

200 self._area_labels = [] 

201 self._gmtversion = gmtversion 

202 

203 def save(self, outpath, resolution=75., oversample=2., size=None, 

204 width=None, height=None, psconvert=False, crop_eps_mode=False): 

205 

206 ''' 

207 Save the image. 

208 

209 Save the image to ``outpath``. The format is determined by the filename 

210 extension. Formats are handled as follows: ``'.eps'`` and ``'.ps'`` 

211 produce EPS and PS, respectively, directly with GMT. If the file name 

212 ends with ``'.pdf'``, GMT output is fed through ``gmtpy-epstopdf`` to 

213 create a PDF file. For any other filename extension, output is first 

214 converted to PDF with ``gmtpy-epstopdf``, then with ``pdftocairo`` to 

215 PNG with a resolution oversampled by the factor ``oversample`` and 

216 finally the PNG is downsampled and converted to the target format with 

217 ``convert``. The resolution of rasterized target image can be 

218 controlled either by ``resolution`` in DPI or by specifying ``width`` 

219 or ``height`` or ``size``, where the latter fits the image into a 

220 square with given side length. To save transparency use 

221 ``psconvert=True``. To crop the output image with a rectangle to the 

222 nearest non-white element set ``crop_eps_mode=True``. 

223 ''' 

224 

225 gmt = self.gmt 

226 self.draw_labels() 

227 self.draw_axes() 

228 if self.show_topo and self.show_topo_scale: 

229 self._draw_topo_scale() 

230 

231 gmt.save(outpath, resolution=resolution, oversample=oversample, 

232 crop_eps_mode=crop_eps_mode, 

233 size=size, width=width, height=height, psconvert=psconvert) 

234 

235 @property 

236 def scaler(self): 

237 if self._scaler is None: 

238 self._setup_geometry() 

239 

240 return self._scaler 

241 

242 @property 

243 def wesn(self): 

244 if self._wesn is None: 

245 self._setup_geometry() 

246 

247 return self._wesn 

248 

249 @property 

250 def widget(self): 

251 if self._widget is None: 

252 self._setup() 

253 

254 return self._widget 

255 

256 @property 

257 def layout(self): 

258 if self._layout is None: 

259 self._setup() 

260 

261 return self._layout 

262 

263 @property 

264 def jxyr(self): 

265 if self._jxyr is None: 

266 self._setup() 

267 

268 return self._jxyr 

269 

270 @property 

271 def pxyr(self): 

272 if self._pxyr is None: 

273 self._setup() 

274 

275 return self._pxyr 

276 

277 @property 

278 def gmt(self): 

279 if self._gmt is None: 

280 self._setup() 

281 

282 if self._have_topo_ocean is None: 

283 self._draw_background() 

284 

285 return self._gmt 

286 

287 def _setup(self): 

288 if not self._widget: 

289 self._setup_geometry() 

290 

291 self._setup_lod() 

292 self._setup_gmt() 

293 

294 def _setup_geometry(self): 

295 wpage, hpage = self.width, self.height 

296 ml, mr, mt, mb = self._expand_margins() 

297 wpage -= ml + mr 

298 hpage -= mt + mb 

299 

300 wreg = self.radius * 2.0 

301 hreg = self.radius * 2.0 

302 if wpage >= hpage: 

303 wreg *= wpage/hpage 

304 else: 

305 hreg *= hpage/wpage 

306 

307 self._wreg = wreg 

308 self._hreg = hreg 

309 

310 self._corners = corners(self.lon, self.lat, wreg, hreg) 

311 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10) 

312 

313 x, y, z = ((west, east), (south, north), (-6000., 4500.)) 

314 

315 xax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks) 

316 yax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks) 

317 zax = gmtpy.Ax(mode='min-max', inc=1000., label='Height', 

318 scaled_unit='km', scaled_unit_factor=0.001) 

319 

320 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax)) 

321 

322 par = scaler.get_params() 

323 

324 west = par['xmin'] 

325 east = par['xmax'] 

326 south = par['ymin'] 

327 north = par['ymax'] 

328 

329 self._wesn = west, east, south, north 

330 self._scaler = scaler 

331 

332 def _setup_lod(self): 

333 w, e, s, n = self._wesn 

334 if self.radius > 1500.*km: 

335 coastline_resolution = 'i' 

336 rivers = False 

337 else: 

338 coastline_resolution = 'f' 

339 rivers = True 

340 

341 self._minarea = (self.skip_feature_factor * self.radius/km)**2 

342 

343 self._coastline_resolution = coastline_resolution 

344 self._rivers = rivers 

345 

346 self._prep_topo_have = {} 

347 self._dems = {} 

348 

349 cm2inch = gmtpy.cm/gmtpy.inch 

350 

351 dmin = 2.0 * self.radius * m2d / (self.topo_resolution_max * 

352 (self.height * cm2inch)) 

353 dmax = 2.0 * self.radius * m2d / (self.topo_resolution_min * 

354 (self.height * cm2inch)) 

355 

356 for k in ['ocean', 'land']: 

357 self._dems[k] = topo.select_dem_names(k, dmin, dmax, self._wesn) 

358 if self._dems[k]: 

359 logger.debug('using topography dataset %s for %s' 

360 % (','.join(self._dems[k]), k)) 

361 

362 def _expand_margins(self): 

363 if len(self.margins) == 0 or len(self.margins) > 4: 

364 ml = mr = mt = mb = 2.0 

365 elif len(self.margins) == 1: 

366 ml = mr = mt = mb = self.margins[0] 

367 elif len(self.margins) == 2: 

368 ml = mr = self.margins[0] 

369 mt = mb = self.margins[1] 

370 elif len(self.margins) == 4: 

371 ml, mr, mt, mb = self.margins 

372 

373 return ml, mr, mt, mb 

374 

375 def _setup_gmt(self): 

376 w, h = self.width, self.height 

377 scaler = self._scaler 

378 

379 if gmtpy.is_gmt5(self._gmtversion): 

380 gmtconf = dict( 

381 MAP_TICK_PEN_PRIMARY='1.25p', 

382 MAP_TICK_PEN_SECONDARY='1.25p', 

383 MAP_TICK_LENGTH_PRIMARY='0.2c', 

384 MAP_TICK_LENGTH_SECONDARY='0.6c', 

385 FONT_ANNOT_PRIMARY='12p,1,black', 

386 FONT_LABEL='12p,1,black', 

387 PS_CHAR_ENCODING='ISOLatin1+', 

388 MAP_FRAME_TYPE='fancy', 

389 FORMAT_GEO_MAP='D', 

390 PS_MEDIA='Custom_%ix%i' % ( 

391 w*gmtpy.cm, 

392 h*gmtpy.cm), 

393 PS_PAGE_ORIENTATION='portrait', 

394 MAP_GRID_PEN_PRIMARY='thinnest,0/50/0', 

395 MAP_ANNOT_OBLIQUE='6') 

396 else: 

397 gmtconf = dict( 

398 TICK_PEN='1.25p', 

399 TICK_LENGTH='0.2c', 

400 ANNOT_FONT_PRIMARY='1', 

401 ANNOT_FONT_SIZE_PRIMARY='12p', 

402 LABEL_FONT='1', 

403 LABEL_FONT_SIZE='12p', 

404 CHAR_ENCODING='ISOLatin1+', 

405 BASEMAP_TYPE='fancy', 

406 PLOT_DEGREE_FORMAT='D', 

407 PAPER_MEDIA='Custom_%ix%i' % ( 

408 w*gmtpy.cm, 

409 h*gmtpy.cm), 

410 GRID_PEN_PRIMARY='thinnest/0/50/0', 

411 DOTS_PR_INCH='1200', 

412 OBLIQUE_ANNOTATION='6') 

413 

414 gmtconf.update( 

415 (k.upper(), v) for (k, v) in self.gmt_config.items()) 

416 

417 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion) 

418 

419 layout = gmt.default_layout() 

420 

421 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()]) 

422 

423 widget = layout.get_widget() 

424 widget['P'] = widget['J'] 

425 widget['J'] = ('-JA%g/%g' % (self.lon, self.lat)) + '/%(width)gp' 

426 scaler['R'] = '-R%g/%g/%g/%gr' % self._corners 

427 

428 # aspect = gmtpy.aspect_for_projection( 

429 # gmt.installation['version'], *(widget.J() + scaler.R())) 

430 

431 aspect = self._map_aspect(jr=widget.J() + scaler.R()) 

432 widget.set_aspect(aspect) 

433 

434 self._gmt = gmt 

435 self._layout = layout 

436 self._widget = widget 

437 self._jxyr = self._widget.JXY() + self._scaler.R() 

438 self._pxyr = self._widget.PXY() + [ 

439 '-R%g/%g/%g/%g' % (0, widget.width(), 0, widget.height())] 

440 self._have_drawn_axes = False 

441 self._have_drawn_labels = False 

442 

443 def _draw_background(self): 

444 self._have_topo_land = False 

445 self._have_topo_ocean = False 

446 if self.show_topo: 

447 self._have_topo = self._draw_topo() 

448 

449 self._draw_basefeatures() 

450 

451 def _get_topo_tile(self, k): 

452 t = None 

453 demname = None 

454 for dem in self._dems[k]: 

455 t = topo.get(dem, self._wesn) 

456 demname = dem 

457 if t is not None: 

458 break 

459 

460 if not t: 

461 raise NoTopo() 

462 

463 return t, demname 

464 

465 def _prep_topo(self, k): 

466 gmt = self._gmt 

467 t, demname = self._get_topo_tile(k) 

468 

469 if demname not in self._prep_topo_have: 

470 

471 grdfile = gmt.tempfilename() 

472 

473 is_flat = num.all(t.data[0] == t.data) 

474 

475 gmtpy.savegrd( 

476 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat') 

477 

478 if self.illuminate and not is_flat: 

479 if k == 'ocean': 

480 factor = self.illuminate_factor_ocean 

481 else: 

482 factor = self.illuminate_factor_land 

483 

484 ilumfn = gmt.tempfilename() 

485 gmt.grdgradient( 

486 grdfile, 

487 N='e%g' % factor, 

488 A=-45, 

489 G=ilumfn, 

490 out_discard=True) 

491 

492 ilumargs = ['-I%s' % ilumfn] 

493 else: 

494 ilumargs = [] 

495 

496 if self.replace_topo_color_only: 

497 t2 = self.replace_topo_color_only 

498 grdfile2 = gmt.tempfilename() 

499 

500 gmtpy.savegrd( 

501 t2.x(), t2.y(), t2.data, filename=grdfile2, 

502 naming='lonlat') 

503 

504 if gmt.is_gmt5(): 

505 gmt.grdsample( 

506 grdfile2, 

507 G=grdfile, 

508 n='l', 

509 I='%g/%g' % (t.dx, t.dy), # noqa 

510 R=grdfile, 

511 out_discard=True) 

512 else: 

513 gmt.grdsample( 

514 grdfile2, 

515 G=grdfile, 

516 Q='l', 

517 I='%g/%g' % (t.dx, t.dy), # noqa 

518 R=grdfile, 

519 out_discard=True) 

520 

521 gmt.grdmath( 

522 grdfile, '0.0', 'AND', '=', grdfile2, 

523 out_discard=True) 

524 

525 grdfile = grdfile2 

526 

527 self._prep_topo_have[demname] = grdfile, ilumargs 

528 

529 return self._prep_topo_have[demname] 

530 

531 def _draw_topo(self): 

532 widget = self._widget 

533 scaler = self._scaler 

534 gmt = self._gmt 

535 cres = self._coastline_resolution 

536 minarea = self._minarea 

537 

538 JXY = widget.JXY() 

539 R = scaler.R() 

540 

541 try: 

542 grdfile, ilumargs = self._prep_topo('ocean') 

543 gmt.pscoast(D=cres, S='c', A=minarea, *(JXY+R)) 

544 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_wet), 

545 *(ilumargs+JXY+R)) 

546 gmt.pscoast(Q=True, *(JXY+R)) 

547 self._have_topo_ocean = True 

548 except NoTopo: 

549 self._have_topo_ocean = False 

550 

551 try: 

552 grdfile, ilumargs = self._prep_topo('land') 

553 gmt.pscoast(D=cres, G='c', A=minarea, *(JXY+R)) 

554 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_dry), 

555 *(ilumargs+JXY+R)) 

556 gmt.pscoast(Q=True, *(JXY+R)) 

557 self._have_topo_land = True 

558 except NoTopo: 

559 self._have_topo_land = False 

560 

561 def _draw_topo_scale(self, label='Elevation [km]'): 

562 dry = read_cpt(topo.cpt(self.topo_cpt_dry)) 

563 wet = read_cpt(topo.cpt(self.topo_cpt_wet)) 

564 combi = cpt_merge_wet_dry(wet, dry) 

565 for level in combi.levels: 

566 level.vmin /= km 

567 level.vmax /= km 

568 

569 topo_cpt = self.gmt.tempfilename() + '.cpt' 

570 write_cpt(combi, topo_cpt) 

571 

572 (w, h), (xo, yo) = self.widget.get_size() 

573 self.gmt.psscale( 

574 D='%gp/%gp/%gp/%gph' % (xo + 0.5*w, yo - 2.0*gmtpy.cm, w, 

575 0.5*gmtpy.cm), 

576 C=topo_cpt, 

577 B='1:%s:' % label) 

578 

579 def _draw_basefeatures(self): 

580 gmt = self._gmt 

581 cres = self._coastline_resolution 

582 rivers = self._rivers 

583 minarea = self._minarea 

584 

585 color_wet = self.color_wet 

586 color_dry = self.color_dry 

587 

588 if self.show_rivers and rivers: 

589 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)] 

590 else: 

591 rivers = [] 

592 

593 fill = {} 

594 if not self._have_topo_land: 

595 fill['G'] = color_dry 

596 

597 if not self._have_topo_ocean: 

598 fill['S'] = color_wet 

599 

600 if self.show_boundaries: 

601 fill['N'] = '1/1p,%s,%s' % ( 

602 gmtpy.color(self.color_boundaries), 'solid') 

603 

604 gmt.pscoast( 

605 D=cres, 

606 W='thinnest,%s' % gmtpy.color(darken(gmtpy.color_tup(color_dry))), 

607 A=minarea, 

608 *(rivers+self._jxyr), **fill) 

609 

610 if self.show_plates: 

611 self.draw_plates() 

612 

613 def _draw_axes(self): 

614 gmt = self._gmt 

615 scaler = self._scaler 

616 widget = self._widget 

617 

618 if self.axes_layout is None: 

619 if self.lat > 0.0: 

620 axes_layout = 'WSen' 

621 else: 

622 axes_layout = 'WseN' 

623 else: 

624 axes_layout = self.axes_layout 

625 

626 scale_km = nice_value(self.radius/5.) / 1000. 

627 

628 if self.show_center_mark: 

629 gmt.psxy( 

630 in_rows=[[self.lon, self.lat]], 

631 S='c20p', W='2p,black', 

632 *self._jxyr) 

633 

634 if self.show_grid: 

635 btmpl = ('%(xinc)gg%(xinc)g:%(xlabel)s:/' 

636 '%(yinc)gg%(yinc)g:%(ylabel)s:') 

637 else: 

638 btmpl = '%(xinc)g:%(xlabel)s:/%(yinc)g:%(ylabel)s:' 

639 

640 if self.show_scale: 

641 scale = 'x%gp/%gp/%g/%g/%gk' % ( 

642 6./7*widget.width(), 

643 widget.height()/7., 

644 self.lon, 

645 self.lat, 

646 scale_km) 

647 else: 

648 scale = False 

649 

650 gmt.psbasemap( 

651 B=(btmpl % scaler.get_params())+axes_layout, 

652 L=scale, 

653 *self._jxyr) 

654 

655 if self.comment: 

656 font_size = self.gmt.label_font_size() 

657 

658 _, east, south, _ = self._wesn 

659 if gmt.is_gmt5(): 

660 row = [ 

661 1, 0, 

662 '%gp,%s,%s' % (font_size, 0, 'black'), 'BR', 

663 self.comment] 

664 

665 farg = ['-F+f+j'] 

666 else: 

667 row = [1, 0, font_size, 0, 0, 'BR', self.comment] 

668 farg = [] 

669 

670 gmt.pstext( 

671 in_rows=[row], 

672 N=True, 

673 R=(0, 1, 0, 1), 

674 D='%gp/%gp' % (-font_size*0.2, font_size*0.3), 

675 *(widget.PXY() + farg)) 

676 

677 def draw_axes(self): 

678 if not self._have_drawn_axes: 

679 self._draw_axes() 

680 self._have_drawn_axes = True 

681 

682 def _have_coastlines(self): 

683 gmt = self._gmt 

684 cres = self._coastline_resolution 

685 minarea = self._minarea 

686 

687 checkfile = gmt.tempfilename() 

688 

689 gmt.pscoast( 

690 M=True, 

691 D=cres, 

692 W='thinnest,black', 

693 A=minarea, 

694 out_filename=checkfile, 

695 *self._jxyr) 

696 

697 points = [] 

698 with open(checkfile, 'r') as f: 

699 for line in f: 

700 ls = line.strip() 

701 if ls.startswith('#') or ls.startswith('>') or ls == '': 

702 continue 

703 plon, plat = [float(x) for x in ls.split()] 

704 points.append((plat, plon)) 

705 

706 points = num.array(points, dtype=float) 

707 return num.any(points_in_region(points, self._wesn)) 

708 

709 def have_coastlines(self): 

710 self.gmt 

711 return self._have_coastlines() 

712 

713 def project(self, lats, lons, jr=None): 

714 onepoint = False 

715 if isinstance(lats, float) and isinstance(lons, float): 

716 lats = [lats] 

717 lons = [lons] 

718 onepoint = True 

719 

720 if jr is not None: 

721 j, r = jr 

722 gmt = gmtpy.GMT(version=self._gmtversion) 

723 else: 

724 j, _, _, r = self.jxyr 

725 gmt = self.gmt 

726 

727 f = BytesIO() 

728 gmt.mapproject(j, r, in_columns=(lons, lats), out_stream=f, D='p') 

729 f.seek(0) 

730 data = num.loadtxt(f, ndmin=2) 

731 xs, ys = data.T 

732 if onepoint: 

733 xs = xs[0] 

734 ys = ys[0] 

735 return xs, ys 

736 

737 def _map_box(self, jr=None): 

738 ll_lon, ll_lat, ur_lon, ur_lat = self._corners 

739 

740 xs_corner, ys_corner = self.project( 

741 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr) 

742 

743 w = xs_corner[1] - xs_corner[0] 

744 h = ys_corner[1] - ys_corner[0] 

745 

746 return w, h 

747 

748 def _map_aspect(self, jr=None): 

749 w, h = self._map_box(jr=jr) 

750 return h/w 

751 

752 def _draw_labels(self): 

753 points_taken = [] 

754 regions_taken = [] 

755 

756 def no_points_in_rect(xs, ys, xmin, ymin, xmax, ymax): 

757 xx = not num.any(la(la(xmin < xs, xs < xmax), 

758 la(ymin < ys, ys < ymax))) 

759 return xx 

760 

761 def roverlaps(a, b): 

762 return (a[0] < b[2] and b[0] < a[2] and 

763 a[1] < b[3] and b[1] < a[3]) 

764 

765 w, h = self._map_box() 

766 

767 label_font_size = self.gmt.label_font_size() 

768 

769 if self._labels: 

770 

771 n = len(self._labels) 

772 

773 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \ 

774 angles, styles = list(zip(*self._labels)) 

775 

776 font_sizes = [ 

777 (font_size or label_font_size) for font_size in font_sizes] 

778 

779 sx = num.array(sx, dtype=float) 

780 sy = num.array(sy, dtype=float) 

781 

782 xs, ys = self.project(lats, lons) 

783 

784 points_taken.append((xs, ys)) 

785 

786 dxs = num.zeros(n) 

787 dys = num.zeros(n) 

788 

789 for i in range(n): 

790 dx, dy = gmtpy.text_box( 

791 texts[i], 

792 font=fonts[i], 

793 font_size=font_sizes[i], 

794 **styles[i]) 

795 

796 dxs[i] = dx 

797 dys[i] = dy 

798 

799 la = num.logical_and 

800 anchors_ok = ( 

801 la(xs + sx + dxs < w, ys + sy + dys < h), 

802 la(xs - sx - dxs > 0., ys - sy - dys > 0.), 

803 la(xs + sx + dxs < w, ys - sy - dys > 0.), 

804 la(xs - sx - dxs > 0., ys + sy + dys < h), 

805 ) 

806 

807 arects = [ 

808 (xs, ys, xs + sx + dxs, ys + sy + dys), 

809 (xs - sx - dxs, ys - sy - dys, xs, ys), 

810 (xs, ys - sy - dys, xs + sx + dxs, ys), 

811 (xs - sx - dxs, ys, xs, ys + sy + dys)] 

812 

813 for i in range(n): 

814 for ianch in range(4): 

815 anchors_ok[ianch][i] &= no_points_in_rect( 

816 xs, ys, *[xxx[i] for xxx in arects[ianch]]) 

817 

818 anchor_choices = [] 

819 anchor_take = [] 

820 for i in range(n): 

821 choices = [ianch for ianch in range(4) 

822 if anchors_ok[ianch][i]] 

823 anchor_choices.append(choices) 

824 if choices: 

825 anchor_take.append(choices[0]) 

826 else: 

827 anchor_take.append(None) 

828 

829 def cost(anchor_take): 

830 noverlaps = 0 

831 for i in range(n): 

832 for j in range(n): 

833 if i != j: 

834 i_take = anchor_take[i] 

835 j_take = anchor_take[j] 

836 if i_take is None or j_take is None: 

837 continue 

838 r_i = [xxx[i] for xxx in arects[i_take]] 

839 r_j = [xxx[j] for xxx in arects[j_take]] 

840 if roverlaps(r_i, r_j): 

841 noverlaps += 1 

842 

843 return noverlaps 

844 

845 cur_cost = cost(anchor_take) 

846 imax = 30 

847 while cur_cost != 0 and imax > 0: 

848 for i in range(n): 

849 for t in anchor_choices[i]: 

850 anchor_take_new = list(anchor_take) 

851 anchor_take_new[i] = t 

852 new_cost = cost(anchor_take_new) 

853 if new_cost < cur_cost: 

854 anchor_take = anchor_take_new 

855 cur_cost = new_cost 

856 

857 imax -= 1 

858 

859 while cur_cost != 0: 

860 for i in range(n): 

861 anchor_take_new = list(anchor_take) 

862 anchor_take_new[i] = None 

863 new_cost = cost(anchor_take_new) 

864 if new_cost < cur_cost: 

865 anchor_take = anchor_take_new 

866 cur_cost = new_cost 

867 break 

868 

869 anchor_strs = ['BL', 'TR', 'TL', 'BR'] 

870 

871 for i in range(n): 

872 ianchor = anchor_take[i] 

873 color = colors[i] 

874 if color is None: 

875 color = 'black' 

876 

877 if ianchor is not None: 

878 regions_taken.append([xxx[i] for xxx in arects[ianchor]]) 

879 

880 anchor = anchor_strs[ianchor] 

881 

882 yoff = [-sy[i], sy[i]][anchor[0] == 'B'] 

883 xoff = [-sx[i], sx[i]][anchor[1] == 'L'] 

884 if self.gmt.is_gmt5(): 

885 row = ( 

886 lons[i], lats[i], 

887 '%i,%s,%s' % (font_sizes[i], fonts[i], color), 

888 anchor, 

889 texts[i]) 

890 

891 farg = ['-F+f+j+a%g' % angles[i]] 

892 else: 

893 row = ( 

894 lons[i], lats[i], 

895 font_sizes[i], angles[i], fonts[i], anchor, 

896 texts[i]) 

897 farg = ['-G%s' % color] 

898 

899 self.gmt.pstext( 

900 in_rows=[row], 

901 D='%gp/%gp' % (xoff, yoff), 

902 *(self.jxyr + farg), 

903 **styles[i]) 

904 

905 if self._area_labels: 

906 

907 for lons, lats, text, color, font, font_size, style in \ 

908 self._area_labels: 

909 

910 if font_size is None: 

911 font_size = label_font_size 

912 

913 if color is None: 

914 color = 'black' 

915 

916 if self.gmt.is_gmt5(): 

917 farg = ['-F+f+j'] 

918 else: 

919 farg = ['-G%s' % color] 

920 

921 xs, ys = self.project(lats, lons) 

922 dx, dy = gmtpy.text_box( 

923 text, font=font, font_size=font_size, **style) 

924 

925 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy] 

926 

927 locs_ok = num.ones(xs.size, dtype=bool) 

928 

929 for iloc in range(xs.size): 

930 rcandi = [xxx[iloc] for xxx in rects] 

931 

932 locs_ok[iloc] = True 

933 locs_ok[iloc] &= ( 

934 0 < rcandi[0] and rcandi[2] < w 

935 and 0 < rcandi[1] and rcandi[3] < h) 

936 

937 overlap = False 

938 for r in regions_taken: 

939 if roverlaps(r, rcandi): 

940 overlap = True 

941 break 

942 

943 locs_ok[iloc] &= not overlap 

944 

945 for xs_taken, ys_taken in points_taken: 

946 locs_ok[iloc] &= no_points_in_rect( 

947 xs_taken, ys_taken, *rcandi) 

948 

949 if not locs_ok[iloc]: 

950 break 

951 

952 rows = [] 

953 for iloc, (lon, lat) in enumerate(zip(lons, lats)): 

954 if not locs_ok[iloc]: 

955 continue 

956 

957 if self.gmt.is_gmt5(): 

958 row = ( 

959 lon, lat, 

960 '%i,%s,%s' % (font_size, font, color), 

961 'MC', 

962 text) 

963 

964 else: 

965 row = ( 

966 lon, lat, 

967 font_size, 0, font, 'MC', 

968 text) 

969 

970 rows.append(row) 

971 

972 regions_taken.append([xxx[iloc] for xxx in rects]) 

973 break 

974 

975 self.gmt.pstext( 

976 in_rows=rows, 

977 *(self.jxyr + farg), 

978 **style) 

979 

980 def draw_labels(self): 

981 self.gmt 

982 if not self._have_drawn_labels: 

983 self._draw_labels() 

984 self._have_drawn_labels = True 

985 

986 def add_label( 

987 self, lat, lon, text, 

988 offset_x=5., offset_y=5., 

989 color=None, 

990 font='1', 

991 font_size=None, 

992 angle=0, 

993 style={}): 

994 

995 if 'G' in style: 

996 style = style.copy() 

997 color = style.pop('G') 

998 

999 self._labels.append( 

1000 (lon, lat, text, offset_x, offset_y, color, font, font_size, 

1001 angle, style)) 

1002 

1003 def add_area_label( 

1004 self, lat, lon, text, 

1005 color=None, 

1006 font='3', 

1007 font_size=None, 

1008 style={}): 

1009 

1010 self._area_labels.append( 

1011 (lon, lat, text, color, font, font_size, style)) 

1012 

1013 def cities_in_region(self): 

1014 from pyrocko.dataset import geonames 

1015 cities = geonames.get_cities_region(region=self.wesn, minpop=0) 

1016 cities.extend(self.custom_cities) 

1017 cities.sort(key=lambda x: x.population) 

1018 return cities 

1019 

1020 def draw_cities(self, 

1021 exact=None, 

1022 include=[], 

1023 exclude=[], 

1024 nmax_soft=10, 

1025 psxy_style=dict(S='s5p', G='black')): 

1026 

1027 cities = self.cities_in_region() 

1028 

1029 if exact is not None: 

1030 cities = [c for c in cities if c.name in exact] 

1031 minpop = None 

1032 else: 

1033 cities = [c for c in cities if c.name not in exclude] 

1034 minpop = 10**3 

1035 for minpop_new in [1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6, 3e6, 1e7]: 

1036 cities_new = [ 

1037 c for c in cities 

1038 if c.population > minpop_new or c.name in include] 

1039 

1040 if len(cities_new) == 0 or ( 

1041 len(cities_new) < 3 and len(cities) < nmax_soft*2): 

1042 break 

1043 

1044 cities = cities_new 

1045 minpop = minpop_new 

1046 if len(cities) <= nmax_soft: 

1047 break 

1048 

1049 if cities: 

1050 lats = [c.lat for c in cities] 

1051 lons = [c.lon for c in cities] 

1052 

1053 self.gmt.psxy( 

1054 in_columns=(lons, lats), 

1055 *self.jxyr, **psxy_style) 

1056 

1057 for c in cities: 

1058 try: 

1059 text = c.name.encode('iso-8859-1').decode('iso-8859-1') 

1060 except UnicodeEncodeError: 

1061 text = c.asciiname 

1062 

1063 self.add_label(c.lat, c.lon, text) 

1064 

1065 self._cities_minpop = minpop 

1066 

1067 def add_stations(self, stations, psxy_style=dict()): 

1068 

1069 default_psxy_style = { 

1070 'S': 't8p', 

1071 'G': 'black' 

1072 } 

1073 default_psxy_style.update(psxy_style) 

1074 

1075 lats, lons = zip(*[s.effective_latlon for s in stations]) 

1076 

1077 self.gmt.psxy( 

1078 in_columns=(lons, lats), 

1079 *self.jxyr, **default_psxy_style) 

1080 

1081 for station in stations: 

1082 self.add_label( 

1083 station.effective_lat, 

1084 station.effective_lon, 

1085 '.'.join(x for x in (station.network, station.station) if x)) 

1086 

1087 def add_kite_scene(self, scene): 

1088 tile = FloatTile( 

1089 scene.frame.llLon, 

1090 scene.frame.llLat, 

1091 scene.frame.dLon, 

1092 scene.frame.dLat, 

1093 scene.displacement) 

1094 

1095 return tile 

1096 

1097 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None, 

1098 labels=True, vertical=False, fontsize=10): 

1099 

1100 stations = campaign.stations 

1101 

1102 if offset_scale is None: 

1103 offset_scale = num.zeros(campaign.nstations) 

1104 for ista, sta in enumerate(stations): 

1105 for comp in sta.components.values(): 

1106 offset_scale[ista] += comp.shift 

1107 offset_scale = num.sqrt(offset_scale**2).max() 

1108 

1109 size = math.sqrt(self.height**2 + self.width**2) 

1110 scale = (size/10.) / offset_scale 

1111 logger.debug('GNSS: Using offset scale %f, map scale %f', 

1112 offset_scale, scale) 

1113 

1114 lats, lons = zip(*[s.effective_latlon for s in stations]) 

1115 

1116 if self.gmt.is_gmt6(): 

1117 sign_factor = 1. 

1118 arrow_head_placement = 'e' 

1119 else: 

1120 sign_factor = -1. 

1121 arrow_head_placement = 'b' 

1122 

1123 if vertical: 

1124 rows = [[lons[ista], lats[ista], 

1125 0., sign_factor * s.up.shift, 

1126 (s.east.sigma + s.north.sigma) if s.east.sigma else 0., 

1127 s.up.sigma, 0., 

1128 s.code if labels else None] 

1129 for ista, s in enumerate(stations) 

1130 if s.up is not None] 

1131 

1132 else: 

1133 rows = [[lons[ista], lats[ista], 

1134 sign_factor * s.east.shift, sign_factor * s.north.shift, 

1135 s.east.sigma, s.north.sigma, s.correlation_ne, 

1136 s.code if labels else None] 

1137 for ista, s in enumerate(stations) 

1138 if s.east is not None or s.north is not None] 

1139 

1140 default_psxy_style = { 

1141 'h': 0, 

1142 'W': '2p,black', 

1143 'A': '+p2p,black+{}+a40'.format(arrow_head_placement), 

1144 'G': 'black', 

1145 'L': True, 

1146 'S': 'e%dc/0.95/%d' % (scale, fontsize), 

1147 } 

1148 

1149 if not labels: 

1150 for row in rows: 

1151 row.pop(-1) 

1152 

1153 if psxy_style is not None: 

1154 default_psxy_style.update(psxy_style) 

1155 

1156 self.gmt.psvelo( 

1157 in_rows=rows, 

1158 *self.jxyr, 

1159 **default_psxy_style) 

1160 

1161 def draw_plates(self): 

1162 from pyrocko.dataset import tectonics 

1163 

1164 neast = 20 

1165 nnorth = max(1, int(round(num.round(self._hreg/self._wreg * neast)))) 

1166 norths = num.linspace(-self._hreg*0.5, self._hreg*0.5, nnorth) 

1167 easts = num.linspace(-self._wreg*0.5, self._wreg*0.5, neast) 

1168 norths2 = num.repeat(norths, neast) 

1169 easts2 = num.tile(easts, nnorth) 

1170 lats, lons = od.ne_to_latlon( 

1171 self.lat, self.lon, norths2, easts2) 

1172 

1173 bird = tectonics.PeterBird2003() 

1174 plates = bird.get_plates() 

1175 

1176 color_plates = gmtpy.color('aluminium5') 

1177 color_velocities = gmtpy.color('skyblue1') 

1178 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1'))) 

1179 

1180 points = num.vstack((lats, lons)).T 

1181 used = [] 

1182 for plate in plates: 

1183 mask = plate.contains_points(points) 

1184 if num.any(mask): 

1185 used.append((plate, mask)) 

1186 

1187 if len(used) > 1: 

1188 

1189 candi_fixed = {} 

1190 

1191 label_data = [] 

1192 for plate, mask in used: 

1193 

1194 mean_north = num.mean(norths2[mask]) 

1195 mean_east = num.mean(easts2[mask]) 

1196 iorder = num.argsort(num.sqrt( 

1197 (norths2[mask] - mean_north)**2 + 

1198 (easts2[mask] - mean_east)**2)) 

1199 

1200 lat_candis = lats[mask][iorder] 

1201 lon_candis = lons[mask][iorder] 

1202 

1203 candi_fixed[plate.name] = lat_candis.size 

1204 

1205 label_data.append(( 

1206 lat_candis, lon_candis, plate, color_plates)) 

1207 

1208 boundaries = bird.get_boundaries() 

1209 

1210 size = 1. 

1211 

1212 psxy_kwargs = [] 

1213 

1214 for boundary in boundaries: 

1215 if num.any(points_in_region(boundary.points, self._wesn)): 

1216 for typ, part in boundary.split_types( 

1217 [['SUB'], 

1218 ['OSR', 'CRB'], 

1219 ['OTF', 'CTF', 'OCB', 'CCB']]): 

1220 

1221 lats, lons = part.T 

1222 

1223 kwargs = {} 

1224 

1225 kwargs['in_columns'] = (lons, lats) 

1226 kwargs['W'] = '%gp,%s' % (size, color_plates) 

1227 

1228 if typ[0] == 'SUB': 

1229 if boundary.kind == '\\': 

1230 kwargs['S'] = 'f%g/%gp+t+r' % ( 

1231 0.45*size, 3.*size) 

1232 elif boundary.kind == '/': 

1233 kwargs['S'] = 'f%g/%gp+t+l' % ( 

1234 0.45*size, 3.*size) 

1235 

1236 kwargs['G'] = color_plates 

1237 

1238 elif typ[0] in ['OSR', 'CRB']: 

1239 kwargs_bg = {} 

1240 kwargs_bg['in_columns'] = (lons, lats) 

1241 kwargs_bg['W'] = '%gp,%s' % ( 

1242 size * 3, color_plates) 

1243 psxy_kwargs.append(kwargs_bg) 

1244 

1245 kwargs['W'] = '%gp,%s' % (size * 2, 'white') 

1246 

1247 psxy_kwargs.append(kwargs) 

1248 

1249 if boundary.kind == '\\': 

1250 if boundary.plate_name2 in candi_fixed: 

1251 candi_fixed[boundary.plate_name2] += \ 

1252 neast*nnorth 

1253 

1254 elif boundary.kind == '/': 

1255 if boundary.plate_name1 in candi_fixed: 

1256 candi_fixed[boundary.plate_name1] += \ 

1257 neast*nnorth 

1258 

1259 candi_fixed = [name for name in sorted( 

1260 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])] 

1261 

1262 candi_fixed.append(None) 

1263 

1264 gsrm = tectonics.GSRM1() 

1265 

1266 for name in candi_fixed: 

1267 if name not in gsrm.plate_names() \ 

1268 and name not in gsrm.plate_alt_names(): 

1269 

1270 continue 

1271 

1272 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \ 

1273 gsrm.get_velocities(name, region=self._wesn) 

1274 

1275 fixed_plate_name = name 

1276 

1277 if self.show_plate_velocities: 

1278 self.gmt.psvelo( 

1279 in_columns=( 

1280 lons, lats, veast, vnorth, veast_err, vnorth_err, 

1281 corr), 

1282 W='0.25p,%s' % color_velocities, 

1283 A='9p+e+g%s' % color_velocities, 

1284 S='e0.2p/0.95/10', 

1285 *self.jxyr) 

1286 

1287 for _ in range(len(lons) // 50 + 1): 

1288 ii = random.randint(0, len(lons)-1) 

1289 v = math.sqrt(vnorth[ii]**2 + veast[ii]**2) 

1290 self.add_label( 

1291 lats[ii], lons[ii], '%.0f' % v, 

1292 font_size=0.7*self.gmt.label_font_size(), 

1293 style=dict( 

1294 G=color_velocities_lab)) 

1295 

1296 break 

1297 

1298 if self.show_plate_names: 

1299 for (lat_candis, lon_candis, plate, color) in label_data: 

1300 full_name = bird.full_name(plate.name) 

1301 if plate.name == fixed_plate_name: 

1302 full_name = '@_' + full_name + '@_' 

1303 

1304 self.add_area_label( 

1305 lat_candis, lon_candis, 

1306 full_name, 

1307 color=color, 

1308 font='3') 

1309 

1310 for kwargs in psxy_kwargs: 

1311 self.gmt.psxy(*self.jxyr, **kwargs) 

1312 

1313 

1314def rand(mi, ma): 

1315 mi = float(mi) 

1316 ma = float(ma) 

1317 return random.random() * (ma-mi) + mi 

1318 

1319 

1320def split_region(region): 

1321 west, east, south, north = topo.positive_region(region) 

1322 if east > 180: 

1323 return [(west, 180., south, north), 

1324 (-180., east-360., south, north)] 

1325 else: 

1326 return [region] 

1327 

1328 

1329class CPTLevel(Object): 

1330 vmin = Float.T() 

1331 vmax = Float.T() 

1332 color_min = Tuple.T(3, Float.T()) 

1333 color_max = Tuple.T(3, Float.T()) 

1334 

1335 

1336class CPT(Object): 

1337 color_below = Tuple.T(3, Float.T(), optional=True) 

1338 color_above = Tuple.T(3, Float.T(), optional=True) 

1339 color_nan = Tuple.T(3, Float.T(), optional=True) 

1340 levels = List.T(CPTLevel.T()) 

1341 

1342 def scale(self, vmin, vmax): 

1343 vmin_old, vmax_old = self.levels[0].vmin, self.levels[-1].vmax 

1344 for level in self.levels: 

1345 level.vmin = (level.vmin - vmin_old) / (vmax_old - vmin_old) * \ 

1346 (vmax - vmin) + vmin 

1347 level.vmax = (level.vmax - vmin_old) / (vmax_old - vmin_old) * \ 

1348 (vmax - vmin) + vmin 

1349 

1350 def discretize(self, nlevels): 

1351 colors = [] 

1352 vals = [] 

1353 for level in self.levels: 

1354 vals.append(level.vmin) 

1355 vals.append(level.vmax) 

1356 colors.append(level.color_min) 

1357 colors.append(level.color_max) 

1358 

1359 r, g, b = num.array(colors, dtype=float).T 

1360 vals = num.array(vals, dtype=float) 

1361 

1362 vmin, vmax = self.levels[0].vmin, self.levels[-1].vmax 

1363 x = num.linspace(vmin, vmax, nlevels+1) 

1364 rd = num.interp(x, vals, r) 

1365 gd = num.interp(x, vals, g) 

1366 bd = num.interp(x, vals, b) 

1367 

1368 levels = [] 

1369 for ilevel in range(nlevels): 

1370 color = ( 

1371 float(0.5*(rd[ilevel]+rd[ilevel+1])), 

1372 float(0.5*(gd[ilevel]+gd[ilevel+1])), 

1373 float(0.5*(bd[ilevel]+bd[ilevel+1]))) 

1374 

1375 levels.append(CPTLevel( 

1376 vmin=x[ilevel], 

1377 vmax=x[ilevel+1], 

1378 color_min=color, 

1379 color_max=color)) 

1380 

1381 cpt = CPT( 

1382 color_below=self.color_below, 

1383 color_above=self.color_above, 

1384 color_nan=self.color_nan, 

1385 levels=levels) 

1386 

1387 return cpt 

1388 

1389 

1390class CPTParseError(Exception): 

1391 pass 

1392 

1393 

1394def read_cpt(filename): 

1395 with open(filename) as f: 

1396 color_below = None 

1397 color_above = None 

1398 color_nan = None 

1399 levels = [] 

1400 try: 

1401 for line in f: 

1402 line = line.strip() 

1403 toks = line.split() 

1404 

1405 if line.startswith('#'): 

1406 continue 

1407 

1408 elif line.startswith('B'): 

1409 color_below = tuple(map(float, toks[1:4])) 

1410 

1411 elif line.startswith('F'): 

1412 color_above = tuple(map(float, toks[1:4])) 

1413 

1414 elif line.startswith('N'): 

1415 color_nan = tuple(map(float, toks[1:4])) 

1416 

1417 else: 

1418 values = list(map(float, line.split())) 

1419 vmin = values[0] 

1420 color_min = tuple(values[1:4]) 

1421 vmax = values[4] 

1422 color_max = tuple(values[5:8]) 

1423 levels.append(CPTLevel( 

1424 vmin=vmin, 

1425 vmax=vmax, 

1426 color_min=color_min, 

1427 color_max=color_max)) 

1428 

1429 except Exception: 

1430 raise CPTParseError() 

1431 

1432 return CPT( 

1433 color_below=color_below, 

1434 color_above=color_above, 

1435 color_nan=color_nan, 

1436 levels=levels) 

1437 

1438 

1439def color_to_int(color): 

1440 return tuple(max(0, min(255, int(round(x)))) for x in color) 

1441 

1442 

1443def write_cpt(cpt, filename): 

1444 with open(filename, 'w') as f: 

1445 for level in cpt.levels: 

1446 f.write( 

1447 '%e %i %i %i %e %i %i %i\n' % 

1448 ((level.vmin, ) + color_to_int(level.color_min) + 

1449 (level.vmax, ) + color_to_int(level.color_max))) 

1450 

1451 if cpt.color_below: 

1452 f.write('B %i %i %i\n' % color_to_int(cpt.color_below)) 

1453 

1454 if cpt.color_above: 

1455 f.write('F %i %i %i\n' % color_to_int(cpt.color_above)) 

1456 

1457 if cpt.color_nan: 

1458 f.write('N %i %i %i\n' % color_to_int(cpt.color_nan)) 

1459 

1460 

1461def cpt_merge_wet_dry(wet, dry): 

1462 levels = [] 

1463 for level in wet.levels: 

1464 if level.vmin < 0.: 

1465 if level.vmax > 0.: 

1466 level.vmax = 0. 

1467 

1468 levels.append(level) 

1469 

1470 for level in dry.levels: 

1471 if level.vmax > 0.: 

1472 if level.vmin < 0.: 

1473 level.vmin = 0. 

1474 

1475 levels.append(level) 

1476 

1477 combi = CPT( 

1478 color_below=wet.color_below, 

1479 color_above=dry.color_above, 

1480 color_nan=dry.color_nan, 

1481 levels=levels) 

1482 

1483 return combi 

1484 

1485 

1486if __name__ == '__main__': 

1487 from pyrocko import util 

1488 util.setup_logging('pyrocko.automap', 'info') 

1489 

1490 import sys 

1491 if len(sys.argv) == 2: 

1492 

1493 n = int(sys.argv[1]) 

1494 

1495 for i in range(n): 

1496 m = Map( 

1497 lat=rand(-60., 60.), 

1498 lon=rand(-180., 180.), 

1499 radius=math.exp(rand(math.log(500*km), math.log(3000*km))), 

1500 width=30., height=30., 

1501 show_grid=True, 

1502 show_topo=True, 

1503 color_dry=(238, 236, 230), 

1504 topo_cpt_wet='light_sea_uniform', 

1505 topo_cpt_dry='light_land_uniform', 

1506 illuminate=True, 

1507 illuminate_factor_ocean=0.15, 

1508 show_rivers=False, 

1509 show_plates=True) 

1510 

1511 m.draw_cities() 

1512 print(m) 

1513 m.save('map_%02i.pdf' % i)