1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
5from __future__ import absolute_import, print_function
7import math
8import random
9import logging
11try:
12 from StringIO import StringIO as BytesIO
13except ImportError:
14 from io import BytesIO
16import numpy as num
18from pyrocko.guts import (Object, Float, Bool, Int, Tuple, String, List,
19 Unicode, Dict)
20from pyrocko.guts_array import Array
21from pyrocko.dataset import topo
22from pyrocko import orthodrome as od
23from . import gmtpy
25try:
26 newstr = unicode
27except NameError:
28 newstr = str
30points_in_region = od.points_in_region
32logger = logging.getLogger('pyrocko.plot.automap')
34earthradius = 6371000.0
35r2d = 180./math.pi
36d2r = 1./r2d
37km = 1000.
38d2m = d2r*earthradius
39m2d = 1./d2m
40cm = gmtpy.cm
43def darken(c, f=0.7):
44 return (c[0]*f, c[1]*f, c[2]*f)
47def corners(lon, lat, w, h):
48 ll_lat, ll_lon = od.ne_to_latlon(lat, lon, -0.5*h, -0.5*w)
49 ur_lat, ur_lon = od.ne_to_latlon(lat, lon, 0.5*h, 0.5*w)
50 return ll_lon, ll_lat, ur_lon, ur_lat
53def extent(lon, lat, w, h, n):
54 x = num.linspace(-0.5*w, 0.5*w, n)
55 y = num.linspace(-0.5*h, 0.5*h, n)
56 slats, slons = od.ne_to_latlon(lat, lon, y[0], x)
57 nlats, nlons = od.ne_to_latlon(lat, lon, y[-1], x)
58 south = slats.min()
59 north = nlats.max()
61 wlats, wlons = od.ne_to_latlon(lat, lon, y, x[0])
62 elats, elons = od.ne_to_latlon(lat, lon, y, x[-1])
63 elons = num.where(elons < wlons, elons + 360., elons)
65 if elons.max() - elons.min() > 180 or wlons.max() - wlons.min() > 180.:
66 west = -180.
67 east = 180.
68 else:
69 west = wlons.min()
70 east = elons.max()
72 return topo.positive_region((west, east, south, north))
75class NoTopo(Exception):
76 pass
79class OutOfBounds(Exception):
80 pass
83class FloatTile(Object):
84 xmin = Float.T()
85 ymin = Float.T()
86 dx = Float.T()
87 dy = Float.T()
88 data = Array.T(shape=(None, None), dtype=float, serialize_as='table')
90 def __init__(self, xmin, ymin, dx, dy, data):
91 Object.__init__(self, init_props=False)
92 self.xmin = float(xmin)
93 self.ymin = float(ymin)
94 self.dx = float(dx)
95 self.dy = float(dy)
96 self.data = data
97 self._set_maxes()
99 def _set_maxes(self):
100 self.ny, self.nx = self.data.shape
101 self.xmax = self.xmin + (self.nx-1) * self.dx
102 self.ymax = self.ymin + (self.ny-1) * self.dy
104 def x(self):
105 return self.xmin + num.arange(self.nx) * self.dx
107 def y(self):
108 return self.ymin + num.arange(self.ny) * self.dy
110 def get(self, x, y):
111 ix = int(round((x - self.xmin) / self.dx))
112 iy = int(round((y - self.ymin) / self.dy))
113 if 0 <= ix < self.nx and 0 <= iy < self.ny:
114 return self.data[iy, ix]
115 else:
116 raise OutOfBounds()
119class City(Object):
120 def __init__(self, name, lat, lon, population=None, asciiname=None):
121 name = newstr(name)
122 lat = float(lat)
123 lon = float(lon)
124 if asciiname is None:
125 asciiname = name.encode('ascii', errors='replace')
127 if population is None:
128 population = 0
129 else:
130 population = int(population)
132 Object.__init__(self, name=name, lat=lat, lon=lon,
133 population=population, asciiname=asciiname)
135 name = Unicode.T()
136 lat = Float.T()
137 lon = Float.T()
138 population = Int.T()
139 asciiname = String.T()
142class Map(Object):
143 lat = Float.T(optional=True)
144 lon = Float.T(optional=True)
145 radius = Float.T(optional=True)
146 width = Float.T(default=20.)
147 height = Float.T(default=14.)
148 margins = List.T(Float.T())
149 illuminate = Bool.T(default=True)
150 skip_feature_factor = Float.T(default=0.02)
151 show_grid = Bool.T(default=False)
152 show_topo = Bool.T(default=True)
153 show_scale = Bool.T(default=False)
154 show_topo_scale = Bool.T(default=False)
155 show_center_mark = Bool.T(default=False)
156 show_rivers = Bool.T(default=True)
157 show_plates = Bool.T(default=False)
158 show_plate_velocities = Bool.T(default=False)
159 show_plate_names = Bool.T(default=False)
160 show_boundaries = Bool.T(default=False)
161 illuminate_factor_land = Float.T(default=0.5)
162 illuminate_factor_ocean = Float.T(default=0.25)
163 color_wet = Tuple.T(3, Int.T(), default=(216, 242, 254))
164 color_dry = Tuple.T(3, Int.T(), default=(172, 208, 165))
165 color_boundaries = Tuple.T(3, Int.T(), default=(1, 1, 1))
166 topo_resolution_min = Float.T(
167 default=40.,
168 help='minimum resolution of topography [dpi]')
169 topo_resolution_max = Float.T(
170 default=200.,
171 help='maximum resolution of topography [dpi]')
172 replace_topo_color_only = FloatTile.T(
173 optional=True,
174 help='replace topo color while keeping topographic shading')
175 topo_cpt_wet = String.T(default='light_sea')
176 topo_cpt_dry = String.T(default='light_land')
177 axes_layout = String.T(optional=True)
178 custom_cities = List.T(City.T())
179 gmt_config = Dict.T(String.T(), String.T())
180 comment = String.T(optional=True)
181 approx_ticks = Int.T(default=4)
183 def __init__(self, gmtversion='newest', **kwargs):
184 Object.__init__(self, **kwargs)
185 self._gmt = None
186 self._scaler = None
187 self._widget = None
188 self._corners = None
189 self._wesn = None
190 self._minarea = None
191 self._coastline_resolution = None
192 self._rivers = None
193 self._dems = None
194 self._have_topo_land = None
195 self._have_topo_ocean = None
196 self._jxyr = None
197 self._prep_topo_have = None
198 self._labels = []
199 self._area_labels = []
200 self._gmtversion = gmtversion
202 def save(self, outpath, resolution=75., oversample=2., size=None,
203 width=None, height=None, psconvert=False, crop_eps_mode=False):
205 '''
206 Save the image.
208 Save the image to ``outpath``. The format is determined by the filename
209 extension. Formats are handled as follows: ``'.eps'`` and ``'.ps'``
210 produce EPS and PS, respectively, directly with GMT. If the file name
211 ends with ``'.pdf'``, GMT output is fed through ``gmtpy-epstopdf`` to
212 create a PDF file. For any other filename extension, output is first
213 converted to PDF with ``gmtpy-epstopdf``, then with ``pdftocairo`` to
214 PNG with a resolution oversampled by the factor ``oversample`` and
215 finally the PNG is downsampled and converted to the target format with
216 ``convert``. The resolution of rasterized target image can be
217 controlled either by ``resolution`` in DPI or by specifying ``width``
218 or ``height`` or ``size``, where the latter fits the image into a
219 square with given side length. To save transparency use
220 ``psconvert=True``. To crop the output image with a rectangle to the
221 nearest non-white element set ``crop_eps_mode=True``.
222 '''
224 gmt = self.gmt
225 self.draw_labels()
226 self.draw_axes()
227 if self.show_topo and self.show_topo_scale:
228 self._draw_topo_scale()
230 gmt.save(outpath, resolution=resolution, oversample=oversample,
231 crop_eps_mode=crop_eps_mode,
232 size=size, width=width, height=height, psconvert=psconvert)
234 @property
235 def scaler(self):
236 if self._scaler is None:
237 self._setup_geometry()
239 return self._scaler
241 @property
242 def wesn(self):
243 if self._wesn is None:
244 self._setup_geometry()
246 return self._wesn
248 @property
249 def widget(self):
250 if self._widget is None:
251 self._setup()
253 return self._widget
255 @property
256 def layout(self):
257 if self._layout is None:
258 self._setup()
260 return self._layout
262 @property
263 def jxyr(self):
264 if self._jxyr is None:
265 self._setup()
267 return self._jxyr
269 @property
270 def pxyr(self):
271 if self._pxyr is None:
272 self._setup()
274 return self._pxyr
276 @property
277 def gmt(self):
278 if self._gmt is None:
279 self._setup()
281 if self._have_topo_ocean is None:
282 self._draw_background()
284 return self._gmt
286 def _setup(self):
287 if not self._widget:
288 self._setup_geometry()
290 self._setup_lod()
291 self._setup_gmt()
293 def _setup_geometry(self):
294 wpage, hpage = self.width, self.height
295 ml, mr, mt, mb = self._expand_margins()
296 wpage -= ml + mr
297 hpage -= mt + mb
299 wreg = self.radius * 2.0
300 hreg = self.radius * 2.0
301 if wpage >= hpage:
302 wreg *= wpage/hpage
303 else:
304 hreg *= hpage/wpage
306 self._wreg = wreg
307 self._hreg = hreg
309 self._corners = corners(self.lon, self.lat, wreg, hreg)
310 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10)
312 x, y, z = ((west, east), (south, north), (-6000., 4500.))
314 xax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
315 yax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
316 zax = gmtpy.Ax(mode='min-max', inc=1000., label='Height',
317 scaled_unit='km', scaled_unit_factor=0.001)
319 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax))
321 par = scaler.get_params()
323 west = par['xmin']
324 east = par['xmax']
325 south = par['ymin']
326 north = par['ymax']
328 self._wesn = west, east, south, north
329 self._scaler = scaler
331 def _setup_lod(self):
332 w, e, s, n = self._wesn
333 if self.radius > 1500.*km:
334 coastline_resolution = 'i'
335 rivers = False
336 else:
337 coastline_resolution = 'f'
338 rivers = True
340 self._minarea = (self.skip_feature_factor * self.radius/km)**2
342 self._coastline_resolution = coastline_resolution
343 self._rivers = rivers
345 self._prep_topo_have = {}
346 self._dems = {}
348 cm2inch = gmtpy.cm/gmtpy.inch
350 dmin = 2.0 * self.radius * m2d / (self.topo_resolution_max *
351 (self.height * cm2inch))
352 dmax = 2.0 * self.radius * m2d / (self.topo_resolution_min *
353 (self.height * cm2inch))
355 for k in ['ocean', 'land']:
356 self._dems[k] = topo.select_dem_names(k, dmin, dmax, self._wesn)
357 if self._dems[k]:
358 logger.debug('using topography dataset %s for %s'
359 % (','.join(self._dems[k]), k))
361 def _expand_margins(self):
362 if len(self.margins) == 0 or len(self.margins) > 4:
363 ml = mr = mt = mb = 2.0
364 elif len(self.margins) == 1:
365 ml = mr = mt = mb = self.margins[0]
366 elif len(self.margins) == 2:
367 ml = mr = self.margins[0]
368 mt = mb = self.margins[1]
369 elif len(self.margins) == 4:
370 ml, mr, mt, mb = self.margins
372 return ml, mr, mt, mb
374 def _setup_gmt(self):
375 w, h = self.width, self.height
376 scaler = self._scaler
378 if gmtpy.is_gmt5(self._gmtversion):
379 gmtconf = dict(
380 MAP_TICK_PEN_PRIMARY='1.25p',
381 MAP_TICK_PEN_SECONDARY='1.25p',
382 MAP_TICK_LENGTH_PRIMARY='0.2c',
383 MAP_TICK_LENGTH_SECONDARY='0.6c',
384 FONT_ANNOT_PRIMARY='12p,1,black',
385 FONT_LABEL='12p,1,black',
386 PS_CHAR_ENCODING='ISOLatin1+',
387 MAP_FRAME_TYPE='fancy',
388 FORMAT_GEO_MAP='D',
389 PS_MEDIA='Custom_%ix%i' % (
390 w*gmtpy.cm,
391 h*gmtpy.cm),
392 PS_PAGE_ORIENTATION='portrait',
393 MAP_GRID_PEN_PRIMARY='thinnest,0/50/0',
394 MAP_ANNOT_OBLIQUE='6')
395 else:
396 gmtconf = dict(
397 TICK_PEN='1.25p',
398 TICK_LENGTH='0.2c',
399 ANNOT_FONT_PRIMARY='1',
400 ANNOT_FONT_SIZE_PRIMARY='12p',
401 LABEL_FONT='1',
402 LABEL_FONT_SIZE='12p',
403 CHAR_ENCODING='ISOLatin1+',
404 BASEMAP_TYPE='fancy',
405 PLOT_DEGREE_FORMAT='D',
406 PAPER_MEDIA='Custom_%ix%i' % (
407 w*gmtpy.cm,
408 h*gmtpy.cm),
409 GRID_PEN_PRIMARY='thinnest/0/50/0',
410 DOTS_PR_INCH='1200',
411 OBLIQUE_ANNOTATION='6')
413 gmtconf.update(
414 (k.upper(), v) for (k, v) in self.gmt_config.items())
416 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion)
418 layout = gmt.default_layout()
420 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()])
422 widget = layout.get_widget()
423 widget['P'] = widget['J']
424 widget['J'] = ('-JA%g/%g' % (self.lon, self.lat)) + '/%(width)gp'
425 scaler['R'] = '-R%g/%g/%g/%gr' % self._corners
427 # aspect = gmtpy.aspect_for_projection(
428 # gmt.installation['version'], *(widget.J() + scaler.R()))
430 aspect = self._map_aspect(jr=widget.J() + scaler.R())
431 widget.set_aspect(aspect)
433 self._gmt = gmt
434 self._layout = layout
435 self._widget = widget
436 self._jxyr = self._widget.JXY() + self._scaler.R()
437 self._pxyr = self._widget.PXY() + [
438 '-R%g/%g/%g/%g' % (0, widget.width(), 0, widget.height())]
439 self._have_drawn_axes = False
440 self._have_drawn_labels = False
442 def _draw_background(self):
443 self._have_topo_land = False
444 self._have_topo_ocean = False
445 if self.show_topo:
446 self._have_topo = self._draw_topo()
448 self._draw_basefeatures()
450 def _get_topo_tile(self, k):
451 t = None
452 demname = None
453 for dem in self._dems[k]:
454 t = topo.get(dem, self._wesn)
455 demname = dem
456 if t is not None:
457 break
459 if not t:
460 raise NoTopo()
462 return t, demname
464 def _prep_topo(self, k):
465 gmt = self._gmt
466 t, demname = self._get_topo_tile(k)
468 if demname not in self._prep_topo_have:
470 grdfile = gmt.tempfilename()
472 is_flat = num.all(t.data[0] == t.data)
474 gmtpy.savegrd(
475 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat')
477 if self.illuminate and not is_flat:
478 if k == 'ocean':
479 factor = self.illuminate_factor_ocean
480 else:
481 factor = self.illuminate_factor_land
483 ilumfn = gmt.tempfilename()
484 gmt.grdgradient(
485 grdfile,
486 N='e%g' % factor,
487 A=-45,
488 G=ilumfn,
489 out_discard=True)
491 ilumargs = ['-I%s' % ilumfn]
492 else:
493 ilumargs = []
495 if self.replace_topo_color_only:
496 t2 = self.replace_topo_color_only
497 grdfile2 = gmt.tempfilename()
499 gmtpy.savegrd(
500 t2.x(), t2.y(), t2.data, filename=grdfile2,
501 naming='lonlat')
503 if gmt.is_gmt5():
504 gmt.grdsample(
505 grdfile2,
506 G=grdfile,
507 n='l',
508 I='%g/%g' % (t.dx, t.dy), # noqa
509 R=grdfile,
510 out_discard=True)
511 else:
512 gmt.grdsample(
513 grdfile2,
514 G=grdfile,
515 Q='l',
516 I='%g/%g' % (t.dx, t.dy), # noqa
517 R=grdfile,
518 out_discard=True)
520 gmt.grdmath(
521 grdfile, '0.0', 'AND', '=', grdfile2,
522 out_discard=True)
524 grdfile = grdfile2
526 self._prep_topo_have[demname] = grdfile, ilumargs
528 return self._prep_topo_have[demname]
530 def _draw_topo(self):
531 widget = self._widget
532 scaler = self._scaler
533 gmt = self._gmt
534 cres = self._coastline_resolution
535 minarea = self._minarea
537 JXY = widget.JXY()
538 R = scaler.R()
540 try:
541 grdfile, ilumargs = self._prep_topo('ocean')
542 gmt.pscoast(D=cres, S='c', A=minarea, *(JXY+R))
543 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_wet),
544 *(ilumargs+JXY+R))
545 gmt.pscoast(Q=True, *(JXY+R))
546 self._have_topo_ocean = True
547 except NoTopo:
548 self._have_topo_ocean = False
550 try:
551 grdfile, ilumargs = self._prep_topo('land')
552 gmt.pscoast(D=cres, G='c', A=minarea, *(JXY+R))
553 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_dry),
554 *(ilumargs+JXY+R))
555 gmt.pscoast(Q=True, *(JXY+R))
556 self._have_topo_land = True
557 except NoTopo:
558 self._have_topo_land = False
560 def _draw_topo_scale(self, label='Elevation [km]'):
561 dry = read_cpt(topo.cpt(self.topo_cpt_dry))
562 wet = read_cpt(topo.cpt(self.topo_cpt_wet))
563 combi = cpt_merge_wet_dry(wet, dry)
564 for level in combi.levels:
565 level.vmin /= km
566 level.vmax /= km
568 topo_cpt = self.gmt.tempfilename() + '.cpt'
569 write_cpt(combi, topo_cpt)
571 (w, h), (xo, yo) = self.widget.get_size()
572 self.gmt.psscale(
573 D='%gp/%gp/%gp/%gph' % (xo + 0.5*w, yo - 2.0*gmtpy.cm, w,
574 0.5*gmtpy.cm),
575 C=topo_cpt,
576 B='1:%s:' % label)
578 def _draw_basefeatures(self):
579 gmt = self._gmt
580 cres = self._coastline_resolution
581 rivers = self._rivers
582 minarea = self._minarea
584 color_wet = self.color_wet
585 color_dry = self.color_dry
587 if self.show_rivers and rivers:
588 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)]
589 else:
590 rivers = []
592 fill = {}
593 if not self._have_topo_land:
594 fill['G'] = color_dry
596 if not self._have_topo_ocean:
597 fill['S'] = color_wet
599 if self.show_boundaries:
600 fill['N'] = '1/1p,%s,%s' % (
601 gmtpy.color(self.color_boundaries), 'solid')
603 gmt.pscoast(
604 D=cres,
605 W='thinnest,%s' % gmtpy.color(darken(gmtpy.color_tup(color_dry))),
606 A=minarea,
607 *(rivers+self._jxyr), **fill)
609 if self.show_plates:
610 self.draw_plates()
612 def _draw_axes(self):
613 gmt = self._gmt
614 scaler = self._scaler
615 widget = self._widget
617 if self.axes_layout is None:
618 if self.lat > 0.0:
619 axes_layout = 'WSen'
620 else:
621 axes_layout = 'WseN'
622 else:
623 axes_layout = self.axes_layout
625 scale_km = gmtpy.nice_value(self.radius/5.) / 1000.
627 if self.show_center_mark:
628 gmt.psxy(
629 in_rows=[[self.lon, self.lat]],
630 S='c20p', W='2p,black',
631 *self._jxyr)
633 if self.show_grid:
634 btmpl = ('%(xinc)gg%(xinc)g:%(xlabel)s:/'
635 '%(yinc)gg%(yinc)g:%(ylabel)s:')
636 else:
637 btmpl = '%(xinc)g:%(xlabel)s:/%(yinc)g:%(ylabel)s:'
639 if self.show_scale:
640 scale = 'x%gp/%gp/%g/%g/%gk' % (
641 6./7*widget.width(),
642 widget.height()/7.,
643 self.lon,
644 self.lat,
645 scale_km)
646 else:
647 scale = False
649 gmt.psbasemap(
650 B=(btmpl % scaler.get_params())+axes_layout,
651 L=scale,
652 *self._jxyr)
654 if self.comment:
655 font_size = self.gmt.label_font_size()
657 _, east, south, _ = self._wesn
658 if gmt.is_gmt5():
659 row = [
660 1, 0,
661 '%gp,%s,%s' % (font_size, 0, 'black'), 'BR',
662 self.comment]
664 farg = ['-F+f+j']
665 else:
666 row = [1, 0, font_size, 0, 0, 'BR', self.comment]
667 farg = []
669 gmt.pstext(
670 in_rows=[row],
671 N=True,
672 R=(0, 1, 0, 1),
673 D='%gp/%gp' % (-font_size*0.2, font_size*0.3),
674 *(widget.PXY() + farg))
676 def draw_axes(self):
677 if not self._have_drawn_axes:
678 self._draw_axes()
679 self._have_drawn_axes = True
681 def _have_coastlines(self):
682 gmt = self._gmt
683 cres = self._coastline_resolution
684 minarea = self._minarea
686 checkfile = gmt.tempfilename()
688 gmt.pscoast(
689 M=True,
690 D=cres,
691 W='thinnest,black',
692 A=minarea,
693 out_filename=checkfile,
694 *self._jxyr)
696 points = []
697 with open(checkfile, 'r') as f:
698 for line in f:
699 ls = line.strip()
700 if ls.startswith('#') or ls.startswith('>') or ls == '':
701 continue
702 plon, plat = [float(x) for x in ls.split()]
703 points.append((plat, plon))
705 points = num.array(points, dtype=float)
706 return num.any(points_in_region(points, self._wesn))
708 def have_coastlines(self):
709 self.gmt
710 return self._have_coastlines()
712 def project(self, lats, lons, jr=None):
713 onepoint = False
714 if isinstance(lats, float) and isinstance(lons, float):
715 lats = [lats]
716 lons = [lons]
717 onepoint = True
719 if jr is not None:
720 j, r = jr
721 gmt = gmtpy.GMT(version=self._gmtversion)
722 else:
723 j, _, _, r = self.jxyr
724 gmt = self.gmt
726 f = BytesIO()
727 gmt.mapproject(j, r, in_columns=(lons, lats), out_stream=f, D='p')
728 f.seek(0)
729 data = num.loadtxt(f, ndmin=2)
730 xs, ys = data.T
731 if onepoint:
732 xs = xs[0]
733 ys = ys[0]
734 return xs, ys
736 def _map_box(self, jr=None):
737 ll_lon, ll_lat, ur_lon, ur_lat = self._corners
739 xs_corner, ys_corner = self.project(
740 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr)
742 w = xs_corner[1] - xs_corner[0]
743 h = ys_corner[1] - ys_corner[0]
745 return w, h
747 def _map_aspect(self, jr=None):
748 w, h = self._map_box(jr=jr)
749 return h/w
751 def _draw_labels(self):
752 points_taken = []
753 regions_taken = []
755 def no_points_in_rect(xs, ys, xmin, ymin, xmax, ymax):
756 xx = not num.any(la(la(xmin < xs, xs < xmax),
757 la(ymin < ys, ys < ymax)))
758 return xx
760 def roverlaps(a, b):
761 return (a[0] < b[2] and b[0] < a[2] and
762 a[1] < b[3] and b[1] < a[3])
764 w, h = self._map_box()
766 label_font_size = self.gmt.label_font_size()
768 if self._labels:
770 n = len(self._labels)
772 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \
773 angles, styles = list(zip(*self._labels))
775 font_sizes = [
776 (font_size or label_font_size) for font_size in font_sizes]
778 sx = num.array(sx, dtype=float)
779 sy = num.array(sy, dtype=float)
781 xs, ys = self.project(lats, lons)
783 points_taken.append((xs, ys))
785 dxs = num.zeros(n)
786 dys = num.zeros(n)
788 for i in range(n):
789 dx, dy = gmtpy.text_box(
790 texts[i],
791 font=fonts[i],
792 font_size=font_sizes[i],
793 **styles[i])
795 dxs[i] = dx
796 dys[i] = dy
798 la = num.logical_and
799 anchors_ok = (
800 la(xs + sx + dxs < w, ys + sy + dys < h),
801 la(xs - sx - dxs > 0., ys - sy - dys > 0.),
802 la(xs + sx + dxs < w, ys - sy - dys > 0.),
803 la(xs - sx - dxs > 0., ys + sy + dys < h),
804 )
806 arects = [
807 (xs, ys, xs + sx + dxs, ys + sy + dys),
808 (xs - sx - dxs, ys - sy - dys, xs, ys),
809 (xs, ys - sy - dys, xs + sx + dxs, ys),
810 (xs - sx - dxs, ys, xs, ys + sy + dys)]
812 for i in range(n):
813 for ianch in range(4):
814 anchors_ok[ianch][i] &= no_points_in_rect(
815 xs, ys, *[xxx[i] for xxx in arects[ianch]])
817 anchor_choices = []
818 anchor_take = []
819 for i in range(n):
820 choices = [ianch for ianch in range(4)
821 if anchors_ok[ianch][i]]
822 anchor_choices.append(choices)
823 if choices:
824 anchor_take.append(choices[0])
825 else:
826 anchor_take.append(None)
828 def cost(anchor_take):
829 noverlaps = 0
830 for i in range(n):
831 for j in range(n):
832 if i != j:
833 i_take = anchor_take[i]
834 j_take = anchor_take[j]
835 if i_take is None or j_take is None:
836 continue
837 r_i = [xxx[i] for xxx in arects[i_take]]
838 r_j = [xxx[j] for xxx in arects[j_take]]
839 if roverlaps(r_i, r_j):
840 noverlaps += 1
842 return noverlaps
844 cur_cost = cost(anchor_take)
845 imax = 30
846 while cur_cost != 0 and imax > 0:
847 for i in range(n):
848 for t in anchor_choices[i]:
849 anchor_take_new = list(anchor_take)
850 anchor_take_new[i] = t
851 new_cost = cost(anchor_take_new)
852 if new_cost < cur_cost:
853 anchor_take = anchor_take_new
854 cur_cost = new_cost
856 imax -= 1
858 while cur_cost != 0:
859 for i in range(n):
860 anchor_take_new = list(anchor_take)
861 anchor_take_new[i] = None
862 new_cost = cost(anchor_take_new)
863 if new_cost < cur_cost:
864 anchor_take = anchor_take_new
865 cur_cost = new_cost
866 break
868 anchor_strs = ['BL', 'TR', 'TL', 'BR']
870 for i in range(n):
871 ianchor = anchor_take[i]
872 color = colors[i]
873 if color is None:
874 color = 'black'
876 if ianchor is not None:
877 regions_taken.append([xxx[i] for xxx in arects[ianchor]])
879 anchor = anchor_strs[ianchor]
881 yoff = [-sy[i], sy[i]][anchor[0] == 'B']
882 xoff = [-sx[i], sx[i]][anchor[1] == 'L']
883 if self.gmt.is_gmt5():
884 row = (
885 lons[i], lats[i],
886 '%i,%s,%s' % (font_sizes[i], fonts[i], color),
887 anchor,
888 texts[i])
890 farg = ['-F+f+j+a%g' % angles[i]]
891 else:
892 row = (
893 lons[i], lats[i],
894 font_sizes[i], angles[i], fonts[i], anchor,
895 texts[i])
896 farg = ['-G%s' % color]
898 self.gmt.pstext(
899 in_rows=[row],
900 D='%gp/%gp' % (xoff, yoff),
901 *(self.jxyr + farg),
902 **styles[i])
904 if self._area_labels:
906 for lons, lats, text, color, font, font_size, style in \
907 self._area_labels:
909 if font_size is None:
910 font_size = label_font_size
912 if color is None:
913 color = 'black'
915 if self.gmt.is_gmt5():
916 farg = ['-F+f+j']
917 else:
918 farg = ['-G%s' % color]
920 xs, ys = self.project(lats, lons)
921 dx, dy = gmtpy.text_box(
922 text, font=font, font_size=font_size, **style)
924 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy]
926 locs_ok = num.ones(xs.size, dtype=num.bool)
928 for iloc in range(xs.size):
929 rcandi = [xxx[iloc] for xxx in rects]
931 locs_ok[iloc] = True
932 locs_ok[iloc] &= (
933 0 < rcandi[0] and rcandi[2] < w
934 and 0 < rcandi[1] and rcandi[3] < h)
936 overlap = False
937 for r in regions_taken:
938 if roverlaps(r, rcandi):
939 overlap = True
940 break
942 locs_ok[iloc] &= not overlap
944 for xs_taken, ys_taken in points_taken:
945 locs_ok[iloc] &= no_points_in_rect(
946 xs_taken, ys_taken, *rcandi)
948 if not locs_ok[iloc]:
949 break
951 rows = []
952 for iloc, (lon, lat) in enumerate(zip(lons, lats)):
953 if not locs_ok[iloc]:
954 continue
956 if self.gmt.is_gmt5():
957 row = (
958 lon, lat,
959 '%i,%s,%s' % (font_size, font, color),
960 'MC',
961 text)
963 else:
964 row = (
965 lon, lat,
966 font_size, 0, font, 'MC',
967 text)
969 rows.append(row)
971 regions_taken.append([xxx[iloc] for xxx in rects])
972 break
974 self.gmt.pstext(
975 in_rows=rows,
976 *(self.jxyr + farg),
977 **style)
979 def draw_labels(self):
980 self.gmt
981 if not self._have_drawn_labels:
982 self._draw_labels()
983 self._have_drawn_labels = True
985 def add_label(
986 self, lat, lon, text,
987 offset_x=5., offset_y=5.,
988 color=None,
989 font='1',
990 font_size=None,
991 angle=0,
992 style={}):
994 if 'G' in style:
995 style = style.copy()
996 color = style.pop('G')
998 self._labels.append(
999 (lon, lat, text, offset_x, offset_y, color, font, font_size,
1000 angle, style))
1002 def add_area_label(
1003 self, lat, lon, text,
1004 color=None,
1005 font='3',
1006 font_size=None,
1007 style={}):
1009 self._area_labels.append(
1010 (lon, lat, text, color, font, font_size, style))
1012 def cities_in_region(self):
1013 from pyrocko.dataset import geonames
1014 cities = geonames.get_cities_region(region=self.wesn, minpop=0)
1015 cities.extend(self.custom_cities)
1016 cities.sort(key=lambda x: x.population)
1017 return cities
1019 def draw_cities(self,
1020 exact=None,
1021 include=[],
1022 exclude=[],
1023 nmax_soft=10,
1024 psxy_style=dict(S='s5p', G='black')):
1026 cities = self.cities_in_region()
1028 if exact is not None:
1029 cities = [c for c in cities if c.name in exact]
1030 minpop = None
1031 else:
1032 cities = [c for c in cities if c.name not in exclude]
1033 minpop = 10**3
1034 for minpop_new in [1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6, 3e6, 1e7]:
1035 cities_new = [
1036 c for c in cities
1037 if c.population > minpop_new or c.name in include]
1039 if len(cities_new) == 0 or (
1040 len(cities_new) < 3 and len(cities) < nmax_soft*2):
1041 break
1043 cities = cities_new
1044 minpop = minpop_new
1045 if len(cities) <= nmax_soft:
1046 break
1048 if cities:
1049 lats = [c.lat for c in cities]
1050 lons = [c.lon for c in cities]
1052 self.gmt.psxy(
1053 in_columns=(lons, lats),
1054 *self.jxyr, **psxy_style)
1056 for c in cities:
1057 try:
1058 text = c.name.encode('iso-8859-1').decode('iso-8859-1')
1059 except UnicodeEncodeError:
1060 text = c.asciiname
1062 self.add_label(c.lat, c.lon, text)
1064 self._cities_minpop = minpop
1066 def add_stations(self, stations, psxy_style=dict()):
1068 default_psxy_style = {
1069 'S': 't8p',
1070 'G': 'black'
1071 }
1072 default_psxy_style.update(psxy_style)
1074 lats, lons = zip(*[s.effective_latlon for s in stations])
1076 self.gmt.psxy(
1077 in_columns=(lons, lats),
1078 *self.jxyr, **default_psxy_style)
1080 for station in stations:
1081 self.add_label(
1082 station.effective_lat,
1083 station.effective_lon,
1084 '.'.join(x for x in (station.network, station.station) if x))
1086 def add_kite_scene(self, scene):
1087 tile = FloatTile(
1088 scene.frame.llLon,
1089 scene.frame.llLat,
1090 scene.frame.dLon,
1091 scene.frame.dLat,
1092 scene.displacement)
1094 return tile
1096 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None,
1097 labels=True, vertical=False, fontsize=10):
1099 stations = campaign.stations
1101 if offset_scale is None:
1102 offset_scale = num.zeros(campaign.nstations)
1103 for ista, sta in enumerate(stations):
1104 for comp in sta.components.values():
1105 offset_scale[ista] += comp.shift
1106 offset_scale = num.sqrt(offset_scale**2).max()
1108 size = math.sqrt(self.height**2 + self.width**2)
1109 scale = (size/10.) / offset_scale
1110 logger.debug('GNSS: Using offset scale %f, map scale %f',
1111 offset_scale, scale)
1113 lats, lons = zip(*[s.effective_latlon for s in stations])
1115 if vertical:
1116 rows = [[lons[ista], lats[ista],
1117 0., -s.up.shift,
1118 (s.east.sigma + s.north.sigma) if s.east.sigma else 0.,
1119 s.up.sigma, 0.,
1120 s.code if labels else None]
1121 for ista, s in enumerate(stations)
1122 if s.up is not None]
1124 else:
1125 rows = [[lons[ista], lats[ista],
1126 -s.east.shift, -s.north.shift,
1127 s.east.sigma, s.north.sigma, s.correlation_ne,
1128 s.code if labels else None]
1129 for ista, s in enumerate(stations)
1130 if s.east is not None or s.north is not None]
1132 default_psxy_style = {
1133 'h': 0,
1134 'W': '2p,black',
1135 'A': '+p2p,black+b+a40',
1136 'G': 'black',
1137 'L': True,
1138 'S': 'e%dc/0.95/%d' % (scale, fontsize),
1139 }
1141 if not labels:
1142 for row in rows:
1143 row.pop(-1)
1145 if psxy_style is not None:
1146 default_psxy_style.update(psxy_style)
1148 self.gmt.psvelo(
1149 in_rows=rows,
1150 *self.jxyr,
1151 **default_psxy_style)
1153 def draw_plates(self):
1154 from pyrocko.dataset import tectonics
1156 neast = 20
1157 nnorth = max(1, int(round(num.round(self._hreg/self._wreg * neast))))
1158 norths = num.linspace(-self._hreg*0.5, self._hreg*0.5, nnorth)
1159 easts = num.linspace(-self._wreg*0.5, self._wreg*0.5, neast)
1160 norths2 = num.repeat(norths, neast)
1161 easts2 = num.tile(easts, nnorth)
1162 lats, lons = od.ne_to_latlon(
1163 self.lat, self.lon, norths2, easts2)
1165 bird = tectonics.PeterBird2003()
1166 plates = bird.get_plates()
1168 color_plates = gmtpy.color('aluminium5')
1169 color_velocities = gmtpy.color('skyblue1')
1170 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1')))
1172 points = num.vstack((lats, lons)).T
1173 used = []
1174 for plate in plates:
1175 mask = plate.contains_points(points)
1176 if num.any(mask):
1177 used.append((plate, mask))
1179 if len(used) > 1:
1181 candi_fixed = {}
1183 label_data = []
1184 for plate, mask in used:
1186 mean_north = num.mean(norths2[mask])
1187 mean_east = num.mean(easts2[mask])
1188 iorder = num.argsort(num.sqrt(
1189 (norths2[mask] - mean_north)**2 +
1190 (easts2[mask] - mean_east)**2))
1192 lat_candis = lats[mask][iorder]
1193 lon_candis = lons[mask][iorder]
1195 candi_fixed[plate.name] = lat_candis.size
1197 label_data.append((
1198 lat_candis, lon_candis, plate, color_plates))
1200 boundaries = bird.get_boundaries()
1202 size = 1.
1204 psxy_kwargs = []
1206 for boundary in boundaries:
1207 if num.any(points_in_region(boundary.points, self._wesn)):
1208 for typ, part in boundary.split_types(
1209 [['SUB'],
1210 ['OSR', 'OTF', 'OCB', 'CTF', 'CCB', 'CRB']]):
1212 lats, lons = part.T
1214 kwargs = {}
1215 if typ[0] == 'SUB':
1216 if boundary.kind == '\\':
1217 kwargs['S'] = 'f%g/%gp+t+r' % (
1218 0.45*size, 3.*size)
1219 elif boundary.kind == '/':
1220 kwargs['S'] = 'f%g/%gp+t+l' % (
1221 0.45*size, 3.*size)
1223 kwargs['G'] = color_plates
1225 kwargs['in_columns'] = (lons, lats)
1226 kwargs['W'] = '%gp,%s' % (size, color_plates),
1228 psxy_kwargs.append(kwargs)
1230 if boundary.kind == '\\':
1231 if boundary.plate_name2 in candi_fixed:
1232 candi_fixed[boundary.plate_name2] += \
1233 neast*nnorth
1235 elif boundary.kind == '/':
1236 if boundary.plate_name1 in candi_fixed:
1237 candi_fixed[boundary.plate_name1] += \
1238 neast*nnorth
1240 candi_fixed = [name for name in sorted(
1241 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])]
1243 candi_fixed.append(None)
1245 gsrm = tectonics.GSRM1()
1247 for name in candi_fixed:
1248 if name not in gsrm.plate_names() \
1249 and name not in gsrm.plate_alt_names():
1251 continue
1253 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \
1254 gsrm.get_velocities(name, region=self._wesn)
1256 fixed_plate_name = name
1258 if self.show_plate_velocities:
1259 self.gmt.psvelo(
1260 in_columns=(
1261 lons, lats, veast, vnorth, veast_err, vnorth_err,
1262 corr),
1263 W='0.25p,%s' % color_velocities,
1264 A='9p+e+g%s' % color_velocities,
1265 S='e0.2p/0.95/10',
1266 *self.jxyr)
1268 for _ in range(len(lons) // 50 + 1):
1269 ii = random.randint(0, len(lons)-1)
1270 v = math.sqrt(vnorth[ii]**2 + veast[ii]**2)
1271 self.add_label(
1272 lats[ii], lons[ii], '%.0f' % v,
1273 font_size=0.7*self.gmt.label_font_size(),
1274 style=dict(
1275 G=color_velocities_lab))
1277 break
1279 if self.show_plate_names:
1280 for (lat_candis, lon_candis, plate, color) in label_data:
1281 full_name = bird.full_name(plate.name)
1282 if plate.name == fixed_plate_name:
1283 full_name = '@_' + full_name + '@_'
1285 self.add_area_label(
1286 lat_candis, lon_candis,
1287 full_name,
1288 color=color,
1289 font='3')
1291 for kwargs in psxy_kwargs:
1292 self.gmt.psxy(*self.jxyr, **kwargs)
1295def rand(mi, ma):
1296 mi = float(mi)
1297 ma = float(ma)
1298 return random.random() * (ma-mi) + mi
1301def split_region(region):
1302 west, east, south, north = topo.positive_region(region)
1303 if east > 180:
1304 return [(west, 180., south, north),
1305 (-180., east-360., south, north)]
1306 else:
1307 return [region]
1310class CPTLevel(Object):
1311 vmin = Float.T()
1312 vmax = Float.T()
1313 color_min = Tuple.T(3, Float.T())
1314 color_max = Tuple.T(3, Float.T())
1317class CPT(Object):
1318 color_below = Tuple.T(3, Float.T(), optional=True)
1319 color_above = Tuple.T(3, Float.T(), optional=True)
1320 color_nan = Tuple.T(3, Float.T(), optional=True)
1321 levels = List.T(CPTLevel.T())
1323 def scale(self, vmin, vmax):
1324 vmin_old, vmax_old = self.levels[0].vmin, self.levels[-1].vmax
1325 for level in self.levels:
1326 level.vmin = (level.vmin - vmin_old) / (vmax_old - vmin_old) * \
1327 (vmax - vmin) + vmin
1328 level.vmax = (level.vmax - vmin_old) / (vmax_old - vmin_old) * \
1329 (vmax - vmin) + vmin
1331 def discretize(self, nlevels):
1332 colors = []
1333 vals = []
1334 for level in self.levels:
1335 vals.append(level.vmin)
1336 vals.append(level.vmax)
1337 colors.append(level.color_min)
1338 colors.append(level.color_max)
1340 r, g, b = num.array(colors, dtype=float).T
1341 vals = num.array(vals, dtype=float)
1343 vmin, vmax = self.levels[0].vmin, self.levels[-1].vmax
1344 x = num.linspace(vmin, vmax, nlevels+1)
1345 rd = num.interp(x, vals, r)
1346 gd = num.interp(x, vals, g)
1347 bd = num.interp(x, vals, b)
1349 levels = []
1350 for ilevel in range(nlevels):
1351 color = (
1352 float(0.5*(rd[ilevel]+rd[ilevel+1])),
1353 float(0.5*(gd[ilevel]+gd[ilevel+1])),
1354 float(0.5*(bd[ilevel]+bd[ilevel+1])))
1356 levels.append(CPTLevel(
1357 vmin=x[ilevel],
1358 vmax=x[ilevel+1],
1359 color_min=color,
1360 color_max=color))
1362 cpt = CPT(
1363 color_below=self.color_below,
1364 color_above=self.color_above,
1365 color_nan=self.color_nan,
1366 levels=levels)
1368 return cpt
1371class CPTParseError(Exception):
1372 pass
1375def read_cpt(filename):
1376 with open(filename) as f:
1377 color_below = None
1378 color_above = None
1379 color_nan = None
1380 levels = []
1381 try:
1382 for line in f:
1383 line = line.strip()
1384 toks = line.split()
1386 if line.startswith('#'):
1387 continue
1389 elif line.startswith('B'):
1390 color_below = tuple(map(float, toks[1:4]))
1392 elif line.startswith('F'):
1393 color_above = tuple(map(float, toks[1:4]))
1395 elif line.startswith('N'):
1396 color_nan = tuple(map(float, toks[1:4]))
1398 else:
1399 values = list(map(float, line.split()))
1400 vmin = values[0]
1401 color_min = tuple(values[1:4])
1402 vmax = values[4]
1403 color_max = tuple(values[5:8])
1404 levels.append(CPTLevel(
1405 vmin=vmin,
1406 vmax=vmax,
1407 color_min=color_min,
1408 color_max=color_max))
1410 except Exception:
1411 raise CPTParseError()
1413 return CPT(
1414 color_below=color_below,
1415 color_above=color_above,
1416 color_nan=color_nan,
1417 levels=levels)
1420def color_to_int(color):
1421 return tuple(max(0, min(255, int(round(x)))) for x in color)
1424def write_cpt(cpt, filename):
1425 with open(filename, 'w') as f:
1426 for level in cpt.levels:
1427 f.write(
1428 '%e %i %i %i %e %i %i %i\n' %
1429 ((level.vmin, ) + color_to_int(level.color_min) +
1430 (level.vmax, ) + color_to_int(level.color_max)))
1432 if cpt.color_below:
1433 f.write('B %i %i %i\n' % color_to_int(cpt.color_below))
1435 if cpt.color_above:
1436 f.write('F %i %i %i\n' % color_to_int(cpt.color_above))
1438 if cpt.color_nan:
1439 f.write('N %i %i %i\n' % color_to_int(cpt.color_nan))
1442def cpt_merge_wet_dry(wet, dry):
1443 levels = []
1444 for level in wet.levels:
1445 if level.vmin < 0.:
1446 if level.vmax > 0.:
1447 level.vmax = 0.
1449 levels.append(level)
1451 for level in dry.levels:
1452 if level.vmax > 0.:
1453 if level.vmin < 0.:
1454 level.vmin = 0.
1456 levels.append(level)
1458 combi = CPT(
1459 color_below=wet.color_below,
1460 color_above=dry.color_above,
1461 color_nan=dry.color_nan,
1462 levels=levels)
1464 return combi
1467if __name__ == '__main__':
1468 from pyrocko import util
1469 util.setup_logging('pyrocko.automap', 'info')
1471 import sys
1472 if len(sys.argv) == 2:
1474 n = int(sys.argv[1])
1476 for i in range(n):
1477 m = Map(
1478 lat=rand(-60., 60.),
1479 lon=rand(-180., 180.),
1480 radius=math.exp(rand(math.log(500*km), math.log(3000*km))),
1481 width=30., height=30.,
1482 show_grid=True,
1483 show_topo=True,
1484 color_dry=(238, 236, 230),
1485 topo_cpt_wet='light_sea_uniform',
1486 topo_cpt_dry='light_land_uniform',
1487 illuminate=True,
1488 illuminate_factor_ocean=0.15,
1489 show_rivers=False,
1490 show_plates=True)
1492 m.draw_cities()
1493 print(m)
1494 m.save('map_%02i.pdf' % i)