Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/plot/automap.py: 83%

768 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-04 09:52 +0000

1# http://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

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

5 

6import math 

7import random 

8import logging 

9 

10try: 

11 from StringIO import StringIO as BytesIO 

12except ImportError: 

13 from io import BytesIO 

14 

15import numpy as num 

16 

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

18 Unicode, Dict) 

19from pyrocko.guts_array import Array 

20from pyrocko.dataset import topo 

21from pyrocko import orthodrome as od 

22from .cpt import ( 

23 get_cpt_path, CPT, CPTLevel, CPTParseError, read_cpt, write_cpt, 

24 cpt_merge_wet_dry) 

25from . import gmtpy 

26from . import nice_value 

27 

28 

29CPT, CPTLevel, CPTParseError 

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 = str(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 topo_dry_path = String.T(optional=True) 

179 topo_wet_path = String.T(optional=True) 

180 axes_layout = String.T(optional=True) 

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

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

183 comment = String.T(optional=True) 

184 approx_ticks = Int.T(default=4) 

185 

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

187 Object.__init__(self, **kwargs) 

188 self._gmt = None 

189 self._scaler = None 

190 self._widget = None 

191 self._corners = None 

192 self._wesn = None 

193 self._minarea = None 

194 self._coastline_resolution = None 

195 self._rivers = None 

196 self._dems = None 

197 self._have_topo_land = None 

198 self._have_topo_ocean = None 

199 self._jxyr = None 

200 self._prep_topo_have = None 

201 self._labels = [] 

202 self._area_labels = [] 

203 self._gmtversion = gmtversion 

204 

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

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

207 

208 ''' 

209 Save the image. 

210 

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

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

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

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

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

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

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

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

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

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

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

222 square with given side length. To save transparency use 

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

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

225 ''' 

226 

227 gmt = self.gmt 

228 self.draw_labels() 

229 self.draw_axes() 

230 if self.show_topo and self.show_topo_scale: 

231 self._draw_topo_scale() 

232 

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

234 crop_eps_mode=crop_eps_mode, 

235 size=size, width=width, height=height, psconvert=psconvert) 

236 

237 @property 

238 def scaler(self): 

239 if self._scaler is None: 

240 self._setup_geometry() 

241 

242 return self._scaler 

243 

244 @property 

245 def wesn(self): 

246 if self._wesn is None: 

247 self._setup_geometry() 

248 

249 return self._wesn 

250 

251 @property 

252 def widget(self): 

253 if self._widget is None: 

254 self._setup() 

255 

256 return self._widget 

257 

258 @property 

259 def layout(self): 

260 if self._layout is None: 

261 self._setup() 

262 

263 return self._layout 

264 

265 @property 

266 def jxyr(self): 

267 if self._jxyr is None: 

268 self._setup() 

269 

270 return self._jxyr 

271 

272 @property 

273 def pxyr(self): 

274 if self._pxyr is None: 

275 self._setup() 

276 

277 return self._pxyr 

278 

279 @property 

280 def gmt(self): 

281 if self._gmt is None: 

282 self._setup() 

283 

284 if self._have_topo_ocean is None: 

285 self._draw_background() 

286 

287 return self._gmt 

288 

289 def _setup(self): 

290 if not self._widget: 

291 self._setup_geometry() 

292 

293 self._setup_lod() 

294 self._setup_gmt() 

295 

296 def _setup_geometry(self): 

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

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

299 wpage -= ml + mr 

300 hpage -= mt + mb 

301 

302 wreg = self.radius * 2.0 

303 hreg = self.radius * 2.0 

304 if wpage >= hpage: 

305 wreg *= wpage/hpage 

306 else: 

307 hreg *= hpage/wpage 

308 

309 self._wreg = wreg 

310 self._hreg = hreg 

311 

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

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

314 

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

316 

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

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

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

320 scaled_unit='km', scaled_unit_factor=0.001) 

321 

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

323 

324 par = scaler.get_params() 

325 

326 west = par['xmin'] 

327 east = par['xmax'] 

328 south = par['ymin'] 

329 north = par['ymax'] 

330 

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

332 self._scaler = scaler 

333 

334 def _setup_lod(self): 

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

336 if self.radius > 1500.*km: 

337 coastline_resolution = 'i' 

338 rivers = False 

339 else: 

340 coastline_resolution = 'f' 

341 rivers = True 

342 

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

344 

345 self._coastline_resolution = coastline_resolution 

346 self._rivers = rivers 

347 

348 self._prep_topo_have = {} 

349 self._dems = {} 

350 

351 cm2inch = gmtpy.cm/gmtpy.inch 

352 

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

354 (self.height * cm2inch)) 

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

356 (self.height * cm2inch)) 

357 

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

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

360 if self._dems[k]: 

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

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

363 

364 def _expand_margins(self): 

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

366 ml = mr = mt = mb = 2.0 

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

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

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

370 ml = mr = self.margins[0] 

371 mt = mb = self.margins[1] 

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

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

374 

375 return ml, mr, mt, mb 

376 

377 def _setup_gmt(self): 

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

379 scaler = self._scaler 

380 

381 if gmtpy.is_gmt5(self._gmtversion): 

382 gmtconf = dict( 

383 MAP_TICK_PEN_PRIMARY='1.25p', 

384 MAP_TICK_PEN_SECONDARY='1.25p', 

385 MAP_TICK_LENGTH_PRIMARY='0.2c', 

386 MAP_TICK_LENGTH_SECONDARY='0.6c', 

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

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

389 PS_CHAR_ENCODING='ISOLatin1+', 

390 MAP_FRAME_TYPE='fancy', 

391 FORMAT_GEO_MAP='D', 

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

393 w*gmtpy.cm, 

394 h*gmtpy.cm), 

395 PS_PAGE_ORIENTATION='portrait', 

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

397 MAP_ANNOT_OBLIQUE='6') 

398 else: 

399 gmtconf = dict( 

400 TICK_PEN='1.25p', 

401 TICK_LENGTH='0.2c', 

402 ANNOT_FONT_PRIMARY='1', 

403 ANNOT_FONT_SIZE_PRIMARY='12p', 

404 LABEL_FONT='1', 

405 LABEL_FONT_SIZE='12p', 

406 CHAR_ENCODING='ISOLatin1+', 

407 BASEMAP_TYPE='fancy', 

408 PLOT_DEGREE_FORMAT='D', 

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

410 w*gmtpy.cm, 

411 h*gmtpy.cm), 

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

413 DOTS_PR_INCH='1200', 

414 OBLIQUE_ANNOTATION='6') 

415 

416 gmtconf.update( 

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

418 

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

420 

421 layout = gmt.default_layout() 

422 

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

424 

425 widget = layout.get_widget() 

426 widget['P'] = widget['J'] 

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

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

429 

430 # aspect = gmtpy.aspect_for_projection( 

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

432 

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

434 widget.set_aspect(aspect) 

435 

436 self._gmt = gmt 

437 self._layout = layout 

438 self._widget = widget 

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

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

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

442 self._have_drawn_axes = False 

443 self._have_drawn_labels = False 

444 

445 def _draw_background(self): 

446 self._have_topo_land = False 

447 self._have_topo_ocean = False 

448 if self.show_topo: 

449 self._have_topo = self._draw_topo() 

450 

451 self._draw_basefeatures() 

452 

453 def _read_grd(self, path): 

454 lons, lats, z = gmtpy.loadgrd(path) 

455 dlons = num.diff(lons) 

456 dlats = num.diff(lats) 

457 dlon = dlons[0] 

458 dlat = dlats[0] 

459 eps = 1e-5 

460 assert num.all(num.abs(dlons - dlon) < dlon*eps) 

461 assert num.all(num.abs(dlats - dlat) < dlat*eps) 

462 return topo.tile.Tile( 

463 lons[0], lats[0], dlon, dlat, z) 

464 

465 def _get_topo_tile(self, k): 

466 if k == 'land' and self.topo_dry_path is not None: 

467 return self._read_grd(self.topo_dry_path), 'custom' 

468 

469 if k == 'ocean' and self.topo_wet_path is not None: 

470 return self._read_grd(self.topo_wet_path), 'custom' 

471 

472 t = None 

473 demname = None 

474 for dem in self._dems[k]: 

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

476 demname = dem 

477 if t is not None: 

478 break 

479 

480 if not t: 

481 raise NoTopo() 

482 

483 return t, demname 

484 

485 def _prep_topo(self, k): 

486 gmt = self._gmt 

487 t, demname = self._get_topo_tile(k) 

488 

489 if demname not in self._prep_topo_have: 

490 

491 grdfile = gmt.tempfilename() 

492 

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

494 

495 gmtpy.savegrd( 

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

497 

498 if self.illuminate and not is_flat: 

499 if k == 'ocean': 

500 factor = self.illuminate_factor_ocean 

501 else: 

502 factor = self.illuminate_factor_land 

503 

504 ilumfn = gmt.tempfilename() 

505 gmt.grdgradient( 

506 grdfile, 

507 N='e%g' % factor, 

508 A=-45, 

509 G=ilumfn, 

510 out_discard=True) 

511 

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

513 else: 

514 ilumargs = [] 

515 

516 if self.replace_topo_color_only: 

517 t2 = self.replace_topo_color_only 

518 grdfile2 = gmt.tempfilename() 

519 

520 gmtpy.savegrd( 

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

522 naming='lonlat') 

523 

524 if gmt.is_gmt5(): 

525 gmt.grdsample( 

526 grdfile2, 

527 G=grdfile, 

528 n='l', 

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

530 R=grdfile, 

531 out_discard=True) 

532 else: 

533 gmt.grdsample( 

534 grdfile2, 

535 G=grdfile, 

536 Q='l', 

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

538 R=grdfile, 

539 out_discard=True) 

540 

541 gmt.grdmath( 

542 grdfile, '0.0', 'AND', '=', grdfile2, 

543 out_discard=True) 

544 

545 grdfile = grdfile2 

546 

547 self._prep_topo_have[demname] = grdfile, ilumargs 

548 

549 return self._prep_topo_have[demname] 

550 

551 def _draw_topo(self): 

552 widget = self._widget 

553 scaler = self._scaler 

554 gmt = self._gmt 

555 cres = self._coastline_resolution 

556 minarea = self._minarea 

557 

558 JXY = widget.JXY() 

559 R = scaler.R() 

560 

561 try: 

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

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

564 gmt.grdimage(grdfile, C=get_cpt_path(self.topo_cpt_wet), 

565 *(ilumargs+JXY+R)) 

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

567 self._have_topo_ocean = True 

568 except NoTopo: 

569 self._have_topo_ocean = False 

570 

571 try: 

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

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

574 gmt.grdimage(grdfile, C=get_cpt_path(self.topo_cpt_dry), 

575 *(ilumargs+JXY+R)) 

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

577 self._have_topo_land = True 

578 except NoTopo: 

579 self._have_topo_land = False 

580 

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

582 dry = read_cpt(get_cpt_path(self.topo_cpt_dry)) 

583 wet = read_cpt(get_cpt_path(self.topo_cpt_wet)) 

584 combi = cpt_merge_wet_dry(wet, dry) 

585 for level in combi.levels: 

586 level.vmin /= km 

587 level.vmax /= km 

588 

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

590 write_cpt(combi, topo_cpt) 

591 

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

593 self.gmt.psscale( 

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

595 0.5*gmtpy.cm), 

596 C=topo_cpt, 

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

598 

599 def _draw_basefeatures(self): 

600 gmt = self._gmt 

601 cres = self._coastline_resolution 

602 rivers = self._rivers 

603 minarea = self._minarea 

604 

605 color_wet = self.color_wet 

606 color_dry = self.color_dry 

607 

608 if self.show_rivers and rivers: 

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

610 else: 

611 rivers = [] 

612 

613 fill = {} 

614 if not self._have_topo_land: 

615 fill['G'] = color_dry 

616 

617 if not self._have_topo_ocean: 

618 fill['S'] = color_wet 

619 

620 if self.show_boundaries: 

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

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

623 

624 gmt.pscoast( 

625 D=cres, 

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

627 A=minarea, 

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

629 

630 if self.show_plates: 

631 self.draw_plates() 

632 

633 def _draw_axes(self): 

634 gmt = self._gmt 

635 scaler = self._scaler 

636 widget = self._widget 

637 

638 if self.axes_layout is None: 

639 if self.lat > 0.0: 

640 axes_layout = 'WSen' 

641 else: 

642 axes_layout = 'WseN' 

643 else: 

644 axes_layout = self.axes_layout 

645 

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

647 

648 if self.show_center_mark: 

649 gmt.psxy( 

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

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

652 *self._jxyr) 

653 

654 if self.show_grid: 

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

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

657 else: 

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

659 

660 if self.show_scale: 

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

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

663 widget.height()/7., 

664 self.lon, 

665 self.lat, 

666 scale_km) 

667 else: 

668 scale = False 

669 

670 gmt.psbasemap( 

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

672 L=scale, 

673 *self._jxyr) 

674 

675 if self.comment: 

676 font_size = self.gmt.label_font_size() 

677 

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

679 if gmt.is_gmt5(): 

680 row = [ 

681 1, 0, 

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

683 self.comment] 

684 

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

686 else: 

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

688 farg = [] 

689 

690 gmt.pstext( 

691 in_rows=[row], 

692 N=True, 

693 R=(0, 1, 0, 1), 

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

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

696 

697 def draw_axes(self): 

698 if not self._have_drawn_axes: 

699 self._draw_axes() 

700 self._have_drawn_axes = True 

701 

702 def _have_coastlines(self): 

703 gmt = self._gmt 

704 cres = self._coastline_resolution 

705 minarea = self._minarea 

706 

707 checkfile = gmt.tempfilename() 

708 

709 gmt.pscoast( 

710 M=True, 

711 D=cres, 

712 W='thinnest,black', 

713 A=minarea, 

714 out_filename=checkfile, 

715 *self._jxyr) 

716 

717 points = [] 

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

719 for line in f: 

720 ls = line.strip() 

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

722 continue 

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

724 points.append((plat, plon)) 

725 

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

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

728 

729 def have_coastlines(self): 

730 self.gmt 

731 return self._have_coastlines() 

732 

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

734 onepoint = False 

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

736 lats = [lats] 

737 lons = [lons] 

738 onepoint = True 

739 

740 if jr is not None: 

741 j, r = jr 

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

743 else: 

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

745 gmt = self.gmt 

746 

747 f = BytesIO() 

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

749 f.seek(0) 

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

751 xs, ys = data.T 

752 if onepoint: 

753 xs = xs[0] 

754 ys = ys[0] 

755 return xs, ys 

756 

757 def _map_box(self, jr=None): 

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

759 

760 xs_corner, ys_corner = self.project( 

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

762 

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

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

765 

766 return w, h 

767 

768 def _map_aspect(self, jr=None): 

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

770 return h/w 

771 

772 def _draw_labels(self): 

773 points_taken = [] 

774 regions_taken = [] 

775 

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

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

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

779 return xx 

780 

781 def roverlaps(a, b): 

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

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

784 

785 w, h = self._map_box() 

786 

787 label_font_size = self.gmt.label_font_size() 

788 

789 if self._labels: 

790 

791 n = len(self._labels) 

792 

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

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

795 

796 font_sizes = [ 

797 (font_size or label_font_size) for font_size in font_sizes] 

798 

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

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

801 

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

803 

804 points_taken.append((xs, ys)) 

805 

806 dxs = num.zeros(n) 

807 dys = num.zeros(n) 

808 

809 for i in range(n): 

810 dx, dy = gmtpy.text_box( 

811 texts[i], 

812 font=fonts[i], 

813 font_size=font_sizes[i], 

814 **styles[i]) 

815 

816 dxs[i] = dx 

817 dys[i] = dy 

818 

819 la = num.logical_and 

820 anchors_ok = ( 

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

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

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

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

825 ) 

826 

827 arects = [ 

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

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

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

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

832 

833 for i in range(n): 

834 for ianch in range(4): 

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

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

837 

838 anchor_choices = [] 

839 anchor_take = [] 

840 for i in range(n): 

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

842 if anchors_ok[ianch][i]] 

843 anchor_choices.append(choices) 

844 if choices: 

845 anchor_take.append(choices[0]) 

846 else: 

847 anchor_take.append(None) 

848 

849 def cost(anchor_take): 

850 noverlaps = 0 

851 for i in range(n): 

852 for j in range(n): 

853 if i != j: 

854 i_take = anchor_take[i] 

855 j_take = anchor_take[j] 

856 if i_take is None or j_take is None: 

857 continue 

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

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

860 if roverlaps(r_i, r_j): 

861 noverlaps += 1 

862 

863 return noverlaps 

864 

865 cur_cost = cost(anchor_take) 

866 imax = 30 

867 while cur_cost != 0 and imax > 0: 

868 for i in range(n): 

869 for t in anchor_choices[i]: 

870 anchor_take_new = list(anchor_take) 

871 anchor_take_new[i] = t 

872 new_cost = cost(anchor_take_new) 

873 if new_cost < cur_cost: 

874 anchor_take = anchor_take_new 

875 cur_cost = new_cost 

876 

877 imax -= 1 

878 

879 while cur_cost != 0: 

880 for i in range(n): 

881 anchor_take_new = list(anchor_take) 

882 anchor_take_new[i] = None 

883 new_cost = cost(anchor_take_new) 

884 if new_cost < cur_cost: 

885 anchor_take = anchor_take_new 

886 cur_cost = new_cost 

887 break 

888 

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

890 

891 for i in range(n): 

892 ianchor = anchor_take[i] 

893 color = colors[i] 

894 if color is None: 

895 color = 'black' 

896 

897 if ianchor is not None: 

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

899 

900 anchor = anchor_strs[ianchor] 

901 

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

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

904 if self.gmt.is_gmt5(): 

905 row = ( 

906 lons[i], lats[i], 

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

908 anchor, 

909 texts[i]) 

910 

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

912 else: 

913 row = ( 

914 lons[i], lats[i], 

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

916 texts[i]) 

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

918 

919 self.gmt.pstext( 

920 in_rows=[row], 

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

922 *(self.jxyr + farg), 

923 **styles[i]) 

924 

925 if self._area_labels: 

926 

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

928 self._area_labels: 

929 

930 if font_size is None: 

931 font_size = label_font_size 

932 

933 if color is None: 

934 color = 'black' 

935 

936 if self.gmt.is_gmt5(): 

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

938 else: 

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

940 

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

942 dx, dy = gmtpy.text_box( 

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

944 

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

946 

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

948 

949 for iloc in range(xs.size): 

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

951 

952 locs_ok[iloc] = True 

953 locs_ok[iloc] &= ( 

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

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

956 

957 overlap = False 

958 for r in regions_taken: 

959 if roverlaps(r, rcandi): 

960 overlap = True 

961 break 

962 

963 locs_ok[iloc] &= not overlap 

964 

965 for xs_taken, ys_taken in points_taken: 

966 locs_ok[iloc] &= no_points_in_rect( 

967 xs_taken, ys_taken, *rcandi) 

968 

969 if not locs_ok[iloc]: 

970 break 

971 

972 rows = [] 

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

974 if not locs_ok[iloc]: 

975 continue 

976 

977 if self.gmt.is_gmt5(): 

978 row = ( 

979 lon, lat, 

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

981 'MC', 

982 text) 

983 

984 else: 

985 row = ( 

986 lon, lat, 

987 font_size, 0, font, 'MC', 

988 text) 

989 

990 rows.append(row) 

991 

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

993 break 

994 

995 self.gmt.pstext( 

996 in_rows=rows, 

997 *(self.jxyr + farg), 

998 **style) 

999 

1000 def draw_labels(self): 

1001 self.gmt 

1002 if not self._have_drawn_labels: 

1003 self._draw_labels() 

1004 self._have_drawn_labels = True 

1005 

1006 def add_label( 

1007 self, lat, lon, text, 

1008 offset_x=5., offset_y=5., 

1009 color=None, 

1010 font='1', 

1011 font_size=None, 

1012 angle=0, 

1013 style={}): 

1014 

1015 if 'G' in style: 

1016 style = style.copy() 

1017 color = style.pop('G') 

1018 

1019 self._labels.append( 

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

1021 angle, style)) 

1022 

1023 def add_area_label( 

1024 self, lat, lon, text, 

1025 color=None, 

1026 font='3', 

1027 font_size=None, 

1028 style={}): 

1029 

1030 self._area_labels.append( 

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

1032 

1033 def cities_in_region(self): 

1034 from pyrocko.dataset import geonames 

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

1036 cities.extend(self.custom_cities) 

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

1038 return cities 

1039 

1040 def draw_cities(self, 

1041 exact=None, 

1042 include=[], 

1043 exclude=[], 

1044 nmax_soft=10, 

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

1046 

1047 cities = self.cities_in_region() 

1048 

1049 if exact is not None: 

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

1051 minpop = None 

1052 else: 

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

1054 minpop = 10**3 

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

1056 cities_new = [ 

1057 c for c in cities 

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

1059 

1060 if len(cities_new) == 0 or ( 

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

1062 break 

1063 

1064 cities = cities_new 

1065 minpop = minpop_new 

1066 if len(cities) <= nmax_soft: 

1067 break 

1068 

1069 if cities: 

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

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

1072 

1073 self.gmt.psxy( 

1074 in_columns=(lons, lats), 

1075 *self.jxyr, **psxy_style) 

1076 

1077 for c in cities: 

1078 try: 

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

1080 except UnicodeEncodeError: 

1081 text = c.asciiname 

1082 

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

1084 

1085 self._cities_minpop = minpop 

1086 

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

1088 

1089 default_psxy_style = { 

1090 'S': 't8p', 

1091 'G': 'black' 

1092 } 

1093 default_psxy_style.update(psxy_style) 

1094 

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

1096 

1097 self.gmt.psxy( 

1098 in_columns=(lons, lats), 

1099 *self.jxyr, **default_psxy_style) 

1100 

1101 for station in stations: 

1102 self.add_label( 

1103 station.effective_lat, 

1104 station.effective_lon, 

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

1106 

1107 def add_kite_scene(self, scene): 

1108 tile = FloatTile( 

1109 scene.frame.llLon, 

1110 scene.frame.llLat, 

1111 scene.frame.dLon, 

1112 scene.frame.dLat, 

1113 scene.displacement) 

1114 

1115 return tile 

1116 

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

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

1119 

1120 stations = campaign.stations 

1121 

1122 if offset_scale is None: 

1123 offset_scale = num.zeros(campaign.nstations) 

1124 for ista, sta in enumerate(stations): 

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

1126 offset_scale[ista] += comp.shift 

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

1128 

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

1130 scale = (size/10.) / offset_scale 

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

1132 offset_scale, scale) 

1133 

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

1135 

1136 if self.gmt.is_gmt6(): 

1137 sign_factor = 1. 

1138 arrow_head_placement = 'e' 

1139 else: 

1140 sign_factor = -1. 

1141 arrow_head_placement = 'b' 

1142 

1143 if vertical: 

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

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

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

1147 s.up.sigma, 0., 

1148 s.code if labels else None] 

1149 for ista, s in enumerate(stations) 

1150 if s.up is not None] 

1151 

1152 else: 

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

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

1155 s.east.sigma, s.north.sigma, s.correlation_ne, 

1156 s.code if labels else None] 

1157 for ista, s in enumerate(stations) 

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

1159 

1160 default_psxy_style = { 

1161 'h': 0, 

1162 'W': '2p,black', 

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

1164 'G': 'black', 

1165 'L': True, 

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

1167 } 

1168 

1169 if not labels: 

1170 for row in rows: 

1171 row.pop(-1) 

1172 

1173 if psxy_style is not None: 

1174 default_psxy_style.update(psxy_style) 

1175 

1176 self.gmt.psvelo( 

1177 in_rows=rows, 

1178 *self.jxyr, 

1179 **default_psxy_style) 

1180 

1181 def draw_plates(self): 

1182 from pyrocko.dataset import tectonics 

1183 

1184 neast = 20 

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

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

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

1188 norths2 = num.repeat(norths, neast) 

1189 easts2 = num.tile(easts, nnorth) 

1190 lats, lons = od.ne_to_latlon( 

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

1192 

1193 bird = tectonics.PeterBird2003() 

1194 plates = bird.get_plates() 

1195 

1196 color_plates = gmtpy.color('aluminium5') 

1197 color_velocities = gmtpy.color('skyblue1') 

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

1199 

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

1201 used = [] 

1202 for plate in plates: 

1203 mask = plate.contains_points(points) 

1204 if num.any(mask): 

1205 used.append((plate, mask)) 

1206 

1207 if len(used) > 1: 

1208 

1209 candi_fixed = {} 

1210 

1211 label_data = [] 

1212 for plate, mask in used: 

1213 

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

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

1216 iorder = num.argsort(num.sqrt( 

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

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

1219 

1220 lat_candis = lats[mask][iorder] 

1221 lon_candis = lons[mask][iorder] 

1222 

1223 candi_fixed[plate.name] = lat_candis.size 

1224 

1225 label_data.append(( 

1226 lat_candis, lon_candis, plate, color_plates)) 

1227 

1228 boundaries = bird.get_boundaries() 

1229 

1230 size = 1. 

1231 

1232 psxy_kwargs = [] 

1233 

1234 for boundary in boundaries: 

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

1236 for typ, part in boundary.split_types( 

1237 [['SUB'], 

1238 ['OSR', 'CRB'], 

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

1240 

1241 lats, lons = part.T 

1242 

1243 kwargs = {} 

1244 

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

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

1247 

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

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

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

1251 0.45*size, 3.*size) 

1252 elif boundary.kind == '/': 

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

1254 0.45*size, 3.*size) 

1255 

1256 kwargs['G'] = color_plates 

1257 

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

1259 kwargs_bg = {} 

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

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

1262 size * 3, color_plates) 

1263 psxy_kwargs.append(kwargs_bg) 

1264 

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

1266 

1267 psxy_kwargs.append(kwargs) 

1268 

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

1270 if boundary.plate_name2 in candi_fixed: 

1271 candi_fixed[boundary.plate_name2] += \ 

1272 neast*nnorth 

1273 

1274 elif boundary.kind == '/': 

1275 if boundary.plate_name1 in candi_fixed: 

1276 candi_fixed[boundary.plate_name1] += \ 

1277 neast*nnorth 

1278 

1279 candi_fixed = [name for name in sorted( 

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

1281 

1282 candi_fixed.append(None) 

1283 

1284 gsrm = tectonics.GSRM1() 

1285 

1286 for name in candi_fixed: 

1287 if name not in gsrm.plate_names() \ 

1288 and name not in gsrm.plate_alt_names(): 

1289 

1290 continue 

1291 

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

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

1294 

1295 fixed_plate_name = name 

1296 

1297 if self.show_plate_velocities: 

1298 self.gmt.psvelo( 

1299 in_columns=( 

1300 lons, lats, veast, vnorth, veast_err, vnorth_err, 

1301 corr), 

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

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

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

1305 *self.jxyr) 

1306 

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

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

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

1310 self.add_label( 

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

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

1313 style=dict( 

1314 G=color_velocities_lab)) 

1315 

1316 break 

1317 

1318 if self.show_plate_names: 

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

1320 full_name = bird.full_name(plate.name) 

1321 if plate.name == fixed_plate_name: 

1322 full_name = '@_' + full_name + '@_' 

1323 

1324 self.add_area_label( 

1325 lat_candis, lon_candis, 

1326 full_name, 

1327 color=color, 

1328 font='3') 

1329 

1330 for kwargs in psxy_kwargs: 

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

1332 

1333 

1334def rand(mi, ma): 

1335 mi = float(mi) 

1336 ma = float(ma) 

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

1338 

1339 

1340def split_region(region): 

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

1342 if east > 180: 

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

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

1345 else: 

1346 return [region] 

1347 

1348 

1349if __name__ == '__main__': 

1350 from pyrocko import util 

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

1352 

1353 import sys 

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

1355 

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

1357 

1358 for i in range(n): 

1359 m = Map( 

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

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

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

1363 width=30., height=30., 

1364 show_grid=True, 

1365 show_topo=True, 

1366 color_dry=(238, 236, 230), 

1367 topo_cpt_wet='light_sea_uniform', 

1368 topo_cpt_dry='light_land_uniform', 

1369 illuminate=True, 

1370 illuminate_factor_ocean=0.15, 

1371 show_rivers=False, 

1372 show_plates=True) 

1373 

1374 m.draw_cities() 

1375 print(m) 

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