1# http://pyrocko.org - GPLv3
2#
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
105 def x(self):
106 return self.xmin + num.arange(self.nx) * self.dx
108 def y(self):
109 return self.ymin + num.arange(self.ny) * self.dy
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()
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')
128 if population is None:
129 population = 0
130 else:
131 population = int(population)
133 Object.__init__(self, name=name, lat=lat, lon=lon,
134 population=population, asciiname=asciiname)
136 name = Unicode.T()
137 lat = Float.T()
138 lon = Float.T()
139 population = Int.T()
140 asciiname = String.T()
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)
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
205 def save(self, outpath, resolution=75., oversample=2., size=None,
206 width=None, height=None, psconvert=False, crop_eps_mode=False):
208 '''
209 Save the image.
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 '''
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()
233 gmt.save(outpath, resolution=resolution, oversample=oversample,
234 crop_eps_mode=crop_eps_mode,
235 size=size, width=width, height=height, psconvert=psconvert)
237 @property
238 def scaler(self):
239 if self._scaler is None:
240 self._setup_geometry()
242 return self._scaler
244 @property
245 def wesn(self):
246 if self._wesn is None:
247 self._setup_geometry()
249 return self._wesn
251 @property
252 def widget(self):
253 if self._widget is None:
254 self._setup()
256 return self._widget
258 @property
259 def layout(self):
260 if self._layout is None:
261 self._setup()
263 return self._layout
265 @property
266 def jxyr(self):
267 if self._jxyr is None:
268 self._setup()
270 return self._jxyr
272 @property
273 def pxyr(self):
274 if self._pxyr is None:
275 self._setup()
277 return self._pxyr
279 @property
280 def gmt(self):
281 if self._gmt is None:
282 self._setup()
284 if self._have_topo_ocean is None:
285 self._draw_background()
287 return self._gmt
289 def _setup(self):
290 if not self._widget:
291 self._setup_geometry()
293 self._setup_lod()
294 self._setup_gmt()
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
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
309 self._wreg = wreg
310 self._hreg = hreg
312 self._corners = corners(self.lon, self.lat, wreg, hreg)
313 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10)
315 x, y, z = ((west, east), (south, north), (-6000., 4500.))
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)
322 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax))
324 par = scaler.get_params()
326 west = par['xmin']
327 east = par['xmax']
328 south = par['ymin']
329 north = par['ymax']
331 self._wesn = west, east, south, north
332 self._scaler = scaler
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
343 self._minarea = (self.skip_feature_factor * self.radius/km)**2
345 self._coastline_resolution = coastline_resolution
346 self._rivers = rivers
348 self._prep_topo_have = {}
349 self._dems = {}
351 cm2inch = gmtpy.cm/gmtpy.inch
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))
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))
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
375 return ml, mr, mt, mb
377 def _setup_gmt(self):
378 w, h = self.width, self.height
379 scaler = self._scaler
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')
416 gmtconf.update(
417 (k.upper(), v) for (k, v) in self.gmt_config.items())
419 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion)
421 layout = gmt.default_layout()
423 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()])
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
430 # aspect = gmtpy.aspect_for_projection(
431 # gmt.installation['version'], *(widget.J() + scaler.R()))
433 aspect = self._map_aspect(jr=widget.J() + scaler.R())
434 widget.set_aspect(aspect)
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
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()
451 self._draw_basefeatures()
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)
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'
469 if k == 'ocean' and self.topo_wet_path is not None:
470 return self._read_grd(self.topo_wet_path), 'custom'
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
480 if not t:
481 raise NoTopo()
483 return t, demname
485 def _prep_topo(self, k):
486 gmt = self._gmt
487 t, demname = self._get_topo_tile(k)
489 if demname not in self._prep_topo_have:
491 grdfile = gmt.tempfilename()
493 is_flat = num.all(t.data[0] == t.data)
495 gmtpy.savegrd(
496 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat')
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
504 ilumfn = gmt.tempfilename()
505 gmt.grdgradient(
506 grdfile,
507 N='e%g' % factor,
508 A=-45,
509 G=ilumfn,
510 out_discard=True)
512 ilumargs = ['-I%s' % ilumfn]
513 else:
514 ilumargs = []
516 if self.replace_topo_color_only:
517 t2 = self.replace_topo_color_only
518 grdfile2 = gmt.tempfilename()
520 gmtpy.savegrd(
521 t2.x(), t2.y(), t2.data, filename=grdfile2,
522 naming='lonlat')
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)
541 gmt.grdmath(
542 grdfile, '0.0', 'AND', '=', grdfile2,
543 out_discard=True)
545 grdfile = grdfile2
547 self._prep_topo_have[demname] = grdfile, ilumargs
549 return self._prep_topo_have[demname]
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
558 JXY = widget.JXY()
559 R = scaler.R()
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
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
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
589 topo_cpt = self.gmt.tempfilename() + '.cpt'
590 write_cpt(combi, topo_cpt)
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)
599 def _draw_basefeatures(self):
600 gmt = self._gmt
601 cres = self._coastline_resolution
602 rivers = self._rivers
603 minarea = self._minarea
605 color_wet = self.color_wet
606 color_dry = self.color_dry
608 if self.show_rivers and rivers:
609 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)]
610 else:
611 rivers = []
613 fill = {}
614 if not self._have_topo_land:
615 fill['G'] = color_dry
617 if not self._have_topo_ocean:
618 fill['S'] = color_wet
620 if self.show_boundaries:
621 fill['N'] = '1/1p,%s,%s' % (
622 gmtpy.color(self.color_boundaries), 'solid')
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)
630 if self.show_plates:
631 self.draw_plates()
633 def _draw_axes(self):
634 gmt = self._gmt
635 scaler = self._scaler
636 widget = self._widget
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
646 scale_km = nice_value(self.radius/5.) / 1000.
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)
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:'
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
670 gmt.psbasemap(
671 B=(btmpl % scaler.get_params())+axes_layout,
672 L=scale,
673 *self._jxyr)
675 if self.comment:
676 font_size = self.gmt.label_font_size()
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]
685 farg = ['-F+f+j']
686 else:
687 row = [1, 0, font_size, 0, 0, 'BR', self.comment]
688 farg = []
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))
697 def draw_axes(self):
698 if not self._have_drawn_axes:
699 self._draw_axes()
700 self._have_drawn_axes = True
702 def _have_coastlines(self):
703 gmt = self._gmt
704 cres = self._coastline_resolution
705 minarea = self._minarea
707 checkfile = gmt.tempfilename()
709 gmt.pscoast(
710 M=True,
711 D=cres,
712 W='thinnest,black',
713 A=minarea,
714 out_filename=checkfile,
715 *self._jxyr)
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))
726 points = num.array(points, dtype=float)
727 return num.any(points_in_region(points, self._wesn))
729 def have_coastlines(self):
730 self.gmt
731 return self._have_coastlines()
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
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
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
757 def _map_box(self, jr=None):
758 ll_lon, ll_lat, ur_lon, ur_lat = self._corners
760 xs_corner, ys_corner = self.project(
761 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr)
763 w = xs_corner[1] - xs_corner[0]
764 h = ys_corner[1] - ys_corner[0]
766 return w, h
768 def _map_aspect(self, jr=None):
769 w, h = self._map_box(jr=jr)
770 return h/w
772 def _draw_labels(self):
773 points_taken = []
774 regions_taken = []
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
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])
785 w, h = self._map_box()
787 label_font_size = self.gmt.label_font_size()
789 if self._labels:
791 n = len(self._labels)
793 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \
794 angles, styles = list(zip(*self._labels))
796 font_sizes = [
797 (font_size or label_font_size) for font_size in font_sizes]
799 sx = num.array(sx, dtype=float)
800 sy = num.array(sy, dtype=float)
802 xs, ys = self.project(lats, lons)
804 points_taken.append((xs, ys))
806 dxs = num.zeros(n)
807 dys = num.zeros(n)
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])
816 dxs[i] = dx
817 dys[i] = dy
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 )
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)]
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]])
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)
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
863 return noverlaps
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
877 imax -= 1
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
889 anchor_strs = ['BL', 'TR', 'TL', 'BR']
891 for i in range(n):
892 ianchor = anchor_take[i]
893 color = colors[i]
894 if color is None:
895 color = 'black'
897 if ianchor is not None:
898 regions_taken.append([xxx[i] for xxx in arects[ianchor]])
900 anchor = anchor_strs[ianchor]
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])
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]
919 self.gmt.pstext(
920 in_rows=[row],
921 D='%gp/%gp' % (xoff, yoff),
922 *(self.jxyr + farg),
923 **styles[i])
925 if self._area_labels:
927 for lons, lats, text, color, font, font_size, style in \
928 self._area_labels:
930 if font_size is None:
931 font_size = label_font_size
933 if color is None:
934 color = 'black'
936 if self.gmt.is_gmt5():
937 farg = ['-F+f+j']
938 else:
939 farg = ['-G%s' % color]
941 xs, ys = self.project(lats, lons)
942 dx, dy = gmtpy.text_box(
943 text, font=font, font_size=font_size, **style)
945 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy]
947 locs_ok = num.ones(xs.size, dtype=bool)
949 for iloc in range(xs.size):
950 rcandi = [xxx[iloc] for xxx in rects]
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)
957 overlap = False
958 for r in regions_taken:
959 if roverlaps(r, rcandi):
960 overlap = True
961 break
963 locs_ok[iloc] &= not overlap
965 for xs_taken, ys_taken in points_taken:
966 locs_ok[iloc] &= no_points_in_rect(
967 xs_taken, ys_taken, *rcandi)
969 if not locs_ok[iloc]:
970 break
972 rows = []
973 for iloc, (lon, lat) in enumerate(zip(lons, lats)):
974 if not locs_ok[iloc]:
975 continue
977 if self.gmt.is_gmt5():
978 row = (
979 lon, lat,
980 '%i,%s,%s' % (font_size, font, color),
981 'MC',
982 text)
984 else:
985 row = (
986 lon, lat,
987 font_size, 0, font, 'MC',
988 text)
990 rows.append(row)
992 regions_taken.append([xxx[iloc] for xxx in rects])
993 break
995 self.gmt.pstext(
996 in_rows=rows,
997 *(self.jxyr + farg),
998 **style)
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
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={}):
1015 if 'G' in style:
1016 style = style.copy()
1017 color = style.pop('G')
1019 self._labels.append(
1020 (lon, lat, text, offset_x, offset_y, color, font, font_size,
1021 angle, style))
1023 def add_area_label(
1024 self, lat, lon, text,
1025 color=None,
1026 font='3',
1027 font_size=None,
1028 style={}):
1030 self._area_labels.append(
1031 (lon, lat, text, color, font, font_size, style))
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
1040 def draw_cities(self,
1041 exact=None,
1042 include=[],
1043 exclude=[],
1044 nmax_soft=10,
1045 psxy_style=dict(S='s5p', G='black')):
1047 cities = self.cities_in_region()
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]
1060 if len(cities_new) == 0 or (
1061 len(cities_new) < 3 and len(cities) < nmax_soft*2):
1062 break
1064 cities = cities_new
1065 minpop = minpop_new
1066 if len(cities) <= nmax_soft:
1067 break
1069 if cities:
1070 lats = [c.lat for c in cities]
1071 lons = [c.lon for c in cities]
1073 self.gmt.psxy(
1074 in_columns=(lons, lats),
1075 *self.jxyr, **psxy_style)
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
1083 self.add_label(c.lat, c.lon, text)
1085 self._cities_minpop = minpop
1087 def add_stations(self, stations, psxy_style=dict()):
1089 default_psxy_style = {
1090 'S': 't8p',
1091 'G': 'black'
1092 }
1093 default_psxy_style.update(psxy_style)
1095 lats, lons = zip(*[s.effective_latlon for s in stations])
1097 self.gmt.psxy(
1098 in_columns=(lons, lats),
1099 *self.jxyr, **default_psxy_style)
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))
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)
1115 return tile
1117 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None,
1118 labels=True, vertical=False, fontsize=10):
1120 stations = campaign.stations
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()
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)
1134 lats, lons = zip(*[s.effective_latlon for s in stations])
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'
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]
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]
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 }
1169 if not labels:
1170 for row in rows:
1171 row.pop(-1)
1173 if psxy_style is not None:
1174 default_psxy_style.update(psxy_style)
1176 self.gmt.psvelo(
1177 in_rows=rows,
1178 *self.jxyr,
1179 **default_psxy_style)
1181 def draw_plates(self):
1182 from pyrocko.dataset import tectonics
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)
1193 bird = tectonics.PeterBird2003()
1194 plates = bird.get_plates()
1196 color_plates = gmtpy.color('aluminium5')
1197 color_velocities = gmtpy.color('skyblue1')
1198 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1')))
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))
1207 if len(used) > 1:
1209 candi_fixed = {}
1211 label_data = []
1212 for plate, mask in used:
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))
1220 lat_candis = lats[mask][iorder]
1221 lon_candis = lons[mask][iorder]
1223 candi_fixed[plate.name] = lat_candis.size
1225 label_data.append((
1226 lat_candis, lon_candis, plate, color_plates))
1228 boundaries = bird.get_boundaries()
1230 size = 1.
1232 psxy_kwargs = []
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']]):
1241 lats, lons = part.T
1243 kwargs = {}
1245 kwargs['in_columns'] = (lons, lats)
1246 kwargs['W'] = '%gp,%s' % (size, color_plates)
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)
1256 kwargs['G'] = color_plates
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)
1265 kwargs['W'] = '%gp,%s' % (size * 2, 'white')
1267 psxy_kwargs.append(kwargs)
1269 if boundary.kind == '\\':
1270 if boundary.plate_name2 in candi_fixed:
1271 candi_fixed[boundary.plate_name2] += \
1272 neast*nnorth
1274 elif boundary.kind == '/':
1275 if boundary.plate_name1 in candi_fixed:
1276 candi_fixed[boundary.plate_name1] += \
1277 neast*nnorth
1279 candi_fixed = [name for name in sorted(
1280 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])]
1282 candi_fixed.append(None)
1284 gsrm = tectonics.GSRM1()
1286 for name in candi_fixed:
1287 if name not in gsrm.plate_names() \
1288 and name not in gsrm.plate_alt_names():
1290 continue
1292 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \
1293 gsrm.get_velocities(name, region=self._wesn)
1295 fixed_plate_name = name
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)
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))
1316 break
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 + '@_'
1324 self.add_area_label(
1325 lat_candis, lon_candis,
1326 full_name,
1327 color=color,
1328 font='3')
1330 for kwargs in psxy_kwargs:
1331 self.gmt.psxy(*self.jxyr, **kwargs)
1334def rand(mi, ma):
1335 mi = float(mi)
1336 ma = float(ma)
1337 return random.random() * (ma-mi) + mi
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]
1349if __name__ == '__main__':
1350 from pyrocko import util
1351 util.setup_logging('pyrocko.automap', 'info')
1353 import sys
1354 if len(sys.argv) == 2:
1356 n = int(sys.argv[1])
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)
1374 m.draw_cities()
1375 print(m)
1376 m.save('map_%02i.pdf' % i)