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
25points_in_region = od.points_in_region
27logger = logging.getLogger('pyrocko.plot.automap')
29earthradius = 6371000.0
30r2d = 180./math.pi
31d2r = 1./r2d
32km = 1000.
33d2m = d2r*earthradius
34m2d = 1./d2m
35cm = gmtpy.cm
38def darken(c, f=0.7):
39 return (c[0]*f, c[1]*f, c[2]*f)
42def corners(lon, lat, w, h):
43 ll_lat, ll_lon = od.ne_to_latlon(lat, lon, -0.5*h, -0.5*w)
44 ur_lat, ur_lon = od.ne_to_latlon(lat, lon, 0.5*h, 0.5*w)
45 return ll_lon, ll_lat, ur_lon, ur_lat
48def extent(lon, lat, w, h, n):
49 x = num.linspace(-0.5*w, 0.5*w, n)
50 y = num.linspace(-0.5*h, 0.5*h, n)
51 slats, slons = od.ne_to_latlon(lat, lon, y[0], x)
52 nlats, nlons = od.ne_to_latlon(lat, lon, y[-1], x)
53 south = slats.min()
54 north = nlats.max()
56 wlats, wlons = od.ne_to_latlon(lat, lon, y, x[0])
57 elats, elons = od.ne_to_latlon(lat, lon, y, x[-1])
58 elons = num.where(elons < wlons, elons + 360., elons)
60 if elons.max() - elons.min() > 180 or wlons.max() - wlons.min() > 180.:
61 west = -180.
62 east = 180.
63 else:
64 west = wlons.min()
65 east = elons.max()
67 return topo.positive_region((west, east, south, north))
70class NoTopo(Exception):
71 pass
74class OutOfBounds(Exception):
75 pass
78class FloatTile(Object):
79 xmin = Float.T()
80 ymin = Float.T()
81 dx = Float.T()
82 dy = Float.T()
83 data = Array.T(shape=(None, None), dtype=float, serialize_as='table')
85 def __init__(self, xmin, ymin, dx, dy, data):
86 Object.__init__(self, init_props=False)
87 self.xmin = float(xmin)
88 self.ymin = float(ymin)
89 self.dx = float(dx)
90 self.dy = float(dy)
91 self.data = data
92 self._set_maxes()
94 def _set_maxes(self):
95 self.ny, self.nx = self.data.shape
96 self.xmax = self.xmin + (self.nx-1) * self.dx
97 self.ymax = self.ymin + (self.ny-1) * self.dy
99 def x(self):
100 return self.xmin + num.arange(self.nx) * self.dx
102 def y(self):
103 return self.ymin + num.arange(self.ny) * self.dy
105 def get(self, x, y):
106 ix = int(round((x - self.xmin) / self.dx))
107 iy = int(round((y - self.ymin) / self.dy))
108 if 0 <= ix < self.nx and 0 <= iy < self.ny:
109 return self.data[iy, ix]
110 else:
111 raise OutOfBounds()
114class City(Object):
115 def __init__(self, name, lat, lon, population=None, asciiname=None):
116 name = str(name)
117 lat = float(lat)
118 lon = float(lon)
119 if asciiname is None:
120 asciiname = name.encode('ascii', errors='replace')
122 if population is None:
123 population = 0
124 else:
125 population = int(population)
127 Object.__init__(self, name=name, lat=lat, lon=lon,
128 population=population, asciiname=asciiname)
130 name = Unicode.T()
131 lat = Float.T()
132 lon = Float.T()
133 population = Int.T()
134 asciiname = String.T()
137class Map(Object):
138 lat = Float.T(optional=True)
139 lon = Float.T(optional=True)
140 radius = Float.T(optional=True)
141 width = Float.T(default=20.)
142 height = Float.T(default=14.)
143 margins = List.T(Float.T())
144 illuminate = Bool.T(default=True)
145 skip_feature_factor = Float.T(default=0.02)
146 show_grid = Bool.T(default=False)
147 show_topo = Bool.T(default=True)
148 show_scale = Bool.T(default=False)
149 show_topo_scale = Bool.T(default=False)
150 show_center_mark = Bool.T(default=False)
151 show_rivers = Bool.T(default=True)
152 show_plates = Bool.T(default=False)
153 show_plate_velocities = Bool.T(default=False)
154 show_plate_names = Bool.T(default=False)
155 show_boundaries = Bool.T(default=False)
156 illuminate_factor_land = Float.T(default=0.5)
157 illuminate_factor_ocean = Float.T(default=0.25)
158 color_wet = Tuple.T(3, Int.T(), default=(216, 242, 254))
159 color_dry = Tuple.T(3, Int.T(), default=(172, 208, 165))
160 color_boundaries = Tuple.T(3, Int.T(), default=(1, 1, 1))
161 topo_resolution_min = Float.T(
162 default=40.,
163 help='minimum resolution of topography [dpi]')
164 topo_resolution_max = Float.T(
165 default=200.,
166 help='maximum resolution of topography [dpi]')
167 replace_topo_color_only = FloatTile.T(
168 optional=True,
169 help='replace topo color while keeping topographic shading')
170 topo_cpt_wet = String.T(default='light_sea')
171 topo_cpt_dry = String.T(default='light_land')
172 axes_layout = String.T(optional=True)
173 custom_cities = List.T(City.T())
174 gmt_config = Dict.T(String.T(), String.T())
175 comment = String.T(optional=True)
176 approx_ticks = Int.T(default=4)
178 def __init__(self, gmtversion='newest', **kwargs):
179 Object.__init__(self, **kwargs)
180 self._gmt = None
181 self._scaler = None
182 self._widget = None
183 self._corners = None
184 self._wesn = None
185 self._minarea = None
186 self._coastline_resolution = None
187 self._rivers = None
188 self._dems = None
189 self._have_topo_land = None
190 self._have_topo_ocean = None
191 self._jxyr = None
192 self._prep_topo_have = None
193 self._labels = []
194 self._area_labels = []
195 self._gmtversion = gmtversion
197 def save(self, outpath, resolution=75., oversample=2., size=None,
198 width=None, height=None, psconvert=False, crop_eps_mode=False):
200 '''
201 Save the image.
203 Save the image to ``outpath``. The format is determined by the filename
204 extension. Formats are handled as follows: ``'.eps'`` and ``'.ps'``
205 produce EPS and PS, respectively, directly with GMT. If the file name
206 ends with ``'.pdf'``, GMT output is fed through ``gmtpy-epstopdf`` to
207 create a PDF file. For any other filename extension, output is first
208 converted to PDF with ``gmtpy-epstopdf``, then with ``pdftocairo`` to
209 PNG with a resolution oversampled by the factor ``oversample`` and
210 finally the PNG is downsampled and converted to the target format with
211 ``convert``. The resolution of rasterized target image can be
212 controlled either by ``resolution`` in DPI or by specifying ``width``
213 or ``height`` or ``size``, where the latter fits the image into a
214 square with given side length. To save transparency use
215 ``psconvert=True``. To crop the output image with a rectangle to the
216 nearest non-white element set ``crop_eps_mode=True``.
217 '''
219 gmt = self.gmt
220 self.draw_labels()
221 self.draw_axes()
222 if self.show_topo and self.show_topo_scale:
223 self._draw_topo_scale()
225 gmt.save(outpath, resolution=resolution, oversample=oversample,
226 crop_eps_mode=crop_eps_mode,
227 size=size, width=width, height=height, psconvert=psconvert)
229 @property
230 def scaler(self):
231 if self._scaler is None:
232 self._setup_geometry()
234 return self._scaler
236 @property
237 def wesn(self):
238 if self._wesn is None:
239 self._setup_geometry()
241 return self._wesn
243 @property
244 def widget(self):
245 if self._widget is None:
246 self._setup()
248 return self._widget
250 @property
251 def layout(self):
252 if self._layout is None:
253 self._setup()
255 return self._layout
257 @property
258 def jxyr(self):
259 if self._jxyr is None:
260 self._setup()
262 return self._jxyr
264 @property
265 def pxyr(self):
266 if self._pxyr is None:
267 self._setup()
269 return self._pxyr
271 @property
272 def gmt(self):
273 if self._gmt is None:
274 self._setup()
276 if self._have_topo_ocean is None:
277 self._draw_background()
279 return self._gmt
281 def _setup(self):
282 if not self._widget:
283 self._setup_geometry()
285 self._setup_lod()
286 self._setup_gmt()
288 def _setup_geometry(self):
289 wpage, hpage = self.width, self.height
290 ml, mr, mt, mb = self._expand_margins()
291 wpage -= ml + mr
292 hpage -= mt + mb
294 wreg = self.radius * 2.0
295 hreg = self.radius * 2.0
296 if wpage >= hpage:
297 wreg *= wpage/hpage
298 else:
299 hreg *= hpage/wpage
301 self._wreg = wreg
302 self._hreg = hreg
304 self._corners = corners(self.lon, self.lat, wreg, hreg)
305 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10)
307 x, y, z = ((west, east), (south, north), (-6000., 4500.))
309 xax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
310 yax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
311 zax = gmtpy.Ax(mode='min-max', inc=1000., label='Height',
312 scaled_unit='km', scaled_unit_factor=0.001)
314 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax))
316 par = scaler.get_params()
318 west = par['xmin']
319 east = par['xmax']
320 south = par['ymin']
321 north = par['ymax']
323 self._wesn = west, east, south, north
324 self._scaler = scaler
326 def _setup_lod(self):
327 w, e, s, n = self._wesn
328 if self.radius > 1500.*km:
329 coastline_resolution = 'i'
330 rivers = False
331 else:
332 coastline_resolution = 'f'
333 rivers = True
335 self._minarea = (self.skip_feature_factor * self.radius/km)**2
337 self._coastline_resolution = coastline_resolution
338 self._rivers = rivers
340 self._prep_topo_have = {}
341 self._dems = {}
343 cm2inch = gmtpy.cm/gmtpy.inch
345 dmin = 2.0 * self.radius * m2d / (self.topo_resolution_max *
346 (self.height * cm2inch))
347 dmax = 2.0 * self.radius * m2d / (self.topo_resolution_min *
348 (self.height * cm2inch))
350 for k in ['ocean', 'land']:
351 self._dems[k] = topo.select_dem_names(k, dmin, dmax, self._wesn)
352 if self._dems[k]:
353 logger.debug('using topography dataset %s for %s'
354 % (','.join(self._dems[k]), k))
356 def _expand_margins(self):
357 if len(self.margins) == 0 or len(self.margins) > 4:
358 ml = mr = mt = mb = 2.0
359 elif len(self.margins) == 1:
360 ml = mr = mt = mb = self.margins[0]
361 elif len(self.margins) == 2:
362 ml = mr = self.margins[0]
363 mt = mb = self.margins[1]
364 elif len(self.margins) == 4:
365 ml, mr, mt, mb = self.margins
367 return ml, mr, mt, mb
369 def _setup_gmt(self):
370 w, h = self.width, self.height
371 scaler = self._scaler
373 if gmtpy.is_gmt5(self._gmtversion):
374 gmtconf = dict(
375 MAP_TICK_PEN_PRIMARY='1.25p',
376 MAP_TICK_PEN_SECONDARY='1.25p',
377 MAP_TICK_LENGTH_PRIMARY='0.2c',
378 MAP_TICK_LENGTH_SECONDARY='0.6c',
379 FONT_ANNOT_PRIMARY='12p,1,black',
380 FONT_LABEL='12p,1,black',
381 PS_CHAR_ENCODING='ISOLatin1+',
382 MAP_FRAME_TYPE='fancy',
383 FORMAT_GEO_MAP='D',
384 PS_MEDIA='Custom_%ix%i' % (
385 w*gmtpy.cm,
386 h*gmtpy.cm),
387 PS_PAGE_ORIENTATION='portrait',
388 MAP_GRID_PEN_PRIMARY='thinnest,0/50/0',
389 MAP_ANNOT_OBLIQUE='6')
390 else:
391 gmtconf = dict(
392 TICK_PEN='1.25p',
393 TICK_LENGTH='0.2c',
394 ANNOT_FONT_PRIMARY='1',
395 ANNOT_FONT_SIZE_PRIMARY='12p',
396 LABEL_FONT='1',
397 LABEL_FONT_SIZE='12p',
398 CHAR_ENCODING='ISOLatin1+',
399 BASEMAP_TYPE='fancy',
400 PLOT_DEGREE_FORMAT='D',
401 PAPER_MEDIA='Custom_%ix%i' % (
402 w*gmtpy.cm,
403 h*gmtpy.cm),
404 GRID_PEN_PRIMARY='thinnest/0/50/0',
405 DOTS_PR_INCH='1200',
406 OBLIQUE_ANNOTATION='6')
408 gmtconf.update(
409 (k.upper(), v) for (k, v) in self.gmt_config.items())
411 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion)
413 layout = gmt.default_layout()
415 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()])
417 widget = layout.get_widget()
418 widget['P'] = widget['J']
419 widget['J'] = ('-JA%g/%g' % (self.lon, self.lat)) + '/%(width)gp'
420 scaler['R'] = '-R%g/%g/%g/%gr' % self._corners
422 # aspect = gmtpy.aspect_for_projection(
423 # gmt.installation['version'], *(widget.J() + scaler.R()))
425 aspect = self._map_aspect(jr=widget.J() + scaler.R())
426 widget.set_aspect(aspect)
428 self._gmt = gmt
429 self._layout = layout
430 self._widget = widget
431 self._jxyr = self._widget.JXY() + self._scaler.R()
432 self._pxyr = self._widget.PXY() + [
433 '-R%g/%g/%g/%g' % (0, widget.width(), 0, widget.height())]
434 self._have_drawn_axes = False
435 self._have_drawn_labels = False
437 def _draw_background(self):
438 self._have_topo_land = False
439 self._have_topo_ocean = False
440 if self.show_topo:
441 self._have_topo = self._draw_topo()
443 self._draw_basefeatures()
445 def _get_topo_tile(self, k):
446 t = None
447 demname = None
448 for dem in self._dems[k]:
449 t = topo.get(dem, self._wesn)
450 demname = dem
451 if t is not None:
452 break
454 if not t:
455 raise NoTopo()
457 return t, demname
459 def _prep_topo(self, k):
460 gmt = self._gmt
461 t, demname = self._get_topo_tile(k)
463 if demname not in self._prep_topo_have:
465 grdfile = gmt.tempfilename()
467 is_flat = num.all(t.data[0] == t.data)
469 gmtpy.savegrd(
470 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat')
472 if self.illuminate and not is_flat:
473 if k == 'ocean':
474 factor = self.illuminate_factor_ocean
475 else:
476 factor = self.illuminate_factor_land
478 ilumfn = gmt.tempfilename()
479 gmt.grdgradient(
480 grdfile,
481 N='e%g' % factor,
482 A=-45,
483 G=ilumfn,
484 out_discard=True)
486 ilumargs = ['-I%s' % ilumfn]
487 else:
488 ilumargs = []
490 if self.replace_topo_color_only:
491 t2 = self.replace_topo_color_only
492 grdfile2 = gmt.tempfilename()
494 gmtpy.savegrd(
495 t2.x(), t2.y(), t2.data, filename=grdfile2,
496 naming='lonlat')
498 if gmt.is_gmt5():
499 gmt.grdsample(
500 grdfile2,
501 G=grdfile,
502 n='l',
503 I='%g/%g' % (t.dx, t.dy), # noqa
504 R=grdfile,
505 out_discard=True)
506 else:
507 gmt.grdsample(
508 grdfile2,
509 G=grdfile,
510 Q='l',
511 I='%g/%g' % (t.dx, t.dy), # noqa
512 R=grdfile,
513 out_discard=True)
515 gmt.grdmath(
516 grdfile, '0.0', 'AND', '=', grdfile2,
517 out_discard=True)
519 grdfile = grdfile2
521 self._prep_topo_have[demname] = grdfile, ilumargs
523 return self._prep_topo_have[demname]
525 def _draw_topo(self):
526 widget = self._widget
527 scaler = self._scaler
528 gmt = self._gmt
529 cres = self._coastline_resolution
530 minarea = self._minarea
532 JXY = widget.JXY()
533 R = scaler.R()
535 try:
536 grdfile, ilumargs = self._prep_topo('ocean')
537 gmt.pscoast(D=cres, S='c', A=minarea, *(JXY+R))
538 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_wet),
539 *(ilumargs+JXY+R))
540 gmt.pscoast(Q=True, *(JXY+R))
541 self._have_topo_ocean = True
542 except NoTopo:
543 self._have_topo_ocean = False
545 try:
546 grdfile, ilumargs = self._prep_topo('land')
547 gmt.pscoast(D=cres, G='c', A=minarea, *(JXY+R))
548 gmt.grdimage(grdfile, C=topo.cpt(self.topo_cpt_dry),
549 *(ilumargs+JXY+R))
550 gmt.pscoast(Q=True, *(JXY+R))
551 self._have_topo_land = True
552 except NoTopo:
553 self._have_topo_land = False
555 def _draw_topo_scale(self, label='Elevation [km]'):
556 dry = read_cpt(topo.cpt(self.topo_cpt_dry))
557 wet = read_cpt(topo.cpt(self.topo_cpt_wet))
558 combi = cpt_merge_wet_dry(wet, dry)
559 for level in combi.levels:
560 level.vmin /= km
561 level.vmax /= km
563 topo_cpt = self.gmt.tempfilename() + '.cpt'
564 write_cpt(combi, topo_cpt)
566 (w, h), (xo, yo) = self.widget.get_size()
567 self.gmt.psscale(
568 D='%gp/%gp/%gp/%gph' % (xo + 0.5*w, yo - 2.0*gmtpy.cm, w,
569 0.5*gmtpy.cm),
570 C=topo_cpt,
571 B='1:%s:' % label)
573 def _draw_basefeatures(self):
574 gmt = self._gmt
575 cres = self._coastline_resolution
576 rivers = self._rivers
577 minarea = self._minarea
579 color_wet = self.color_wet
580 color_dry = self.color_dry
582 if self.show_rivers and rivers:
583 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)]
584 else:
585 rivers = []
587 fill = {}
588 if not self._have_topo_land:
589 fill['G'] = color_dry
591 if not self._have_topo_ocean:
592 fill['S'] = color_wet
594 if self.show_boundaries:
595 fill['N'] = '1/1p,%s,%s' % (
596 gmtpy.color(self.color_boundaries), 'solid')
598 gmt.pscoast(
599 D=cres,
600 W='thinnest,%s' % gmtpy.color(darken(gmtpy.color_tup(color_dry))),
601 A=minarea,
602 *(rivers+self._jxyr), **fill)
604 if self.show_plates:
605 self.draw_plates()
607 def _draw_axes(self):
608 gmt = self._gmt
609 scaler = self._scaler
610 widget = self._widget
612 if self.axes_layout is None:
613 if self.lat > 0.0:
614 axes_layout = 'WSen'
615 else:
616 axes_layout = 'WseN'
617 else:
618 axes_layout = self.axes_layout
620 scale_km = gmtpy.nice_value(self.radius/5.) / 1000.
622 if self.show_center_mark:
623 gmt.psxy(
624 in_rows=[[self.lon, self.lat]],
625 S='c20p', W='2p,black',
626 *self._jxyr)
628 if self.show_grid:
629 btmpl = ('%(xinc)gg%(xinc)g:%(xlabel)s:/'
630 '%(yinc)gg%(yinc)g:%(ylabel)s:')
631 else:
632 btmpl = '%(xinc)g:%(xlabel)s:/%(yinc)g:%(ylabel)s:'
634 if self.show_scale:
635 scale = 'x%gp/%gp/%g/%g/%gk' % (
636 6./7*widget.width(),
637 widget.height()/7.,
638 self.lon,
639 self.lat,
640 scale_km)
641 else:
642 scale = False
644 gmt.psbasemap(
645 B=(btmpl % scaler.get_params())+axes_layout,
646 L=scale,
647 *self._jxyr)
649 if self.comment:
650 font_size = self.gmt.label_font_size()
652 _, east, south, _ = self._wesn
653 if gmt.is_gmt5():
654 row = [
655 1, 0,
656 '%gp,%s,%s' % (font_size, 0, 'black'), 'BR',
657 self.comment]
659 farg = ['-F+f+j']
660 else:
661 row = [1, 0, font_size, 0, 0, 'BR', self.comment]
662 farg = []
664 gmt.pstext(
665 in_rows=[row],
666 N=True,
667 R=(0, 1, 0, 1),
668 D='%gp/%gp' % (-font_size*0.2, font_size*0.3),
669 *(widget.PXY() + farg))
671 def draw_axes(self):
672 if not self._have_drawn_axes:
673 self._draw_axes()
674 self._have_drawn_axes = True
676 def _have_coastlines(self):
677 gmt = self._gmt
678 cres = self._coastline_resolution
679 minarea = self._minarea
681 checkfile = gmt.tempfilename()
683 gmt.pscoast(
684 M=True,
685 D=cres,
686 W='thinnest,black',
687 A=minarea,
688 out_filename=checkfile,
689 *self._jxyr)
691 points = []
692 with open(checkfile, 'r') as f:
693 for line in f:
694 ls = line.strip()
695 if ls.startswith('#') or ls.startswith('>') or ls == '':
696 continue
697 plon, plat = [float(x) for x in ls.split()]
698 points.append((plat, plon))
700 points = num.array(points, dtype=float)
701 return num.any(points_in_region(points, self._wesn))
703 def have_coastlines(self):
704 self.gmt
705 return self._have_coastlines()
707 def project(self, lats, lons, jr=None):
708 onepoint = False
709 if isinstance(lats, float) and isinstance(lons, float):
710 lats = [lats]
711 lons = [lons]
712 onepoint = True
714 if jr is not None:
715 j, r = jr
716 gmt = gmtpy.GMT(version=self._gmtversion)
717 else:
718 j, _, _, r = self.jxyr
719 gmt = self.gmt
721 f = BytesIO()
722 gmt.mapproject(j, r, in_columns=(lons, lats), out_stream=f, D='p')
723 f.seek(0)
724 data = num.loadtxt(f, ndmin=2)
725 xs, ys = data.T
726 if onepoint:
727 xs = xs[0]
728 ys = ys[0]
729 return xs, ys
731 def _map_box(self, jr=None):
732 ll_lon, ll_lat, ur_lon, ur_lat = self._corners
734 xs_corner, ys_corner = self.project(
735 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr)
737 w = xs_corner[1] - xs_corner[0]
738 h = ys_corner[1] - ys_corner[0]
740 return w, h
742 def _map_aspect(self, jr=None):
743 w, h = self._map_box(jr=jr)
744 return h/w
746 def _draw_labels(self):
747 points_taken = []
748 regions_taken = []
750 def no_points_in_rect(xs, ys, xmin, ymin, xmax, ymax):
751 xx = not num.any(la(la(xmin < xs, xs < xmax),
752 la(ymin < ys, ys < ymax)))
753 return xx
755 def roverlaps(a, b):
756 return (a[0] < b[2] and b[0] < a[2] and
757 a[1] < b[3] and b[1] < a[3])
759 w, h = self._map_box()
761 label_font_size = self.gmt.label_font_size()
763 if self._labels:
765 n = len(self._labels)
767 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \
768 angles, styles = list(zip(*self._labels))
770 font_sizes = [
771 (font_size or label_font_size) for font_size in font_sizes]
773 sx = num.array(sx, dtype=float)
774 sy = num.array(sy, dtype=float)
776 xs, ys = self.project(lats, lons)
778 points_taken.append((xs, ys))
780 dxs = num.zeros(n)
781 dys = num.zeros(n)
783 for i in range(n):
784 dx, dy = gmtpy.text_box(
785 texts[i],
786 font=fonts[i],
787 font_size=font_sizes[i],
788 **styles[i])
790 dxs[i] = dx
791 dys[i] = dy
793 la = num.logical_and
794 anchors_ok = (
795 la(xs + sx + dxs < w, ys + sy + dys < h),
796 la(xs - sx - dxs > 0., ys - sy - dys > 0.),
797 la(xs + sx + dxs < w, ys - sy - dys > 0.),
798 la(xs - sx - dxs > 0., ys + sy + dys < h),
799 )
801 arects = [
802 (xs, ys, xs + sx + dxs, ys + sy + dys),
803 (xs - sx - dxs, ys - sy - dys, xs, ys),
804 (xs, ys - sy - dys, xs + sx + dxs, ys),
805 (xs - sx - dxs, ys, xs, ys + sy + dys)]
807 for i in range(n):
808 for ianch in range(4):
809 anchors_ok[ianch][i] &= no_points_in_rect(
810 xs, ys, *[xxx[i] for xxx in arects[ianch]])
812 anchor_choices = []
813 anchor_take = []
814 for i in range(n):
815 choices = [ianch for ianch in range(4)
816 if anchors_ok[ianch][i]]
817 anchor_choices.append(choices)
818 if choices:
819 anchor_take.append(choices[0])
820 else:
821 anchor_take.append(None)
823 def cost(anchor_take):
824 noverlaps = 0
825 for i in range(n):
826 for j in range(n):
827 if i != j:
828 i_take = anchor_take[i]
829 j_take = anchor_take[j]
830 if i_take is None or j_take is None:
831 continue
832 r_i = [xxx[i] for xxx in arects[i_take]]
833 r_j = [xxx[j] for xxx in arects[j_take]]
834 if roverlaps(r_i, r_j):
835 noverlaps += 1
837 return noverlaps
839 cur_cost = cost(anchor_take)
840 imax = 30
841 while cur_cost != 0 and imax > 0:
842 for i in range(n):
843 for t in anchor_choices[i]:
844 anchor_take_new = list(anchor_take)
845 anchor_take_new[i] = t
846 new_cost = cost(anchor_take_new)
847 if new_cost < cur_cost:
848 anchor_take = anchor_take_new
849 cur_cost = new_cost
851 imax -= 1
853 while cur_cost != 0:
854 for i in range(n):
855 anchor_take_new = list(anchor_take)
856 anchor_take_new[i] = None
857 new_cost = cost(anchor_take_new)
858 if new_cost < cur_cost:
859 anchor_take = anchor_take_new
860 cur_cost = new_cost
861 break
863 anchor_strs = ['BL', 'TR', 'TL', 'BR']
865 for i in range(n):
866 ianchor = anchor_take[i]
867 color = colors[i]
868 if color is None:
869 color = 'black'
871 if ianchor is not None:
872 regions_taken.append([xxx[i] for xxx in arects[ianchor]])
874 anchor = anchor_strs[ianchor]
876 yoff = [-sy[i], sy[i]][anchor[0] == 'B']
877 xoff = [-sx[i], sx[i]][anchor[1] == 'L']
878 if self.gmt.is_gmt5():
879 row = (
880 lons[i], lats[i],
881 '%i,%s,%s' % (font_sizes[i], fonts[i], color),
882 anchor,
883 texts[i])
885 farg = ['-F+f+j+a%g' % angles[i]]
886 else:
887 row = (
888 lons[i], lats[i],
889 font_sizes[i], angles[i], fonts[i], anchor,
890 texts[i])
891 farg = ['-G%s' % color]
893 self.gmt.pstext(
894 in_rows=[row],
895 D='%gp/%gp' % (xoff, yoff),
896 *(self.jxyr + farg),
897 **styles[i])
899 if self._area_labels:
901 for lons, lats, text, color, font, font_size, style in \
902 self._area_labels:
904 if font_size is None:
905 font_size = label_font_size
907 if color is None:
908 color = 'black'
910 if self.gmt.is_gmt5():
911 farg = ['-F+f+j']
912 else:
913 farg = ['-G%s' % color]
915 xs, ys = self.project(lats, lons)
916 dx, dy = gmtpy.text_box(
917 text, font=font, font_size=font_size, **style)
919 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy]
921 locs_ok = num.ones(xs.size, dtype=bool)
923 for iloc in range(xs.size):
924 rcandi = [xxx[iloc] for xxx in rects]
926 locs_ok[iloc] = True
927 locs_ok[iloc] &= (
928 0 < rcandi[0] and rcandi[2] < w
929 and 0 < rcandi[1] and rcandi[3] < h)
931 overlap = False
932 for r in regions_taken:
933 if roverlaps(r, rcandi):
934 overlap = True
935 break
937 locs_ok[iloc] &= not overlap
939 for xs_taken, ys_taken in points_taken:
940 locs_ok[iloc] &= no_points_in_rect(
941 xs_taken, ys_taken, *rcandi)
943 if not locs_ok[iloc]:
944 break
946 rows = []
947 for iloc, (lon, lat) in enumerate(zip(lons, lats)):
948 if not locs_ok[iloc]:
949 continue
951 if self.gmt.is_gmt5():
952 row = (
953 lon, lat,
954 '%i,%s,%s' % (font_size, font, color),
955 'MC',
956 text)
958 else:
959 row = (
960 lon, lat,
961 font_size, 0, font, 'MC',
962 text)
964 rows.append(row)
966 regions_taken.append([xxx[iloc] for xxx in rects])
967 break
969 self.gmt.pstext(
970 in_rows=rows,
971 *(self.jxyr + farg),
972 **style)
974 def draw_labels(self):
975 self.gmt
976 if not self._have_drawn_labels:
977 self._draw_labels()
978 self._have_drawn_labels = True
980 def add_label(
981 self, lat, lon, text,
982 offset_x=5., offset_y=5.,
983 color=None,
984 font='1',
985 font_size=None,
986 angle=0,
987 style={}):
989 if 'G' in style:
990 style = style.copy()
991 color = style.pop('G')
993 self._labels.append(
994 (lon, lat, text, offset_x, offset_y, color, font, font_size,
995 angle, style))
997 def add_area_label(
998 self, lat, lon, text,
999 color=None,
1000 font='3',
1001 font_size=None,
1002 style={}):
1004 self._area_labels.append(
1005 (lon, lat, text, color, font, font_size, style))
1007 def cities_in_region(self):
1008 from pyrocko.dataset import geonames
1009 cities = geonames.get_cities_region(region=self.wesn, minpop=0)
1010 cities.extend(self.custom_cities)
1011 cities.sort(key=lambda x: x.population)
1012 return cities
1014 def draw_cities(self,
1015 exact=None,
1016 include=[],
1017 exclude=[],
1018 nmax_soft=10,
1019 psxy_style=dict(S='s5p', G='black')):
1021 cities = self.cities_in_region()
1023 if exact is not None:
1024 cities = [c for c in cities if c.name in exact]
1025 minpop = None
1026 else:
1027 cities = [c for c in cities if c.name not in exclude]
1028 minpop = 10**3
1029 for minpop_new in [1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6, 3e6, 1e7]:
1030 cities_new = [
1031 c for c in cities
1032 if c.population > minpop_new or c.name in include]
1034 if len(cities_new) == 0 or (
1035 len(cities_new) < 3 and len(cities) < nmax_soft*2):
1036 break
1038 cities = cities_new
1039 minpop = minpop_new
1040 if len(cities) <= nmax_soft:
1041 break
1043 if cities:
1044 lats = [c.lat for c in cities]
1045 lons = [c.lon for c in cities]
1047 self.gmt.psxy(
1048 in_columns=(lons, lats),
1049 *self.jxyr, **psxy_style)
1051 for c in cities:
1052 try:
1053 text = c.name.encode('iso-8859-1').decode('iso-8859-1')
1054 except UnicodeEncodeError:
1055 text = c.asciiname
1057 self.add_label(c.lat, c.lon, text)
1059 self._cities_minpop = minpop
1061 def add_stations(self, stations, psxy_style=dict()):
1063 default_psxy_style = {
1064 'S': 't8p',
1065 'G': 'black'
1066 }
1067 default_psxy_style.update(psxy_style)
1069 lats, lons = zip(*[s.effective_latlon for s in stations])
1071 self.gmt.psxy(
1072 in_columns=(lons, lats),
1073 *self.jxyr, **default_psxy_style)
1075 for station in stations:
1076 self.add_label(
1077 station.effective_lat,
1078 station.effective_lon,
1079 '.'.join(x for x in (station.network, station.station) if x))
1081 def add_kite_scene(self, scene):
1082 tile = FloatTile(
1083 scene.frame.llLon,
1084 scene.frame.llLat,
1085 scene.frame.dLon,
1086 scene.frame.dLat,
1087 scene.displacement)
1089 return tile
1091 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None,
1092 labels=True, vertical=False, fontsize=10):
1094 stations = campaign.stations
1096 if offset_scale is None:
1097 offset_scale = num.zeros(campaign.nstations)
1098 for ista, sta in enumerate(stations):
1099 for comp in sta.components.values():
1100 offset_scale[ista] += comp.shift
1101 offset_scale = num.sqrt(offset_scale**2).max()
1103 size = math.sqrt(self.height**2 + self.width**2)
1104 scale = (size/10.) / offset_scale
1105 logger.debug('GNSS: Using offset scale %f, map scale %f',
1106 offset_scale, scale)
1108 lats, lons = zip(*[s.effective_latlon for s in stations])
1110 if self.gmt.is_gmt6():
1111 sign_factor = 1.
1112 arrow_head_placement = 'e'
1113 else:
1114 sign_factor = -1.
1115 arrow_head_placement = 'b'
1117 if vertical:
1118 rows = [[lons[ista], lats[ista],
1119 0., sign_factor * s.up.shift,
1120 (s.east.sigma + s.north.sigma) if s.east.sigma else 0.,
1121 s.up.sigma, 0.,
1122 s.code if labels else None]
1123 for ista, s in enumerate(stations)
1124 if s.up is not None]
1126 else:
1127 rows = [[lons[ista], lats[ista],
1128 sign_factor * s.east.shift, sign_factor * s.north.shift,
1129 s.east.sigma, s.north.sigma, s.correlation_ne,
1130 s.code if labels else None]
1131 for ista, s in enumerate(stations)
1132 if s.east is not None or s.north is not None]
1134 default_psxy_style = {
1135 'h': 0,
1136 'W': '2p,black',
1137 'A': '+p2p,black+{}+a40'.format(arrow_head_placement),
1138 'G': 'black',
1139 'L': True,
1140 'S': 'e%dc/0.95/%d' % (scale, fontsize),
1141 }
1143 if not labels:
1144 for row in rows:
1145 row.pop(-1)
1147 if psxy_style is not None:
1148 default_psxy_style.update(psxy_style)
1150 self.gmt.psvelo(
1151 in_rows=rows,
1152 *self.jxyr,
1153 **default_psxy_style)
1155 def draw_plates(self):
1156 from pyrocko.dataset import tectonics
1158 neast = 20
1159 nnorth = max(1, int(round(num.round(self._hreg/self._wreg * neast))))
1160 norths = num.linspace(-self._hreg*0.5, self._hreg*0.5, nnorth)
1161 easts = num.linspace(-self._wreg*0.5, self._wreg*0.5, neast)
1162 norths2 = num.repeat(norths, neast)
1163 easts2 = num.tile(easts, nnorth)
1164 lats, lons = od.ne_to_latlon(
1165 self.lat, self.lon, norths2, easts2)
1167 bird = tectonics.PeterBird2003()
1168 plates = bird.get_plates()
1170 color_plates = gmtpy.color('aluminium5')
1171 color_velocities = gmtpy.color('skyblue1')
1172 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1')))
1174 points = num.vstack((lats, lons)).T
1175 used = []
1176 for plate in plates:
1177 mask = plate.contains_points(points)
1178 if num.any(mask):
1179 used.append((plate, mask))
1181 if len(used) > 1:
1183 candi_fixed = {}
1185 label_data = []
1186 for plate, mask in used:
1188 mean_north = num.mean(norths2[mask])
1189 mean_east = num.mean(easts2[mask])
1190 iorder = num.argsort(num.sqrt(
1191 (norths2[mask] - mean_north)**2 +
1192 (easts2[mask] - mean_east)**2))
1194 lat_candis = lats[mask][iorder]
1195 lon_candis = lons[mask][iorder]
1197 candi_fixed[plate.name] = lat_candis.size
1199 label_data.append((
1200 lat_candis, lon_candis, plate, color_plates))
1202 boundaries = bird.get_boundaries()
1204 size = 1.
1206 psxy_kwargs = []
1208 for boundary in boundaries:
1209 if num.any(points_in_region(boundary.points, self._wesn)):
1210 for typ, part in boundary.split_types(
1211 [['SUB'],
1212 ['OSR', 'CRB'],
1213 ['OTF', 'CTF', 'OCB', 'CCB']]):
1215 lats, lons = part.T
1217 kwargs = {}
1219 kwargs['in_columns'] = (lons, lats)
1220 kwargs['W'] = '%gp,%s' % (size, color_plates)
1222 if typ[0] == 'SUB':
1223 if boundary.kind == '\\':
1224 kwargs['S'] = 'f%g/%gp+t+r' % (
1225 0.45*size, 3.*size)
1226 elif boundary.kind == '/':
1227 kwargs['S'] = 'f%g/%gp+t+l' % (
1228 0.45*size, 3.*size)
1230 kwargs['G'] = color_plates
1232 elif typ[0] in ['OSR', 'CRB']:
1233 kwargs_bg = {}
1234 kwargs_bg['in_columns'] = (lons, lats)
1235 kwargs_bg['W'] = '%gp,%s' % (
1236 size * 3, color_plates)
1237 psxy_kwargs.append(kwargs_bg)
1239 kwargs['W'] = '%gp,%s' % (size * 2, 'white')
1241 psxy_kwargs.append(kwargs)
1243 if boundary.kind == '\\':
1244 if boundary.plate_name2 in candi_fixed:
1245 candi_fixed[boundary.plate_name2] += \
1246 neast*nnorth
1248 elif boundary.kind == '/':
1249 if boundary.plate_name1 in candi_fixed:
1250 candi_fixed[boundary.plate_name1] += \
1251 neast*nnorth
1253 candi_fixed = [name for name in sorted(
1254 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])]
1256 candi_fixed.append(None)
1258 gsrm = tectonics.GSRM1()
1260 for name in candi_fixed:
1261 if name not in gsrm.plate_names() \
1262 and name not in gsrm.plate_alt_names():
1264 continue
1266 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \
1267 gsrm.get_velocities(name, region=self._wesn)
1269 fixed_plate_name = name
1271 if self.show_plate_velocities:
1272 self.gmt.psvelo(
1273 in_columns=(
1274 lons, lats, veast, vnorth, veast_err, vnorth_err,
1275 corr),
1276 W='0.25p,%s' % color_velocities,
1277 A='9p+e+g%s' % color_velocities,
1278 S='e0.2p/0.95/10',
1279 *self.jxyr)
1281 for _ in range(len(lons) // 50 + 1):
1282 ii = random.randint(0, len(lons)-1)
1283 v = math.sqrt(vnorth[ii]**2 + veast[ii]**2)
1284 self.add_label(
1285 lats[ii], lons[ii], '%.0f' % v,
1286 font_size=0.7*self.gmt.label_font_size(),
1287 style=dict(
1288 G=color_velocities_lab))
1290 break
1292 if self.show_plate_names:
1293 for (lat_candis, lon_candis, plate, color) in label_data:
1294 full_name = bird.full_name(plate.name)
1295 if plate.name == fixed_plate_name:
1296 full_name = '@_' + full_name + '@_'
1298 self.add_area_label(
1299 lat_candis, lon_candis,
1300 full_name,
1301 color=color,
1302 font='3')
1304 for kwargs in psxy_kwargs:
1305 self.gmt.psxy(*self.jxyr, **kwargs)
1308def rand(mi, ma):
1309 mi = float(mi)
1310 ma = float(ma)
1311 return random.random() * (ma-mi) + mi
1314def split_region(region):
1315 west, east, south, north = topo.positive_region(region)
1316 if east > 180:
1317 return [(west, 180., south, north),
1318 (-180., east-360., south, north)]
1319 else:
1320 return [region]
1323class CPTLevel(Object):
1324 vmin = Float.T()
1325 vmax = Float.T()
1326 color_min = Tuple.T(3, Float.T())
1327 color_max = Tuple.T(3, Float.T())
1330class CPT(Object):
1331 color_below = Tuple.T(3, Float.T(), optional=True)
1332 color_above = Tuple.T(3, Float.T(), optional=True)
1333 color_nan = Tuple.T(3, Float.T(), optional=True)
1334 levels = List.T(CPTLevel.T())
1336 def scale(self, vmin, vmax):
1337 vmin_old, vmax_old = self.levels[0].vmin, self.levels[-1].vmax
1338 for level in self.levels:
1339 level.vmin = (level.vmin - vmin_old) / (vmax_old - vmin_old) * \
1340 (vmax - vmin) + vmin
1341 level.vmax = (level.vmax - vmin_old) / (vmax_old - vmin_old) * \
1342 (vmax - vmin) + vmin
1344 def discretize(self, nlevels):
1345 colors = []
1346 vals = []
1347 for level in self.levels:
1348 vals.append(level.vmin)
1349 vals.append(level.vmax)
1350 colors.append(level.color_min)
1351 colors.append(level.color_max)
1353 r, g, b = num.array(colors, dtype=float).T
1354 vals = num.array(vals, dtype=float)
1356 vmin, vmax = self.levels[0].vmin, self.levels[-1].vmax
1357 x = num.linspace(vmin, vmax, nlevels+1)
1358 rd = num.interp(x, vals, r)
1359 gd = num.interp(x, vals, g)
1360 bd = num.interp(x, vals, b)
1362 levels = []
1363 for ilevel in range(nlevels):
1364 color = (
1365 float(0.5*(rd[ilevel]+rd[ilevel+1])),
1366 float(0.5*(gd[ilevel]+gd[ilevel+1])),
1367 float(0.5*(bd[ilevel]+bd[ilevel+1])))
1369 levels.append(CPTLevel(
1370 vmin=x[ilevel],
1371 vmax=x[ilevel+1],
1372 color_min=color,
1373 color_max=color))
1375 cpt = CPT(
1376 color_below=self.color_below,
1377 color_above=self.color_above,
1378 color_nan=self.color_nan,
1379 levels=levels)
1381 return cpt
1384class CPTParseError(Exception):
1385 pass
1388def read_cpt(filename):
1389 with open(filename) as f:
1390 color_below = None
1391 color_above = None
1392 color_nan = None
1393 levels = []
1394 try:
1395 for line in f:
1396 line = line.strip()
1397 toks = line.split()
1399 if line.startswith('#'):
1400 continue
1402 elif line.startswith('B'):
1403 color_below = tuple(map(float, toks[1:4]))
1405 elif line.startswith('F'):
1406 color_above = tuple(map(float, toks[1:4]))
1408 elif line.startswith('N'):
1409 color_nan = tuple(map(float, toks[1:4]))
1411 else:
1412 values = list(map(float, line.split()))
1413 vmin = values[0]
1414 color_min = tuple(values[1:4])
1415 vmax = values[4]
1416 color_max = tuple(values[5:8])
1417 levels.append(CPTLevel(
1418 vmin=vmin,
1419 vmax=vmax,
1420 color_min=color_min,
1421 color_max=color_max))
1423 except Exception:
1424 raise CPTParseError()
1426 return CPT(
1427 color_below=color_below,
1428 color_above=color_above,
1429 color_nan=color_nan,
1430 levels=levels)
1433def color_to_int(color):
1434 return tuple(max(0, min(255, int(round(x)))) for x in color)
1437def write_cpt(cpt, filename):
1438 with open(filename, 'w') as f:
1439 for level in cpt.levels:
1440 f.write(
1441 '%e %i %i %i %e %i %i %i\n' %
1442 ((level.vmin, ) + color_to_int(level.color_min) +
1443 (level.vmax, ) + color_to_int(level.color_max)))
1445 if cpt.color_below:
1446 f.write('B %i %i %i\n' % color_to_int(cpt.color_below))
1448 if cpt.color_above:
1449 f.write('F %i %i %i\n' % color_to_int(cpt.color_above))
1451 if cpt.color_nan:
1452 f.write('N %i %i %i\n' % color_to_int(cpt.color_nan))
1455def cpt_merge_wet_dry(wet, dry):
1456 levels = []
1457 for level in wet.levels:
1458 if level.vmin < 0.:
1459 if level.vmax > 0.:
1460 level.vmax = 0.
1462 levels.append(level)
1464 for level in dry.levels:
1465 if level.vmax > 0.:
1466 if level.vmin < 0.:
1467 level.vmin = 0.
1469 levels.append(level)
1471 combi = CPT(
1472 color_below=wet.color_below,
1473 color_above=dry.color_above,
1474 color_nan=dry.color_nan,
1475 levels=levels)
1477 return combi
1480if __name__ == '__main__':
1481 from pyrocko import util
1482 util.setup_logging('pyrocko.automap', 'info')
1484 import sys
1485 if len(sys.argv) == 2:
1487 n = int(sys.argv[1])
1489 for i in range(n):
1490 m = Map(
1491 lat=rand(-60., 60.),
1492 lon=rand(-180., 180.),
1493 radius=math.exp(rand(math.log(500*km), math.log(3000*km))),
1494 width=30., height=30.,
1495 show_grid=True,
1496 show_topo=True,
1497 color_dry=(238, 236, 230),
1498 topo_cpt_wet='light_sea_uniform',
1499 topo_cpt_dry='light_land_uniform',
1500 illuminate=True,
1501 illuminate_factor_ocean=0.15,
1502 show_rivers=False,
1503 show_plates=True)
1505 m.draw_cities()
1506 print(m)
1507 m.save('map_%02i.pdf' % i)