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

1# http://pyrocko.org - GPLv3

3# The Pyrocko Developers, 21st Century

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

6import math

7import random

8import logging

10try:

11 from StringIO import StringIO as BytesIO

12except ImportError:

13 from io import BytesIO

15import numpy as num

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

29CPT, CPTLevel, CPTParseError

31points_in_region = od.points_in_region

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

35earthradius = 6371000.0

36r2d = 180./math.pi

37d2r = 1./r2d

38km = 1000.

39d2m = d2r*earthradius

40m2d = 1./d2m

41cm = gmtpy.cm

44def darken(c, f=0.7):

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

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

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()

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)

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()

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

76class NoTopo(Exception):

77 pass

80class OutOfBounds(Exception):

81 pass

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')

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()

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)