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 . import gmtpy
23from . import nice_value
26points_in_region = od.points_in_region
28logger = logging.getLogger('pyrocko.plot.automap')
30earthradius = 6371000.0
31r2d = 180./math.pi
32d2r = 1./r2d
33km = 1000.
34d2m = d2r*earthradius
35m2d = 1./d2m
36cm = gmtpy.cm
39def darken(c, f=0.7):
40 return (c[0]*f, c[1]*f, c[2]*f)
43def corners(lon, lat, w, h):
44 ll_lat, ll_lon = od.ne_to_latlon(lat, lon, -0.5*h, -0.5*w)
45 ur_lat, ur_lon = od.ne_to_latlon(lat, lon, 0.5*h, 0.5*w)
46 return ll_lon, ll_lat, ur_lon, ur_lat
49def extent(lon, lat, w, h, n):
50 x = num.linspace(-0.5*w, 0.5*w, n)
51 y = num.linspace(-0.5*h, 0.5*h, n)
52 slats, slons = od.ne_to_latlon(lat, lon, y[0], x)
53 nlats, nlons = od.ne_to_latlon(lat, lon, y[-1], x)
54 south = slats.min()
55 north = nlats.max()
57 wlats, wlons = od.ne_to_latlon(lat, lon, y, x[0])
58 elats, elons = od.ne_to_latlon(lat, lon, y, x[-1])
59 elons = num.where(elons < wlons, elons + 360., elons)
61 if elons.max() - elons.min() > 180 or wlons.max() - wlons.min() > 180.:
62 west = -180.
63 east = 180.
64 else:
65 west = wlons.min()
66 east = elons.max()
68 return topo.positive_region((west, east, south, north))
71class NoTopo(Exception):
72 pass
75class OutOfBounds(Exception):
76 pass
79class FloatTile(Object):
80 xmin = Float.T()
81 ymin = Float.T()
82 dx = Float.T()
83 dy = Float.T()
84 data = Array.T(shape=(None, None), dtype=float, serialize_as='table')
86 def __init__(self, xmin, ymin, dx, dy, data):
87 Object.__init__(self, init_props=False)
88 self.xmin = float(xmin)
89 self.ymin = float(ymin)
90 self.dx = float(dx)
91 self.dy = float(dy)
92 self.data = data
93 self._set_maxes()
95 def _set_maxes(self):
96 self.ny, self.nx = self.data.shape
97 self.xmax = self.xmin + (self.nx-1) * self.dx
98 self.ymax = self.ymin + (self.ny-1) * self.dy
100 def x(self):
101 return self.xmin + num.arange(self.nx) * self.dx
103 def y(self):
104 return self.ymin + num.arange(self.ny) * self.dy
106 def get(self, x, y):
107 ix = int(round((x - self.xmin) / self.dx))
108 iy = int(round((y - self.ymin) / self.dy))
109 if 0 <= ix < self.nx and 0 <= iy < self.ny:
110 return self.data[iy, ix]
111 else:
112 raise OutOfBounds()
115class City(Object):
116 def __init__(self, name, lat, lon, population=None, asciiname=None):
117 name = str(name)
118 lat = float(lat)
119 lon = float(lon)
120 if asciiname is None:
121 asciiname = name.encode('ascii', errors='replace')
123 if population is None:
124 population = 0
125 else:
126 population = int(population)
128 Object.__init__(self, name=name, lat=lat, lon=lon,
129 population=population, asciiname=asciiname)
131 name = Unicode.T()
132 lat = Float.T()
133 lon = Float.T()
134 population = Int.T()
135 asciiname = String.T()
138class Map(Object):
139 lat = Float.T(optional=True)
140 lon = Float.T(optional=True)
141 radius = Float.T(optional=True)
142 width = Float.T(default=20.)
143 height = Float.T(default=14.)
144 margins = List.T(Float.T())
145 illuminate = Bool.T(default=True)
146 skip_feature_factor = Float.T(default=0.02)
147 show_grid = Bool.T(default=False)
148 show_topo = Bool.T(default=True)
149 show_scale = Bool.T(default=False)
150 show_topo_scale = Bool.T(default=False)
151 show_center_mark = Bool.T(default=False)
152 show_rivers = Bool.T(default=True)
153 show_plates = Bool.T(default=False)
154 show_plate_velocities = Bool.T(default=False)
155 show_plate_names = Bool.T(default=False)
156 show_boundaries = Bool.T(default=False)
157 illuminate_factor_land = Float.T(default=0.5)
158 illuminate_factor_ocean = Float.T(default=0.25)
159 color_wet = Tuple.T(3, Int.T(), default=(216, 242, 254))
160 color_dry = Tuple.T(3, Int.T(), default=(172, 208, 165))
161 color_boundaries = Tuple.T(3, Int.T(), default=(1, 1, 1))
162 topo_resolution_min = Float.T(
163 default=40.,
164 help='minimum resolution of topography [dpi]')
165 topo_resolution_max = Float.T(
166 default=200.,
167 help='maximum resolution of topography [dpi]')
168 replace_topo_color_only = FloatTile.T(
169 optional=True,
170 help='replace topo color while keeping topographic shading')
171 topo_cpt_wet = String.T(default='light_sea')
172 topo_cpt_dry = String.T(default='light_land')
173 topo_dry_path = String.T(optional=True)
174 topo_wet_path = String.T(optional=True)
175 axes_layout = String.T(optional=True)
176 custom_cities = List.T(City.T())
177 gmt_config = Dict.T(String.T(), String.T())
178 comment = String.T(optional=True)
179 approx_ticks = Int.T(default=4)
181 def __init__(self, gmtversion='newest', **kwargs):
182 Object.__init__(self, **kwargs)
183 self._gmt = None
184 self._scaler = None
185 self._widget = None
186 self._corners = None
187 self._wesn = None
188 self._minarea = None
189 self._coastline_resolution = None
190 self._rivers = None
191 self._dems = None
192 self._have_topo_land = None
193 self._have_topo_ocean = None
194 self._jxyr = None
195 self._prep_topo_have = None
196 self._labels = []
197 self._area_labels = []
198 self._gmtversion = gmtversion
200 def save(self, outpath, resolution=75., oversample=2., size=None,
201 width=None, height=None, psconvert=False, crop_eps_mode=False):
203 '''
204 Save the image.
206 Save the image to ``outpath``. The format is determined by the filename
207 extension. Formats are handled as follows: ``'.eps'`` and ``'.ps'``
208 produce EPS and PS, respectively, directly with GMT. If the file name
209 ends with ``'.pdf'``, GMT output is fed through ``gmtpy-epstopdf`` to
210 create a PDF file. For any other filename extension, output is first
211 converted to PDF with ``gmtpy-epstopdf``, then with ``pdftocairo`` to
212 PNG with a resolution oversampled by the factor ``oversample`` and
213 finally the PNG is downsampled and converted to the target format with
214 ``convert``. The resolution of rasterized target image can be
215 controlled either by ``resolution`` in DPI or by specifying ``width``
216 or ``height`` or ``size``, where the latter fits the image into a
217 square with given side length. To save transparency use
218 ``psconvert=True``. To crop the output image with a rectangle to the
219 nearest non-white element set ``crop_eps_mode=True``.
220 '''
222 gmt = self.gmt
223 self.draw_labels()
224 self.draw_axes()
225 if self.show_topo and self.show_topo_scale:
226 self._draw_topo_scale()
228 gmt.save(outpath, resolution=resolution, oversample=oversample,
229 crop_eps_mode=crop_eps_mode,
230 size=size, width=width, height=height, psconvert=psconvert)
232 @property
233 def scaler(self):
234 if self._scaler is None:
235 self._setup_geometry()
237 return self._scaler
239 @property
240 def wesn(self):
241 if self._wesn is None:
242 self._setup_geometry()
244 return self._wesn
246 @property
247 def widget(self):
248 if self._widget is None:
249 self._setup()
251 return self._widget
253 @property
254 def layout(self):
255 if self._layout is None:
256 self._setup()
258 return self._layout
260 @property
261 def jxyr(self):
262 if self._jxyr is None:
263 self._setup()
265 return self._jxyr
267 @property
268 def pxyr(self):
269 if self._pxyr is None:
270 self._setup()
272 return self._pxyr
274 @property
275 def gmt(self):
276 if self._gmt is None:
277 self._setup()
279 if self._have_topo_ocean is None:
280 self._draw_background()
282 return self._gmt
284 def _setup(self):
285 if not self._widget:
286 self._setup_geometry()
288 self._setup_lod()
289 self._setup_gmt()
291 def _setup_geometry(self):
292 wpage, hpage = self.width, self.height
293 ml, mr, mt, mb = self._expand_margins()
294 wpage -= ml + mr
295 hpage -= mt + mb
297 wreg = self.radius * 2.0
298 hreg = self.radius * 2.0
299 if wpage >= hpage:
300 wreg *= wpage/hpage
301 else:
302 hreg *= hpage/wpage
304 self._wreg = wreg
305 self._hreg = hreg
307 self._corners = corners(self.lon, self.lat, wreg, hreg)
308 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10)
310 x, y, z = ((west, east), (south, north), (-6000., 4500.))
312 xax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
313 yax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
314 zax = gmtpy.Ax(mode='min-max', inc=1000., label='Height',
315 scaled_unit='km', scaled_unit_factor=0.001)
317 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax))
319 par = scaler.get_params()
321 west = par['xmin']
322 east = par['xmax']
323 south = par['ymin']
324 north = par['ymax']
326 self._wesn = west, east, south, north
327 self._scaler = scaler
329 def _setup_lod(self):
330 w, e, s, n = self._wesn
331 if self.radius > 1500.*km:
332 coastline_resolution = 'i'
333 rivers = False
334 else:
335 coastline_resolution = 'f'
336 rivers = True
338 self._minarea = (self.skip_feature_factor * self.radius/km)**2
340 self._coastline_resolution = coastline_resolution
341 self._rivers = rivers
343 self._prep_topo_have = {}
344 self._dems = {}
346 cm2inch = gmtpy.cm/gmtpy.inch
348 dmin = 2.0 * self.radius * m2d / (self.topo_resolution_max *
349 (self.height * cm2inch))
350 dmax = 2.0 * self.radius * m2d / (self.topo_resolution_min *
351 (self.height * cm2inch))
353 for k in ['ocean', 'land']:
354 self._dems[k] = topo.select_dem_names(k, dmin, dmax, self._wesn)
355 if self._dems[k]:
356 logger.debug('using topography dataset %s for %s'
357 % (','.join(self._dems[k]), k))
359 def _expand_margins(self):
360 if len(self.margins) == 0 or len(self.margins) > 4:
361 ml = mr = mt = mb = 2.0
362 elif len(self.margins) == 1:
363 ml = mr = mt = mb = self.margins[0]
364 elif len(self.margins) == 2:
365 ml = mr = self.margins[0]
366 mt = mb = self.margins[1]
367 elif len(self.margins) == 4:
368 ml, mr, mt, mb = self.margins
370 return ml, mr, mt, mb
372 def _setup_gmt(self):
373 w, h = self.width, self.height
374 scaler = self._scaler
376 if gmtpy.is_gmt5(self._gmtversion):
377 gmtconf = dict(
378 MAP_TICK_PEN_PRIMARY='1.25p',
379 MAP_TICK_PEN_SECONDARY='1.25p',
380 MAP_TICK_LENGTH_PRIMARY='0.2c',
381 MAP_TICK_LENGTH_SECONDARY='0.6c',
382 FONT_ANNOT_PRIMARY='12p,1,black',
383 FONT_LABEL='12p,1,black',
384 PS_CHAR_ENCODING='ISOLatin1+',
385 MAP_FRAME_TYPE='fancy',
386 FORMAT_GEO_MAP='D',
387 PS_MEDIA='Custom_%ix%i' % (
388 w*gmtpy.cm,
389 h*gmtpy.cm),
390 PS_PAGE_ORIENTATION='portrait',
391 MAP_GRID_PEN_PRIMARY='thinnest,0/50/0',
392 MAP_ANNOT_OBLIQUE='6')
393 else:
394 gmtconf = dict(
395 TICK_PEN='1.25p',
396 TICK_LENGTH='0.2c',
397 ANNOT_FONT_PRIMARY='1',
398 ANNOT_FONT_SIZE_PRIMARY='12p',
399 LABEL_FONT='1',
400 LABEL_FONT_SIZE='12p',
401 CHAR_ENCODING='ISOLatin1+',
402 BASEMAP_TYPE='fancy',
403 PLOT_DEGREE_FORMAT='D',
404 PAPER_MEDIA='Custom_%ix%i' % (
405 w*gmtpy.cm,
406 h*gmtpy.cm),
407 GRID_PEN_PRIMARY='thinnest/0/50/0',
408 DOTS_PR_INCH='1200',
409 OBLIQUE_ANNOTATION='6')
411 gmtconf.update(
412 (k.upper(), v) for (k, v) in self.gmt_config.items())
414 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion)
416 layout = gmt.default_layout()
418 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()])
420 widget = layout.get_widget()
421 widget['P'] = widget['J']
422 widget['J'] = ('-JA%g/%g' % (self.lon, self.lat)) + '/%(width)gp'
423 scaler['R'] = '-R%g/%g/%g/%gr' % self._corners
425 # aspect = gmtpy.aspect_for_projection(
426 # gmt.installation['version'], *(widget.J() + scaler.R()))
428 aspect = self._map_aspect(jr=widget.J() + scaler.R())
429 widget.set_aspect(aspect)
431 self._gmt = gmt
432 self._layout = layout
433 self._widget = widget
434 self._jxyr = self._widget.JXY() + self._scaler.R()
435 self._pxyr = self._widget.PXY() + [
436 '-R%g/%g/%g/%g' % (0, widget.width(), 0, widget.height())]
437 self._have_drawn_axes = False
438 self._have_drawn_labels = False
440 def _draw_background(self):
441 self._have_topo_land = False
442 self._have_topo_ocean = False
443 if self.show_topo:
444 self._have_topo = self._draw_topo()
446 self._draw_basefeatures()
448 def _read_grd(self, path):
449 lons, lats, z = gmtpy.loadgrd(path)
450 dlons = num.diff(lons)
451 dlats = num.diff(lats)
452 dlon = dlons[0]
453 dlat = dlats[0]
454 eps = 1e-5
455 assert num.all(num.abs(dlons - dlon) < dlon*eps)
456 assert num.all(num.abs(dlats - dlat) < dlat*eps)
457 return topo.tile.Tile(
458 lons[0], lats[0], dlon, dlat, z)
460 def _get_topo_tile(self, k):
461 if k == 'land' and self.topo_dry_path is not None:
462 return self._read_grd(self.topo_dry_path), 'custom'
464 if k == 'ocean' and self.topo_wet_path is not None:
465 return self._read_grd(self.topo_wet_path), 'custom'
467 t = None
468 demname = None
469 for dem in self._dems[k]:
470 t = topo.get(dem, self._wesn)
471 demname = dem
472 if t is not None:
473 break
475 if not t:
476 raise NoTopo()
478 return t, demname
480 def _prep_topo(self, k):
481 gmt = self._gmt
482 t, demname = self._get_topo_tile(k)
484 if demname not in self._prep_topo_have:
486 grdfile = gmt.tempfilename()
488 is_flat = num.all(t.data[0] == t.data)
490 gmtpy.savegrd(
491 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat')
493 if self.illuminate and not is_flat:
494 if k == 'ocean':
495 factor = self.illuminate_factor_ocean
496 else:
497 factor = self.illuminate_factor_land
499 ilumfn = gmt.tempfilename()
500 gmt.grdgradient(
501 grdfile,
502 N='e%g' % factor,
503 A=-45,
504 G=ilumfn,
505 out_discard=True)
507 ilumargs = ['-I%s' % ilumfn]
508 else:
509 ilumargs = []
511 if self.replace_topo_color_only:
512 t2 = self.replace_topo_color_only
513 grdfile2 = gmt.tempfilename()
515 gmtpy.savegrd(
516 t2.x(), t2.y(), t2.data, filename=grdfile2,
517 naming='lonlat')
519 if gmt.is_gmt5():
520 gmt.grdsample(
521 grdfile2,
522 G=grdfile,
523 n='l',
524 I='%g/%g' % (t.dx, t.dy), # noqa
525 R=grdfile,
526 out_discard=True)
527 else:
528 gmt.grdsample(
529 grdfile2,
530 G=grdfile,
531 Q='l',
532 I='%g/%g' % (t.dx, t.dy), # noqa
533 R=grdfile,
534 out_discard=True)
536 gmt.grdmath(
537 grdfile, '0.0', 'AND', '=', grdfile2,
538 out_discard=True)
540 grdfile = grdfile2
542 self._prep_topo_have[demname] = grdfile, ilumargs
544 return self._prep_topo_have[demname]
546 def _draw_topo(self):
547 widget = self._widget
548 scaler = self._scaler
549 gmt = self._gmt
550 cres = self._coastline_resolution
551 minarea = self._minarea
553 JXY = widget.JXY()
554 R = scaler.R()
556 try:
557 grdfile, ilumargs = self._prep_topo('ocean')
558 gmt.pscoast(D=cres, S='c', A=minarea, *(JXY+R))
559 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_wet),
560 *(ilumargs+JXY+R))
561 gmt.pscoast(Q=True, *(JXY+R))
562 self._have_topo_ocean = True
563 except NoTopo:
564 self._have_topo_ocean = False
566 try:
567 grdfile, ilumargs = self._prep_topo('land')
568 gmt.pscoast(D=cres, G='c', A=minarea, *(JXY+R))
569 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_dry),
570 *(ilumargs+JXY+R))
571 gmt.pscoast(Q=True, *(JXY+R))
572 self._have_topo_land = True
573 except NoTopo:
574 self._have_topo_land = False
576 def _draw_topo_scale(self, label='Elevation [km]'):
577 dry = read_cpt(topo.cpt(self.topo_cpt_dry))
578 wet = read_cpt(topo.cpt(self.topo_cpt_wet))
579 combi = cpt_merge_wet_dry(wet, dry)
580 for level in combi.levels:
581 level.vmin /= km
582 level.vmax /= km
584 topo_cpt = self.gmt.tempfilename() + '.cpt'
585 write_cpt(combi, topo_cpt)
587 (w, h), (xo, yo) = self.widget.get_size()
588 self.gmt.psscale(
589 D='%gp/%gp/%gp/%gph' % (xo + 0.5*w, yo - 2.0*gmtpy.cm, w,
590 0.5*gmtpy.cm),
591 C=topo_cpt,
592 B='1:%s:' % label)
594 def _draw_basefeatures(self):
595 gmt = self._gmt
596 cres = self._coastline_resolution
597 rivers = self._rivers
598 minarea = self._minarea
600 color_wet = self.color_wet
601 color_dry = self.color_dry
603 if self.show_rivers and rivers:
604 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)]
605 else:
606 rivers = []
608 fill = {}
609 if not self._have_topo_land:
610 fill['G'] = color_dry
612 if not self._have_topo_ocean:
613 fill['S'] = color_wet
615 if self.show_boundaries:
616 fill['N'] = '1/1p,%s,%s' % (
617 gmtpy.color(self.color_boundaries), 'solid')
619 gmt.pscoast(
620 D=cres,
621 W='thinnest,%s' % gmtpy.color(darken(gmtpy.color_tup(color_dry))),
622 A=minarea,
623 *(rivers+self._jxyr), **fill)
625 if self.show_plates:
626 self.draw_plates()
628 def _draw_axes(self):
629 gmt = self._gmt
630 scaler = self._scaler
631 widget = self._widget
633 if self.axes_layout is None:
634 if self.lat > 0.0:
635 axes_layout = 'WSen'
636 else:
637 axes_layout = 'WseN'
638 else:
639 axes_layout = self.axes_layout
641 scale_km = nice_value(self.radius/5.) / 1000.
643 if self.show_center_mark:
644 gmt.psxy(
645 in_rows=[[self.lon, self.lat]],
646 S='c20p', W='2p,black',
647 *self._jxyr)
649 if self.show_grid:
650 btmpl = ('%(xinc)gg%(xinc)g:%(xlabel)s:/'
651 '%(yinc)gg%(yinc)g:%(ylabel)s:')
652 else:
653 btmpl = '%(xinc)g:%(xlabel)s:/%(yinc)g:%(ylabel)s:'
655 if self.show_scale:
656 scale = 'x%gp/%gp/%g/%g/%gk' % (
657 6./7*widget.width(),
658 widget.height()/7.,
659 self.lon,
660 self.lat,
661 scale_km)
662 else:
663 scale = False
665 gmt.psbasemap(
666 B=(btmpl % scaler.get_params())+axes_layout,
667 L=scale,
668 *self._jxyr)
670 if self.comment:
671 font_size = self.gmt.label_font_size()
673 _, east, south, _ = self._wesn
674 if gmt.is_gmt5():
675 row = [
676 1, 0,
677 '%gp,%s,%s' % (font_size, 0, 'black'), 'BR',
678 self.comment]
680 farg = ['-F+f+j']
681 else:
682 row = [1, 0, font_size, 0, 0, 'BR', self.comment]
683 farg = []
685 gmt.pstext(
686 in_rows=[row],
687 N=True,
688 R=(0, 1, 0, 1),
689 D='%gp/%gp' % (-font_size*0.2, font_size*0.3),
690 *(widget.PXY() + farg))
692 def draw_axes(self):
693 if not self._have_drawn_axes:
694 self._draw_axes()
695 self._have_drawn_axes = True
697 def _have_coastlines(self):
698 gmt = self._gmt
699 cres = self._coastline_resolution
700 minarea = self._minarea
702 checkfile = gmt.tempfilename()
704 gmt.pscoast(
705 M=True,
706 D=cres,
707 W='thinnest,black',
708 A=minarea,
709 out_filename=checkfile,
710 *self._jxyr)
712 points = []
713 with open(checkfile, 'r') as f:
714 for line in f:
715 ls = line.strip()
716 if ls.startswith('#') or ls.startswith('>') or ls == '':
717 continue
718 plon, plat = [float(x) for x in ls.split()]
719 points.append((plat, plon))
721 points = num.array(points, dtype=float)
722 return num.any(points_in_region(points, self._wesn))
724 def have_coastlines(self):
725 self.gmt
726 return self._have_coastlines()
728 def project(self, lats, lons, jr=None):
729 onepoint = False
730 if isinstance(lats, float) and isinstance(lons, float):
731 lats = [lats]
732 lons = [lons]
733 onepoint = True
735 if jr is not None:
736 j, r = jr
737 gmt = gmtpy.GMT(version=self._gmtversion)
738 else:
739 j, _, _, r = self.jxyr
740 gmt = self.gmt
742 f = BytesIO()
743 gmt.mapproject(j, r, in_columns=(lons, lats), out_stream=f, D='p')
744 f.seek(0)
745 data = num.loadtxt(f, ndmin=2)
746 xs, ys = data.T
747 if onepoint:
748 xs = xs[0]
749 ys = ys[0]
750 return xs, ys
752 def _map_box(self, jr=None):
753 ll_lon, ll_lat, ur_lon, ur_lat = self._corners
755 xs_corner, ys_corner = self.project(
756 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr)
758 w = xs_corner[1] - xs_corner[0]
759 h = ys_corner[1] - ys_corner[0]
761 return w, h
763 def _map_aspect(self, jr=None):
764 w, h = self._map_box(jr=jr)
765 return h/w
767 def _draw_labels(self):
768 points_taken = []
769 regions_taken = []
771 def no_points_in_rect(xs, ys, xmin, ymin, xmax, ymax):
772 xx = not num.any(la(la(xmin < xs, xs < xmax),
773 la(ymin < ys, ys < ymax)))
774 return xx
776 def roverlaps(a, b):
777 return (a[0] < b[2] and b[0] < a[2] and
778 a[1] < b[3] and b[1] < a[3])
780 w, h = self._map_box()
782 label_font_size = self.gmt.label_font_size()
784 if self._labels:
786 n = len(self._labels)
788 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \
789 angles, styles = list(zip(*self._labels))
791 font_sizes = [
792 (font_size or label_font_size) for font_size in font_sizes]
794 sx = num.array(sx, dtype=float)
795 sy = num.array(sy, dtype=float)
797 xs, ys = self.project(lats, lons)
799 points_taken.append((xs, ys))
801 dxs = num.zeros(n)
802 dys = num.zeros(n)
804 for i in range(n):
805 dx, dy = gmtpy.text_box(
806 texts[i],
807 font=fonts[i],
808 font_size=font_sizes[i],
809 **styles[i])
811 dxs[i] = dx
812 dys[i] = dy
814 la = num.logical_and
815 anchors_ok = (
816 la(xs + sx + dxs < w, ys + sy + dys < h),
817 la(xs - sx - dxs > 0., ys - sy - dys > 0.),
818 la(xs + sx + dxs < w, ys - sy - dys > 0.),
819 la(xs - sx - dxs > 0., ys + sy + dys < h),
820 )
822 arects = [
823 (xs, ys, xs + sx + dxs, ys + sy + dys),
824 (xs - sx - dxs, ys - sy - dys, xs, ys),
825 (xs, ys - sy - dys, xs + sx + dxs, ys),
826 (xs - sx - dxs, ys, xs, ys + sy + dys)]
828 for i in range(n):
829 for ianch in range(4):
830 anchors_ok[ianch][i] &= no_points_in_rect(
831 xs, ys, *[xxx[i] for xxx in arects[ianch]])
833 anchor_choices = []
834 anchor_take = []
835 for i in range(n):
836 choices = [ianch for ianch in range(4)
837 if anchors_ok[ianch][i]]
838 anchor_choices.append(choices)
839 if choices:
840 anchor_take.append(choices[0])
841 else:
842 anchor_take.append(None)
844 def cost(anchor_take):
845 noverlaps = 0
846 for i in range(n):
847 for j in range(n):
848 if i != j:
849 i_take = anchor_take[i]
850 j_take = anchor_take[j]
851 if i_take is None or j_take is None:
852 continue
853 r_i = [xxx[i] for xxx in arects[i_take]]
854 r_j = [xxx[j] for xxx in arects[j_take]]
855 if roverlaps(r_i, r_j):
856 noverlaps += 1
858 return noverlaps
860 cur_cost = cost(anchor_take)
861 imax = 30
862 while cur_cost != 0 and imax > 0:
863 for i in range(n):
864 for t in anchor_choices[i]:
865 anchor_take_new = list(anchor_take)
866 anchor_take_new[i] = t
867 new_cost = cost(anchor_take_new)
868 if new_cost < cur_cost:
869 anchor_take = anchor_take_new
870 cur_cost = new_cost
872 imax -= 1
874 while cur_cost != 0:
875 for i in range(n):
876 anchor_take_new = list(anchor_take)
877 anchor_take_new[i] = None
878 new_cost = cost(anchor_take_new)
879 if new_cost < cur_cost:
880 anchor_take = anchor_take_new
881 cur_cost = new_cost
882 break
884 anchor_strs = ['BL', 'TR', 'TL', 'BR']
886 for i in range(n):
887 ianchor = anchor_take[i]
888 color = colors[i]
889 if color is None:
890 color = 'black'
892 if ianchor is not None:
893 regions_taken.append([xxx[i] for xxx in arects[ianchor]])
895 anchor = anchor_strs[ianchor]
897 yoff = [-sy[i], sy[i]][anchor[0] == 'B']
898 xoff = [-sx[i], sx[i]][anchor[1] == 'L']
899 if self.gmt.is_gmt5():
900 row = (
901 lons[i], lats[i],
902 '%i,%s,%s' % (font_sizes[i], fonts[i], color),
903 anchor,
904 texts[i])
906 farg = ['-F+f+j+a%g' % angles[i]]
907 else:
908 row = (
909 lons[i], lats[i],
910 font_sizes[i], angles[i], fonts[i], anchor,
911 texts[i])
912 farg = ['-G%s' % color]
914 self.gmt.pstext(
915 in_rows=[row],
916 D='%gp/%gp' % (xoff, yoff),
917 *(self.jxyr + farg),
918 **styles[i])
920 if self._area_labels:
922 for lons, lats, text, color, font, font_size, style in \
923 self._area_labels:
925 if font_size is None:
926 font_size = label_font_size
928 if color is None:
929 color = 'black'
931 if self.gmt.is_gmt5():
932 farg = ['-F+f+j']
933 else:
934 farg = ['-G%s' % color]
936 xs, ys = self.project(lats, lons)
937 dx, dy = gmtpy.text_box(
938 text, font=font, font_size=font_size, **style)
940 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy]
942 locs_ok = num.ones(xs.size, dtype=bool)
944 for iloc in range(xs.size):
945 rcandi = [xxx[iloc] for xxx in rects]
947 locs_ok[iloc] = True
948 locs_ok[iloc] &= (
949 0 < rcandi[0] and rcandi[2] < w
950 and 0 < rcandi[1] and rcandi[3] < h)
952 overlap = False
953 for r in regions_taken:
954 if roverlaps(r, rcandi):
955 overlap = True
956 break
958 locs_ok[iloc] &= not overlap
960 for xs_taken, ys_taken in points_taken:
961 locs_ok[iloc] &= no_points_in_rect(
962 xs_taken, ys_taken, *rcandi)
964 if not locs_ok[iloc]:
965 break
967 rows = []
968 for iloc, (lon, lat) in enumerate(zip(lons, lats)):
969 if not locs_ok[iloc]:
970 continue
972 if self.gmt.is_gmt5():
973 row = (
974 lon, lat,
975 '%i,%s,%s' % (font_size, font, color),
976 'MC',
977 text)
979 else:
980 row = (
981 lon, lat,
982 font_size, 0, font, 'MC',
983 text)
985 rows.append(row)
987 regions_taken.append([xxx[iloc] for xxx in rects])
988 break
990 self.gmt.pstext(
991 in_rows=rows,
992 *(self.jxyr + farg),
993 **style)
995 def draw_labels(self):
996 self.gmt
997 if not self._have_drawn_labels:
998 self._draw_labels()
999 self._have_drawn_labels = True
1001 def add_label(
1002 self, lat, lon, text,
1003 offset_x=5., offset_y=5.,
1004 color=None,
1005 font='1',
1006 font_size=None,
1007 angle=0,
1008 style={}):
1010 if 'G' in style:
1011 style = style.copy()
1012 color = style.pop('G')
1014 self._labels.append(
1015 (lon, lat, text, offset_x, offset_y, color, font, font_size,
1016 angle, style))
1018 def add_area_label(
1019 self, lat, lon, text,
1020 color=None,
1021 font='3',
1022 font_size=None,
1023 style={}):
1025 self._area_labels.append(
1026 (lon, lat, text, color, font, font_size, style))
1028 def cities_in_region(self):
1029 from pyrocko.dataset import geonames
1030 cities = geonames.get_cities_region(region=self.wesn, minpop=0)
1031 cities.extend(self.custom_cities)
1032 cities.sort(key=lambda x: x.population)
1033 return cities
1035 def draw_cities(self,
1036 exact=None,
1037 include=[],
1038 exclude=[],
1039 nmax_soft=10,
1040 psxy_style=dict(S='s5p', G='black')):
1042 cities = self.cities_in_region()
1044 if exact is not None:
1045 cities = [c for c in cities if c.name in exact]
1046 minpop = None
1047 else:
1048 cities = [c for c in cities if c.name not in exclude]
1049 minpop = 10**3
1050 for minpop_new in [1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6, 3e6, 1e7]:
1051 cities_new = [
1052 c for c in cities
1053 if c.population > minpop_new or c.name in include]
1055 if len(cities_new) == 0 or (
1056 len(cities_new) < 3 and len(cities) < nmax_soft*2):
1057 break
1059 cities = cities_new
1060 minpop = minpop_new
1061 if len(cities) <= nmax_soft:
1062 break
1064 if cities:
1065 lats = [c.lat for c in cities]
1066 lons = [c.lon for c in cities]
1068 self.gmt.psxy(
1069 in_columns=(lons, lats),
1070 *self.jxyr, **psxy_style)
1072 for c in cities:
1073 try:
1074 text = c.name.encode('iso-8859-1').decode('iso-8859-1')
1075 except UnicodeEncodeError:
1076 text = c.asciiname
1078 self.add_label(c.lat, c.lon, text)
1080 self._cities_minpop = minpop
1082 def add_stations(self, stations, psxy_style=dict()):
1084 default_psxy_style = {
1085 'S': 't8p',
1086 'G': 'black'
1087 }
1088 default_psxy_style.update(psxy_style)
1090 lats, lons = zip(*[s.effective_latlon for s in stations])
1092 self.gmt.psxy(
1093 in_columns=(lons, lats),
1094 *self.jxyr, **default_psxy_style)
1096 for station in stations:
1097 self.add_label(
1098 station.effective_lat,
1099 station.effective_lon,
1100 '.'.join(x for x in (station.network, station.station) if x))
1102 def add_kite_scene(self, scene):
1103 tile = FloatTile(
1104 scene.frame.llLon,
1105 scene.frame.llLat,
1106 scene.frame.dLon,
1107 scene.frame.dLat,
1108 scene.displacement)
1110 return tile
1112 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None,
1113 labels=True, vertical=False, fontsize=10):
1115 stations = campaign.stations
1117 if offset_scale is None:
1118 offset_scale = num.zeros(campaign.nstations)
1119 for ista, sta in enumerate(stations):
1120 for comp in sta.components.values():
1121 offset_scale[ista] += comp.shift
1122 offset_scale = num.sqrt(offset_scale**2).max()
1124 size = math.sqrt(self.height**2 + self.width**2)
1125 scale = (size/10.) / offset_scale
1126 logger.debug('GNSS: Using offset scale %f, map scale %f',
1127 offset_scale, scale)
1129 lats, lons = zip(*[s.effective_latlon for s in stations])
1131 if self.gmt.is_gmt6():
1132 sign_factor = 1.
1133 arrow_head_placement = 'e'
1134 else:
1135 sign_factor = -1.
1136 arrow_head_placement = 'b'
1138 if vertical:
1139 rows = [[lons[ista], lats[ista],
1140 0., sign_factor * s.up.shift,
1141 (s.east.sigma + s.north.sigma) if s.east.sigma else 0.,
1142 s.up.sigma, 0.,
1143 s.code if labels else None]
1144 for ista, s in enumerate(stations)
1145 if s.up is not None]
1147 else:
1148 rows = [[lons[ista], lats[ista],
1149 sign_factor * s.east.shift, sign_factor * s.north.shift,
1150 s.east.sigma, s.north.sigma, s.correlation_ne,
1151 s.code if labels else None]
1152 for ista, s in enumerate(stations)
1153 if s.east is not None or s.north is not None]
1155 default_psxy_style = {
1156 'h': 0,
1157 'W': '2p,black',
1158 'A': '+p2p,black+{}+a40'.format(arrow_head_placement),
1159 'G': 'black',
1160 'L': True,
1161 'S': 'e%dc/0.95/%d' % (scale, fontsize),
1162 }
1164 if not labels:
1165 for row in rows:
1166 row.pop(-1)
1168 if psxy_style is not None:
1169 default_psxy_style.update(psxy_style)
1171 self.gmt.psvelo(
1172 in_rows=rows,
1173 *self.jxyr,
1174 **default_psxy_style)
1176 def draw_plates(self):
1177 from pyrocko.dataset import tectonics
1179 neast = 20
1180 nnorth = max(1, int(round(num.round(self._hreg/self._wreg * neast))))
1181 norths = num.linspace(-self._hreg*0.5, self._hreg*0.5, nnorth)
1182 easts = num.linspace(-self._wreg*0.5, self._wreg*0.5, neast)
1183 norths2 = num.repeat(norths, neast)
1184 easts2 = num.tile(easts, nnorth)
1185 lats, lons = od.ne_to_latlon(
1186 self.lat, self.lon, norths2, easts2)
1188 bird = tectonics.PeterBird2003()
1189 plates = bird.get_plates()
1191 color_plates = gmtpy.color('aluminium5')
1192 color_velocities = gmtpy.color('skyblue1')
1193 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1')))
1195 points = num.vstack((lats, lons)).T
1196 used = []
1197 for plate in plates:
1198 mask = plate.contains_points(points)
1199 if num.any(mask):
1200 used.append((plate, mask))
1202 if len(used) > 1:
1204 candi_fixed = {}
1206 label_data = []
1207 for plate, mask in used:
1209 mean_north = num.mean(norths2[mask])
1210 mean_east = num.mean(easts2[mask])
1211 iorder = num.argsort(num.sqrt(
1212 (norths2[mask] - mean_north)**2 +
1213 (easts2[mask] - mean_east)**2))
1215 lat_candis = lats[mask][iorder]
1216 lon_candis = lons[mask][iorder]
1218 candi_fixed[plate.name] = lat_candis.size
1220 label_data.append((
1221 lat_candis, lon_candis, plate, color_plates))
1223 boundaries = bird.get_boundaries()
1225 size = 1.
1227 psxy_kwargs = []
1229 for boundary in boundaries:
1230 if num.any(points_in_region(boundary.points, self._wesn)):
1231 for typ, part in boundary.split_types(
1232 [['SUB'],
1233 ['OSR', 'CRB'],
1234 ['OTF', 'CTF', 'OCB', 'CCB']]):
1236 lats, lons = part.T
1238 kwargs = {}
1240 kwargs['in_columns'] = (lons, lats)
1241 kwargs['W'] = '%gp,%s' % (size, color_plates)
1243 if typ[0] == 'SUB':
1244 if boundary.kind == '\\':
1245 kwargs['S'] = 'f%g/%gp+t+r' % (
1246 0.45*size, 3.*size)
1247 elif boundary.kind == '/':
1248 kwargs['S'] = 'f%g/%gp+t+l' % (
1249 0.45*size, 3.*size)
1251 kwargs['G'] = color_plates
1253 elif typ[0] in ['OSR', 'CRB']:
1254 kwargs_bg = {}
1255 kwargs_bg['in_columns'] = (lons, lats)
1256 kwargs_bg['W'] = '%gp,%s' % (
1257 size * 3, color_plates)
1258 psxy_kwargs.append(kwargs_bg)
1260 kwargs['W'] = '%gp,%s' % (size * 2, 'white')
1262 psxy_kwargs.append(kwargs)
1264 if boundary.kind == '\\':
1265 if boundary.plate_name2 in candi_fixed:
1266 candi_fixed[boundary.plate_name2] += \
1267 neast*nnorth
1269 elif boundary.kind == '/':
1270 if boundary.plate_name1 in candi_fixed:
1271 candi_fixed[boundary.plate_name1] += \
1272 neast*nnorth
1274 candi_fixed = [name for name in sorted(
1275 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])]
1277 candi_fixed.append(None)
1279 gsrm = tectonics.GSRM1()
1281 for name in candi_fixed:
1282 if name not in gsrm.plate_names() \
1283 and name not in gsrm.plate_alt_names():
1285 continue
1287 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \
1288 gsrm.get_velocities(name, region=self._wesn)
1290 fixed_plate_name = name
1292 if self.show_plate_velocities:
1293 self.gmt.psvelo(
1294 in_columns=(
1295 lons, lats, veast, vnorth, veast_err, vnorth_err,
1296 corr),
1297 W='0.25p,%s' % color_velocities,
1298 A='9p+e+g%s' % color_velocities,
1299 S='e0.2p/0.95/10',
1300 *self.jxyr)
1302 for _ in range(len(lons) // 50 + 1):
1303 ii = random.randint(0, len(lons)-1)
1304 v = math.sqrt(vnorth[ii]**2 + veast[ii]**2)
1305 self.add_label(
1306 lats[ii], lons[ii], '%.0f' % v,
1307 font_size=0.7*self.gmt.label_font_size(),
1308 style=dict(
1309 G=color_velocities_lab))
1311 break
1313 if self.show_plate_names:
1314 for (lat_candis, lon_candis, plate, color) in label_data:
1315 full_name = bird.full_name(plate.name)
1316 if plate.name == fixed_plate_name:
1317 full_name = '@_' + full_name + '@_'
1319 self.add_area_label(
1320 lat_candis, lon_candis,
1321 full_name,
1322 color=color,
1323 font='3')
1325 for kwargs in psxy_kwargs:
1326 self.gmt.psxy(*self.jxyr, **kwargs)
1329def rand(mi, ma):
1330 mi = float(mi)
1331 ma = float(ma)
1332 return random.random() * (ma-mi) + mi
1335def split_region(region):
1336 west, east, south, north = topo.positive_region(region)
1337 if east > 180:
1338 return [(west, 180., south, north),
1339 (-180., east-360., south, north)]
1340 else:
1341 return [region]
1344class CPTLevel(Object):
1345 vmin = Float.T()
1346 vmax = Float.T()
1347 color_min = Tuple.T(3, Float.T())
1348 color_max = Tuple.T(3, Float.T())
1351class CPT(Object):
1352 color_below = Tuple.T(3, Float.T(), optional=True)
1353 color_above = Tuple.T(3, Float.T(), optional=True)
1354 color_nan = Tuple.T(3, Float.T(), optional=True)
1355 levels = List.T(CPTLevel.T())
1357 def scale(self, vmin, vmax):
1358 vmin_old, vmax_old = self.levels[0].vmin, self.levels[-1].vmax
1359 for level in self.levels:
1360 level.vmin = (level.vmin - vmin_old) / (vmax_old - vmin_old) * \
1361 (vmax - vmin) + vmin
1362 level.vmax = (level.vmax - vmin_old) / (vmax_old - vmin_old) * \
1363 (vmax - vmin) + vmin
1365 def discretize(self, nlevels):
1366 colors = []
1367 vals = []
1368 for level in self.levels:
1369 vals.append(level.vmin)
1370 vals.append(level.vmax)
1371 colors.append(level.color_min)
1372 colors.append(level.color_max)
1374 r, g, b = num.array(colors, dtype=float).T
1375 vals = num.array(vals, dtype=float)
1377 vmin, vmax = self.levels[0].vmin, self.levels[-1].vmax
1378 x = num.linspace(vmin, vmax, nlevels+1)
1379 rd = num.interp(x, vals, r)
1380 gd = num.interp(x, vals, g)
1381 bd = num.interp(x, vals, b)
1383 levels = []
1384 for ilevel in range(nlevels):
1385 color = (
1386 float(0.5*(rd[ilevel]+rd[ilevel+1])),
1387 float(0.5*(gd[ilevel]+gd[ilevel+1])),
1388 float(0.5*(bd[ilevel]+bd[ilevel+1])))
1390 levels.append(CPTLevel(
1391 vmin=x[ilevel],
1392 vmax=x[ilevel+1],
1393 color_min=color,
1394 color_max=color))
1396 cpt = CPT(
1397 color_below=self.color_below,
1398 color_above=self.color_above,
1399 color_nan=self.color_nan,
1400 levels=levels)
1402 return cpt
1405class CPTParseError(Exception):
1406 pass
1409def read_cpt(filename):
1410 with open(filename) as f:
1411 color_below = None
1412 color_above = None
1413 color_nan = None
1414 levels = []
1415 try:
1416 for line in f:
1417 line = line.strip()
1418 toks = line.split()
1420 if line.startswith('#'):
1421 continue
1423 elif line.startswith('B'):
1424 color_below = tuple(map(float, toks[1:4]))
1426 elif line.startswith('F'):
1427 color_above = tuple(map(float, toks[1:4]))
1429 elif line.startswith('N'):
1430 color_nan = tuple(map(float, toks[1:4]))
1432 else:
1433 values = list(map(float, line.split()))
1434 vmin = values[0]
1435 color_min = tuple(values[1:4])
1436 vmax = values[4]
1437 color_max = tuple(values[5:8])
1438 levels.append(CPTLevel(
1439 vmin=vmin,
1440 vmax=vmax,
1441 color_min=color_min,
1442 color_max=color_max))
1444 except Exception:
1445 raise CPTParseError()
1447 return CPT(
1448 color_below=color_below,
1449 color_above=color_above,
1450 color_nan=color_nan,
1451 levels=levels)
1454def color_to_int(color):
1455 return tuple(max(0, min(255, int(round(x)))) for x in color)
1458def write_cpt(cpt, filename):
1459 with open(filename, 'w') as f:
1460 for level in cpt.levels:
1461 f.write(
1462 '%e %i %i %i %e %i %i %i\n' %
1463 ((level.vmin, ) + color_to_int(level.color_min) +
1464 (level.vmax, ) + color_to_int(level.color_max)))
1466 if cpt.color_below:
1467 f.write('B %i %i %i\n' % color_to_int(cpt.color_below))
1469 if cpt.color_above:
1470 f.write('F %i %i %i\n' % color_to_int(cpt.color_above))
1472 if cpt.color_nan:
1473 f.write('N %i %i %i\n' % color_to_int(cpt.color_nan))
1476def cpt_merge_wet_dry(wet, dry):
1477 levels = []
1478 for level in wet.levels:
1479 if level.vmin < 0.:
1480 if level.vmax > 0.:
1481 level.vmax = 0.
1483 levels.append(level)
1485 for level in dry.levels:
1486 if level.vmax > 0.:
1487 if level.vmin < 0.:
1488 level.vmin = 0.
1490 levels.append(level)
1492 combi = CPT(
1493 color_below=wet.color_below,
1494 color_above=dry.color_above,
1495 color_nan=dry.color_nan,
1496 levels=levels)
1498 return combi
1501if __name__ == '__main__':
1502 from pyrocko import util
1503 util.setup_logging('pyrocko.automap', 'info')
1505 import sys
1506 if len(sys.argv) == 2:
1508 n = int(sys.argv[1])
1510 for i in range(n):
1511 m = Map(
1512 lat=rand(-60., 60.),
1513 lon=rand(-180., 180.),
1514 radius=math.exp(rand(math.log(500*km), math.log(3000*km))),
1515 width=30., height=30.,
1516 show_grid=True,
1517 show_topo=True,
1518 color_dry=(238, 236, 230),
1519 topo_cpt_wet='light_sea_uniform',
1520 topo_cpt_dry='light_land_uniform',
1521 illuminate=True,
1522 illuminate_factor_ocean=0.15,
1523 show_rivers=False,
1524 show_plates=True)
1526 m.draw_cities()
1527 print(m)
1528 m.save('map_%02i.pdf' % i)