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 

24 

25try: 

26 newstr = unicode 

27except NameError: 

28 newstr = str 

29 

30points_in_region = od.points_in_region 

31 

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

33 

34earthradius = 6371000.0 

35r2d = 180./math.pi 

36d2r = 1./r2d 

37km = 1000. 

38d2m = d2r*earthradius 

39m2d = 1./d2m 

40cm = gmtpy.cm 

41 

42 

43def darken(c, f=0.7): 

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

45 

46 

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

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

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

50 return ll_lon, ll_lat, ur_lon, ur_lat 

51 

52 

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

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

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

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

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

58 south = slats.min() 

59 north = nlats.max() 

60 

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

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

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

64 

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

66 west = -180. 

67 east = 180. 

68 else: 

69 west = wlons.min() 

70 east = elons.max() 

71 

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

73 

74 

75class NoTopo(Exception): 

76 pass 

77 

78 

79class OutOfBounds(Exception): 

80 pass 

81 

82 

83class FloatTile(Object): 

84 xmin = Float.T() 

85 ymin = Float.T() 

86 dx = Float.T() 

87 dy = Float.T() 

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

89 

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

91 Object.__init__(self, init_props=False) 

92 self.xmin = float(xmin) 

93 self.ymin = float(ymin) 

94 self.dx = float(dx) 

95 self.dy = float(dy) 

96 self.data = data 

97 self._set_maxes() 

98 

99 def _set_maxes(self): 

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

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

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

103 

104 def x(self): 

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

106 

107 def y(self): 

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

109 

110 def get(self, x, y): 

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

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

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

114 return self.data[iy, ix] 

115 else: 

116 raise OutOfBounds() 

117 

118 

119class City(Object): 

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

121 name = newstr(name) 

122 lat = float(lat) 

123 lon = float(lon) 

124 if asciiname is None: 

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

126 

127 if population is None: 

128 population = 0 

129 else: 

130 population = int(population) 

131 

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

133 population=population, asciiname=asciiname) 

134 

135 name = Unicode.T() 

136 lat = Float.T() 

137 lon = Float.T() 

138 population = Int.T() 

139 asciiname = String.T() 

140 

141 

142class Map(Object): 

143 lat = Float.T(optional=True) 

144 lon = Float.T(optional=True) 

145 radius = Float.T(optional=True) 

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

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

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

149 illuminate = Bool.T(default=True) 

150 skip_feature_factor = Float.T(default=0.02) 

151 show_grid = Bool.T(default=False) 

152 show_topo = Bool.T(default=True) 

153 show_scale = Bool.T(default=False) 

154 show_topo_scale = Bool.T(default=False) 

155 show_center_mark = Bool.T(default=False) 

156 show_rivers = Bool.T(default=True) 

157 show_plates = Bool.T(default=False) 

158 show_plate_velocities = Bool.T(default=False) 

159 show_plate_names = Bool.T(default=False) 

160 show_boundaries = Bool.T(default=False) 

161 illuminate_factor_land = Float.T(default=0.5) 

162 illuminate_factor_ocean = Float.T(default=0.25) 

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

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

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

166 topo_resolution_min = Float.T( 

167 default=40., 

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

169 topo_resolution_max = Float.T( 

170 default=200., 

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

172 replace_topo_color_only = FloatTile.T( 

173 optional=True, 

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

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

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

177 axes_layout = String.T(optional=True) 

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

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

180 comment = String.T(optional=True) 

181 approx_ticks = Int.T(default=4) 

182 

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

184 Object.__init__(self, **kwargs) 

185 self._gmt = None 

186 self._scaler = None 

187 self._widget = None 

188 self._corners = None 

189 self._wesn = None 

190 self._minarea = None 

191 self._coastline_resolution = None 

192 self._rivers = None 

193 self._dems = None 

194 self._have_topo_land = None 

195 self._have_topo_ocean = None 

196 self._jxyr = None 

197 self._prep_topo_have = None 

198 self._labels = [] 

199 self._area_labels = [] 

200 self._gmtversion = gmtversion 

201 

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

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

204 

205 ''' 

206 Save the image. 

207 

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

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

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

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

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

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

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

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

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

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

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

219 square with given side length. To save transparency use 

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

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

222 ''' 

223 

224 gmt = self.gmt 

225 self.draw_labels() 

226 self.draw_axes() 

227 if self.show_topo and self.show_topo_scale: 

228 self._draw_topo_scale() 

229 

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

231 crop_eps_mode=crop_eps_mode, 

232 size=size, width=width, height=height, psconvert=psconvert) 

233 

234 @property 

235 def scaler(self): 

236 if self._scaler is None: 

237 self._setup_geometry() 

238 

239 return self._scaler 

240 

241 @property 

242 def wesn(self): 

243 if self._wesn is None: 

244 self._setup_geometry() 

245 

246 return self._wesn 

247 

248 @property 

249 def widget(self): 

250 if self._widget is None: 

251 self._setup() 

252 

253 return self._widget 

254 

255 @property 

256 def layout(self): 

257 if self._layout is None: 

258 self._setup() 

259 

260 return self._layout 

261 

262 @property 

263 def jxyr(self): 

264 if self._jxyr is None: 

265 self._setup() 

266 

267 return self._jxyr 

268 

269 @property 

270 def pxyr(self): 

271 if self._pxyr is None: 

272 self._setup() 

273 

274 return self._pxyr 

275 

276 @property 

277 def gmt(self): 

278 if self._gmt is None: 

279 self._setup() 

280 

281 if self._have_topo_ocean is None: 

282 self._draw_background() 

283 

284 return self._gmt 

285 

286 def _setup(self): 

287 if not self._widget: 

288 self._setup_geometry() 

289 

290 self._setup_lod() 

291 self._setup_gmt() 

292 

293 def _setup_geometry(self): 

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

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

296 wpage -= ml + mr 

297 hpage -= mt + mb 

298 

299 wreg = self.radius * 2.0 

300 hreg = self.radius * 2.0 

301 if wpage >= hpage: 

302 wreg *= wpage/hpage 

303 else: 

304 hreg *= hpage/wpage 

305 

306 self._wreg = wreg 

307 self._hreg = hreg 

308 

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

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

311 

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

313 

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

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

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

317 scaled_unit='km', scaled_unit_factor=0.001) 

318 

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

320 

321 par = scaler.get_params() 

322 

323 west = par['xmin'] 

324 east = par['xmax'] 

325 south = par['ymin'] 

326 north = par['ymax'] 

327 

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

329 self._scaler = scaler 

330 

331 def _setup_lod(self): 

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

333 if self.radius > 1500.*km: 

334 coastline_resolution = 'i' 

335 rivers = False 

336 else: 

337 coastline_resolution = 'f' 

338 rivers = True 

339 

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

341 

342 self._coastline_resolution = coastline_resolution 

343 self._rivers = rivers 

344 

345 self._prep_topo_have = {} 

346 self._dems = {} 

347 

348 cm2inch = gmtpy.cm/gmtpy.inch 

349 

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

351 (self.height * cm2inch)) 

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

353 (self.height * cm2inch)) 

354 

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

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

357 if self._dems[k]: 

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

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

360 

361 def _expand_margins(self): 

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

363 ml = mr = mt = mb = 2.0 

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

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

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

367 ml = mr = self.margins[0] 

368 mt = mb = self.margins[1] 

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

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

371 

372 return ml, mr, mt, mb 

373 

374 def _setup_gmt(self): 

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

376 scaler = self._scaler 

377 

378 if gmtpy.is_gmt5(self._gmtversion): 

379 gmtconf = dict( 

380 MAP_TICK_PEN_PRIMARY='1.25p', 

381 MAP_TICK_PEN_SECONDARY='1.25p', 

382 MAP_TICK_LENGTH_PRIMARY='0.2c', 

383 MAP_TICK_LENGTH_SECONDARY='0.6c', 

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

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

386 PS_CHAR_ENCODING='ISOLatin1+', 

387 MAP_FRAME_TYPE='fancy', 

388 FORMAT_GEO_MAP='D', 

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

390 w*gmtpy.cm, 

391 h*gmtpy.cm), 

392 PS_PAGE_ORIENTATION='portrait', 

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

394 MAP_ANNOT_OBLIQUE='6') 

395 else: 

396 gmtconf = dict( 

397 TICK_PEN='1.25p', 

398 TICK_LENGTH='0.2c', 

399 ANNOT_FONT_PRIMARY='1', 

400 ANNOT_FONT_SIZE_PRIMARY='12p', 

401 LABEL_FONT='1', 

402 LABEL_FONT_SIZE='12p', 

403 CHAR_ENCODING='ISOLatin1+', 

404 BASEMAP_TYPE='fancy', 

405 PLOT_DEGREE_FORMAT='D', 

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

407 w*gmtpy.cm, 

408 h*gmtpy.cm), 

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

410 DOTS_PR_INCH='1200', 

411 OBLIQUE_ANNOTATION='6') 

412 

413 gmtconf.update( 

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

415 

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

417 

418 layout = gmt.default_layout() 

419 

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

421 

422 widget = layout.get_widget() 

423 widget['P'] = widget['J'] 

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

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

426 

427 # aspect = gmtpy.aspect_for_projection( 

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

429 

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

431 widget.set_aspect(aspect) 

432 

433 self._gmt = gmt 

434 self._layout = layout 

435 self._widget = widget 

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

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

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

439 self._have_drawn_axes = False 

440 self._have_drawn_labels = False 

441 

442 def _draw_background(self): 

443 self._have_topo_land = False 

444 self._have_topo_ocean = False 

445 if self.show_topo: 

446 self._have_topo = self._draw_topo() 

447 

448 self._draw_basefeatures() 

449 

450 def _get_topo_tile(self, k): 

451 t = None 

452 demname = None 

453 for dem in self._dems[k]: 

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

455 demname = dem 

456 if t is not None: 

457 break 

458 

459 if not t: 

460 raise NoTopo() 

461 

462 return t, demname 

463 

464 def _prep_topo(self, k): 

465 gmt = self._gmt 

466 t, demname = self._get_topo_tile(k) 

467 

468 if demname not in self._prep_topo_have: 

469 

470 grdfile = gmt.tempfilename() 

471 

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

473 

474 gmtpy.savegrd( 

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

476 

477 if self.illuminate and not is_flat: 

478 if k == 'ocean': 

479 factor = self.illuminate_factor_ocean 

480 else: 

481 factor = self.illuminate_factor_land 

482 

483 ilumfn = gmt.tempfilename() 

484 gmt.grdgradient( 

485 grdfile, 

486 N='e%g' % factor, 

487 A=-45, 

488 G=ilumfn, 

489 out_discard=True) 

490 

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

492 else: 

493 ilumargs = [] 

494 

495 if self.replace_topo_color_only: 

496 t2 = self.replace_topo_color_only 

497 grdfile2 = gmt.tempfilename() 

498 

499 gmtpy.savegrd( 

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

501 naming='lonlat') 

502 

503 if gmt.is_gmt5(): 

504 gmt.grdsample( 

505 grdfile2, 

506 G=grdfile, 

507 n='l', 

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

509 R=grdfile, 

510 out_discard=True) 

511 else: 

512 gmt.grdsample( 

513 grdfile2, 

514 G=grdfile, 

515 Q='l', 

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

517 R=grdfile, 

518 out_discard=True) 

519 

520 gmt.grdmath( 

521 grdfile, '0.0', 'AND', '=', grdfile2, 

522 out_discard=True) 

523 

524 grdfile = grdfile2 

525 

526 self._prep_topo_have[demname] = grdfile, ilumargs 

527 

528 return self._prep_topo_have[demname] 

529 

530 def _draw_topo(self): 

531 widget = self._widget 

532 scaler = self._scaler 

533 gmt = self._gmt 

534 cres = self._coastline_resolution 

535 minarea = self._minarea 

536 

537 JXY = widget.JXY() 

538 R = scaler.R() 

539 

540 try: 

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

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

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

544 *(ilumargs+JXY+R)) 

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

546 self._have_topo_ocean = True 

547 except NoTopo: 

548 self._have_topo_ocean = False 

549 

550 try: 

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

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

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

554 *(ilumargs+JXY+R)) 

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

556 self._have_topo_land = True 

557 except NoTopo: 

558 self._have_topo_land = False 

559 

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

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

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

563 combi = cpt_merge_wet_dry(wet, dry) 

564 for level in combi.levels: 

565 level.vmin /= km 

566 level.vmax /= km 

567 

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

569 write_cpt(combi, topo_cpt) 

570 

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

572 self.gmt.psscale( 

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

574 0.5*gmtpy.cm), 

575 C=topo_cpt, 

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

577 

578 def _draw_basefeatures(self): 

579 gmt = self._gmt 

580 cres = self._coastline_resolution 

581 rivers = self._rivers 

582 minarea = self._minarea 

583 

584 color_wet = self.color_wet 

585 color_dry = self.color_dry 

586 

587 if self.show_rivers and rivers: 

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

589 else: 

590 rivers = [] 

591 

592 fill = {} 

593 if not self._have_topo_land: 

594 fill['G'] = color_dry 

595 

596 if not self._have_topo_ocean: 

597 fill['S'] = color_wet 

598 

599 if self.show_boundaries: 

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

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

602 

603 gmt.pscoast( 

604 D=cres, 

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

606 A=minarea, 

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

608 

609 if self.show_plates: 

610 self.draw_plates() 

611 

612 def _draw_axes(self): 

613 gmt = self._gmt 

614 scaler = self._scaler 

615 widget = self._widget 

616 

617 if self.axes_layout is None: 

618 if self.lat > 0.0: 

619 axes_layout = 'WSen' 

620 else: 

621 axes_layout = 'WseN' 

622 else: 

623 axes_layout = self.axes_layout 

624 

625 scale_km = gmtpy.nice_value(self.radius/5.) / 1000. 

626 

627 if self.show_center_mark: 

628 gmt.psxy( 

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

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

631 *self._jxyr) 

632 

633 if self.show_grid: 

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

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

636 else: 

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

638 

639 if self.show_scale: 

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

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

642 widget.height()/7., 

643 self.lon, 

644 self.lat, 

645 scale_km) 

646 else: 

647 scale = False 

648 

649 gmt.psbasemap( 

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

651 L=scale, 

652 *self._jxyr) 

653 

654 if self.comment: 

655 font_size = self.gmt.label_font_size() 

656 

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

658 if gmt.is_gmt5(): 

659 row = [ 

660 1, 0, 

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

662 self.comment] 

663 

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

665 else: 

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

667 farg = [] 

668 

669 gmt.pstext( 

670 in_rows=[row], 

671 N=True, 

672 R=(0, 1, 0, 1), 

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

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

675 

676 def draw_axes(self): 

677 if not self._have_drawn_axes: 

678 self._draw_axes() 

679 self._have_drawn_axes = True 

680 

681 def _have_coastlines(self): 

682 gmt = self._gmt 

683 cres = self._coastline_resolution 

684 minarea = self._minarea 

685 

686 checkfile = gmt.tempfilename() 

687 

688 gmt.pscoast( 

689 M=True, 

690 D=cres, 

691 W='thinnest,black', 

692 A=minarea, 

693 out_filename=checkfile, 

694 *self._jxyr) 

695 

696 points = [] 

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

698 for line in f: 

699 ls = line.strip() 

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

701 continue 

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

703 points.append((plat, plon)) 

704 

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

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

707 

708 def have_coastlines(self): 

709 self.gmt 

710 return self._have_coastlines() 

711 

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

713 onepoint = False 

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

715 lats = [lats] 

716 lons = [lons] 

717 onepoint = True 

718 

719 if jr is not None: 

720 j, r = jr 

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

722 else: 

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

724 gmt = self.gmt 

725 

726 f = BytesIO() 

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

728 f.seek(0) 

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

730 xs, ys = data.T 

731 if onepoint: 

732 xs = xs[0] 

733 ys = ys[0] 

734 return xs, ys 

735 

736 def _map_box(self, jr=None): 

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

738 

739 xs_corner, ys_corner = self.project( 

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

741 

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

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

744 

745 return w, h 

746 

747 def _map_aspect(self, jr=None): 

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

749 return h/w 

750 

751 def _draw_labels(self): 

752 points_taken = [] 

753 regions_taken = [] 

754 

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

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

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

758 return xx 

759 

760 def roverlaps(a, b): 

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

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

763 

764 w, h = self._map_box() 

765 

766 label_font_size = self.gmt.label_font_size() 

767 

768 if self._labels: 

769 

770 n = len(self._labels) 

771 

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

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

774 

775 font_sizes = [ 

776 (font_size or label_font_size) for font_size in font_sizes] 

777 

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

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

780 

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

782 

783 points_taken.append((xs, ys)) 

784 

785 dxs = num.zeros(n) 

786 dys = num.zeros(n) 

787 

788 for i in range(n): 

789 dx, dy = gmtpy.text_box( 

790 texts[i], 

791 font=fonts[i], 

792 font_size=font_sizes[i], 

793 **styles[i]) 

794 

795 dxs[i] = dx 

796 dys[i] = dy 

797 

798 la = num.logical_and 

799 anchors_ok = ( 

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

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

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

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

804 ) 

805 

806 arects = [ 

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

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

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

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

811 

812 for i in range(n): 

813 for ianch in range(4): 

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

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

816 

817 anchor_choices = [] 

818 anchor_take = [] 

819 for i in range(n): 

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

821 if anchors_ok[ianch][i]] 

822 anchor_choices.append(choices) 

823 if choices: 

824 anchor_take.append(choices[0]) 

825 else: 

826 anchor_take.append(None) 

827 

828 def cost(anchor_take): 

829 noverlaps = 0 

830 for i in range(n): 

831 for j in range(n): 

832 if i != j: 

833 i_take = anchor_take[i] 

834 j_take = anchor_take[j] 

835 if i_take is None or j_take is None: 

836 continue 

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

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

839 if roverlaps(r_i, r_j): 

840 noverlaps += 1 

841 

842 return noverlaps 

843 

844 cur_cost = cost(anchor_take) 

845 imax = 30 

846 while cur_cost != 0 and imax > 0: 

847 for i in range(n): 

848 for t in anchor_choices[i]: 

849 anchor_take_new = list(anchor_take) 

850 anchor_take_new[i] = t 

851 new_cost = cost(anchor_take_new) 

852 if new_cost < cur_cost: 

853 anchor_take = anchor_take_new 

854 cur_cost = new_cost 

855 

856 imax -= 1 

857 

858 while cur_cost != 0: 

859 for i in range(n): 

860 anchor_take_new = list(anchor_take) 

861 anchor_take_new[i] = None 

862 new_cost = cost(anchor_take_new) 

863 if new_cost < cur_cost: 

864 anchor_take = anchor_take_new 

865 cur_cost = new_cost 

866 break 

867 

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

869 

870 for i in range(n): 

871 ianchor = anchor_take[i] 

872 color = colors[i] 

873 if color is None: 

874 color = 'black' 

875 

876 if ianchor is not None: 

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

878 

879 anchor = anchor_strs[ianchor] 

880 

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

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

883 if self.gmt.is_gmt5(): 

884 row = ( 

885 lons[i], lats[i], 

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

887 anchor, 

888 texts[i]) 

889 

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

891 else: 

892 row = ( 

893 lons[i], lats[i], 

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

895 texts[i]) 

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

897 

898 self.gmt.pstext( 

899 in_rows=[row], 

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

901 *(self.jxyr + farg), 

902 **styles[i]) 

903 

904 if self._area_labels: 

905 

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

907 self._area_labels: 

908 

909 if font_size is None: 

910 font_size = label_font_size 

911 

912 if color is None: 

913 color = 'black' 

914 

915 if self.gmt.is_gmt5(): 

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

917 else: 

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

919 

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

921 dx, dy = gmtpy.text_box( 

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

923 

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

925 

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

927 

928 for iloc in range(xs.size): 

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

930 

931 locs_ok[iloc] = True 

932 locs_ok[iloc] &= ( 

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

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

935 

936 overlap = False 

937 for r in regions_taken: 

938 if roverlaps(r, rcandi): 

939 overlap = True 

940 break 

941 

942 locs_ok[iloc] &= not overlap 

943 

944 for xs_taken, ys_taken in points_taken: 

945 locs_ok[iloc] &= no_points_in_rect( 

946 xs_taken, ys_taken, *rcandi) 

947 

948 if not locs_ok[iloc]: 

949 break 

950 

951 rows = [] 

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

953 if not locs_ok[iloc]: 

954 continue 

955 

956 if self.gmt.is_gmt5(): 

957 row = ( 

958 lon, lat, 

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

960 'MC', 

961 text) 

962 

963 else: 

964 row = ( 

965 lon, lat, 

966 font_size, 0, font, 'MC', 

967 text) 

968 

969 rows.append(row) 

970 

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

972 break 

973 

974 self.gmt.pstext( 

975 in_rows=rows, 

976 *(self.jxyr + farg), 

977 **style) 

978 

979 def draw_labels(self): 

980 self.gmt 

981 if not self._have_drawn_labels: 

982 self._draw_labels() 

983 self._have_drawn_labels = True 

984 

985 def add_label( 

986 self, lat, lon, text, 

987 offset_x=5., offset_y=5., 

988 color=None, 

989 font='1', 

990 font_size=None, 

991 angle=0, 

992 style={}): 

993 

994 if 'G' in style: 

995 style = style.copy() 

996 color = style.pop('G') 

997 

998 self._labels.append( 

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

1000 angle, style)) 

1001 

1002 def add_area_label( 

1003 self, lat, lon, text, 

1004 color=None, 

1005 font='3', 

1006 font_size=None, 

1007 style={}): 

1008 

1009 self._area_labels.append( 

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

1011 

1012 def cities_in_region(self): 

1013 from pyrocko.dataset import geonames 

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

1015 cities.extend(self.custom_cities) 

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

1017 return cities 

1018 

1019 def draw_cities(self, 

1020 exact=None, 

1021 include=[], 

1022 exclude=[], 

1023 nmax_soft=10, 

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

1025 

1026 cities = self.cities_in_region() 

1027 

1028 if exact is not None: 

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

1030 minpop = None 

1031 else: 

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

1033 minpop = 10**3 

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

1035 cities_new = [ 

1036 c for c in cities 

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

1038 

1039 if len(cities_new) == 0 or ( 

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

1041 break 

1042 

1043 cities = cities_new 

1044 minpop = minpop_new 

1045 if len(cities) <= nmax_soft: 

1046 break 

1047 

1048 if cities: 

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

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

1051 

1052 self.gmt.psxy( 

1053 in_columns=(lons, lats), 

1054 *self.jxyr, **psxy_style) 

1055 

1056 for c in cities: 

1057 try: 

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

1059 except UnicodeEncodeError: 

1060 text = c.asciiname 

1061 

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

1063 

1064 self._cities_minpop = minpop 

1065 

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

1067 

1068 default_psxy_style = { 

1069 'S': 't8p', 

1070 'G': 'black' 

1071 } 

1072 default_psxy_style.update(psxy_style) 

1073 

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

1075 

1076 self.gmt.psxy( 

1077 in_columns=(lons, lats), 

1078 *self.jxyr, **default_psxy_style) 

1079 

1080 for station in stations: 

1081 self.add_label( 

1082 station.effective_lat, 

1083 station.effective_lon, 

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

1085 

1086 def add_kite_scene(self, scene): 

1087 tile = FloatTile( 

1088 scene.frame.llLon, 

1089 scene.frame.llLat, 

1090 scene.frame.dLon, 

1091 scene.frame.dLat, 

1092 scene.displacement) 

1093 

1094 return tile 

1095 

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

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

1098 

1099 stations = campaign.stations 

1100 

1101 if offset_scale is None: 

1102 offset_scale = num.zeros(campaign.nstations) 

1103 for ista, sta in enumerate(stations): 

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

1105 offset_scale[ista] += comp.shift 

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

1107 

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

1109 scale = (size/10.) / offset_scale 

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

1111 offset_scale, scale) 

1112 

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

1114 

1115 if vertical: 

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

1117 0., -s.up.shift, 

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

1119 s.up.sigma, 0., 

1120 s.code if labels else None] 

1121 for ista, s in enumerate(stations) 

1122 if s.up is not None] 

1123 

1124 else: 

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

1126 -s.east.shift, -s.north.shift, 

1127 s.east.sigma, s.north.sigma, s.correlation_ne, 

1128 s.code if labels else None] 

1129 for ista, s in enumerate(stations) 

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

1131 

1132 default_psxy_style = { 

1133 'h': 0, 

1134 'W': '2p,black', 

1135 'A': '+p2p,black+b+a40', 

1136 'G': 'black', 

1137 'L': True, 

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

1139 } 

1140 

1141 if not labels: 

1142 for row in rows: 

1143 row.pop(-1) 

1144 

1145 if psxy_style is not None: 

1146 default_psxy_style.update(psxy_style) 

1147 

1148 self.gmt.psvelo( 

1149 in_rows=rows, 

1150 *self.jxyr, 

1151 **default_psxy_style) 

1152 

1153 def draw_plates(self): 

1154 from pyrocko.dataset import tectonics 

1155 

1156 neast = 20 

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

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

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

1160 norths2 = num.repeat(norths, neast) 

1161 easts2 = num.tile(easts, nnorth) 

1162 lats, lons = od.ne_to_latlon( 

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

1164 

1165 bird = tectonics.PeterBird2003() 

1166 plates = bird.get_plates() 

1167 

1168 color_plates = gmtpy.color('aluminium5') 

1169 color_velocities = gmtpy.color('skyblue1') 

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

1171 

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

1173 used = [] 

1174 for plate in plates: 

1175 mask = plate.contains_points(points) 

1176 if num.any(mask): 

1177 used.append((plate, mask)) 

1178 

1179 if len(used) > 1: 

1180 

1181 candi_fixed = {} 

1182 

1183 label_data = [] 

1184 for plate, mask in used: 

1185 

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

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

1188 iorder = num.argsort(num.sqrt( 

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

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

1191 

1192 lat_candis = lats[mask][iorder] 

1193 lon_candis = lons[mask][iorder] 

1194 

1195 candi_fixed[plate.name] = lat_candis.size 

1196 

1197 label_data.append(( 

1198 lat_candis, lon_candis, plate, color_plates)) 

1199 

1200 boundaries = bird.get_boundaries() 

1201 

1202 size = 1. 

1203 

1204 psxy_kwargs = [] 

1205 

1206 for boundary in boundaries: 

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

1208 for typ, part in boundary.split_types( 

1209 [['SUB'], 

1210 ['OSR', 'CRB'], 

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

1212 

1213 lats, lons = part.T 

1214 

1215 kwargs = {} 

1216 

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

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

1219 

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

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

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

1223 0.45*size, 3.*size) 

1224 elif boundary.kind == '/': 

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

1226 0.45*size, 3.*size) 

1227 

1228 kwargs['G'] = color_plates 

1229 

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

1231 kwargs_bg = {} 

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

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

1234 size * 3, color_plates) 

1235 psxy_kwargs.append(kwargs_bg) 

1236 

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

1238 

1239 psxy_kwargs.append(kwargs) 

1240 

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

1242 if boundary.plate_name2 in candi_fixed: 

1243 candi_fixed[boundary.plate_name2] += \ 

1244 neast*nnorth 

1245 

1246 elif boundary.kind == '/': 

1247 if boundary.plate_name1 in candi_fixed: 

1248 candi_fixed[boundary.plate_name1] += \ 

1249 neast*nnorth 

1250 

1251 candi_fixed = [name for name in sorted( 

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

1253 

1254 candi_fixed.append(None) 

1255 

1256 gsrm = tectonics.GSRM1() 

1257 

1258 for name in candi_fixed: 

1259 if name not in gsrm.plate_names() \ 

1260 and name not in gsrm.plate_alt_names(): 

1261 

1262 continue 

1263 

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

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

1266 

1267 fixed_plate_name = name 

1268 

1269 if self.show_plate_velocities: 

1270 self.gmt.psvelo( 

1271 in_columns=( 

1272 lons, lats, veast, vnorth, veast_err, vnorth_err, 

1273 corr), 

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

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

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

1277 *self.jxyr) 

1278 

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

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

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

1282 self.add_label( 

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

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

1285 style=dict( 

1286 G=color_velocities_lab)) 

1287 

1288 break 

1289 

1290 if self.show_plate_names: 

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

1292 full_name = bird.full_name(plate.name) 

1293 if plate.name == fixed_plate_name: 

1294 full_name = '@_' + full_name + '@_' 

1295 

1296 self.add_area_label( 

1297 lat_candis, lon_candis, 

1298 full_name, 

1299 color=color, 

1300 font='3') 

1301 

1302 for kwargs in psxy_kwargs: 

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

1304 

1305 

1306def rand(mi, ma): 

1307 mi = float(mi) 

1308 ma = float(ma) 

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

1310 

1311 

1312def split_region(region): 

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

1314 if east > 180: 

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

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

1317 else: 

1318 return [region] 

1319 

1320 

1321class CPTLevel(Object): 

1322 vmin = Float.T() 

1323 vmax = Float.T() 

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

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

1326 

1327 

1328class CPT(Object): 

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

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

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

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

1333 

1334 def scale(self, vmin, vmax): 

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

1336 for level in self.levels: 

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

1338 (vmax - vmin) + vmin 

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

1340 (vmax - vmin) + vmin 

1341 

1342 def discretize(self, nlevels): 

1343 colors = [] 

1344 vals = [] 

1345 for level in self.levels: 

1346 vals.append(level.vmin) 

1347 vals.append(level.vmax) 

1348 colors.append(level.color_min) 

1349 colors.append(level.color_max) 

1350 

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

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

1353 

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

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

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

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

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

1359 

1360 levels = [] 

1361 for ilevel in range(nlevels): 

1362 color = ( 

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

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

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

1366 

1367 levels.append(CPTLevel( 

1368 vmin=x[ilevel], 

1369 vmax=x[ilevel+1], 

1370 color_min=color, 

1371 color_max=color)) 

1372 

1373 cpt = CPT( 

1374 color_below=self.color_below, 

1375 color_above=self.color_above, 

1376 color_nan=self.color_nan, 

1377 levels=levels) 

1378 

1379 return cpt 

1380 

1381 

1382class CPTParseError(Exception): 

1383 pass 

1384 

1385 

1386def read_cpt(filename): 

1387 with open(filename) as f: 

1388 color_below = None 

1389 color_above = None 

1390 color_nan = None 

1391 levels = [] 

1392 try: 

1393 for line in f: 

1394 line = line.strip() 

1395 toks = line.split() 

1396 

1397 if line.startswith('#'): 

1398 continue 

1399 

1400 elif line.startswith('B'): 

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

1402 

1403 elif line.startswith('F'): 

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

1405 

1406 elif line.startswith('N'): 

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

1408 

1409 else: 

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

1411 vmin = values[0] 

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

1413 vmax = values[4] 

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

1415 levels.append(CPTLevel( 

1416 vmin=vmin, 

1417 vmax=vmax, 

1418 color_min=color_min, 

1419 color_max=color_max)) 

1420 

1421 except Exception: 

1422 raise CPTParseError() 

1423 

1424 return CPT( 

1425 color_below=color_below, 

1426 color_above=color_above, 

1427 color_nan=color_nan, 

1428 levels=levels) 

1429 

1430 

1431def color_to_int(color): 

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

1433 

1434 

1435def write_cpt(cpt, filename): 

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

1437 for level in cpt.levels: 

1438 f.write( 

1439 '%e %i %i %i %e %i %i %i\n' % 

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

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

1442 

1443 if cpt.color_below: 

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

1445 

1446 if cpt.color_above: 

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

1448 

1449 if cpt.color_nan: 

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

1451 

1452 

1453def cpt_merge_wet_dry(wet, dry): 

1454 levels = [] 

1455 for level in wet.levels: 

1456 if level.vmin < 0.: 

1457 if level.vmax > 0.: 

1458 level.vmax = 0. 

1459 

1460 levels.append(level) 

1461 

1462 for level in dry.levels: 

1463 if level.vmax > 0.: 

1464 if level.vmin < 0.: 

1465 level.vmin = 0. 

1466 

1467 levels.append(level) 

1468 

1469 combi = CPT( 

1470 color_below=wet.color_below, 

1471 color_above=dry.color_above, 

1472 color_nan=dry.color_nan, 

1473 levels=levels) 

1474 

1475 return combi 

1476 

1477 

1478if __name__ == '__main__': 

1479 from pyrocko import util 

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

1481 

1482 import sys 

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

1484 

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

1486 

1487 for i in range(n): 

1488 m = Map( 

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

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

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

1492 width=30., height=30., 

1493 show_grid=True, 

1494 show_topo=True, 

1495 color_dry=(238, 236, 230), 

1496 topo_cpt_wet='light_sea_uniform', 

1497 topo_cpt_dry='light_land_uniform', 

1498 illuminate=True, 

1499 illuminate_factor_ocean=0.15, 

1500 show_rivers=False, 

1501 show_plates=True) 

1502 

1503 m.draw_cities() 

1504 print(m) 

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