1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6import math
7import random
8import logging
10try:
11 from StringIO import StringIO as BytesIO
12except ImportError:
13 from io import BytesIO
15import numpy as num
17from pyrocko.guts import (Object, Float, Bool, Int, Tuple, String, List,
18 Unicode, Dict)
19from pyrocko.guts_array import Array
20from pyrocko.dataset import topo
21from pyrocko import orthodrome as od
22from .cpt import (
23 get_cpt_path, CPT, CPTLevel, CPTParseError, read_cpt, write_cpt,
24 cpt_merge_wet_dry)
25from . import gmtpy
26from . import nice_value
29CPT, CPTLevel, CPTParseError
31points_in_region = od.points_in_region
33logger = logging.getLogger('pyrocko.plot.automap')
35earthradius = 6371000.0
36r2d = 180./math.pi
37d2r = 1./r2d
38km = 1000.
39d2m = d2r*earthradius
40m2d = 1./d2m
41cm = gmtpy.cm
44def darken(c, f=0.7):
45 return (c[0]*f, c[1]*f, c[2]*f)
48def corners(lon, lat, w, h):
49 ll_lat, ll_lon = od.ne_to_latlon(lat, lon, -0.5*h, -0.5*w)
50 ur_lat, ur_lon = od.ne_to_latlon(lat, lon, 0.5*h, 0.5*w)
51 return ll_lon, ll_lat, ur_lon, ur_lat
54def extent(lon, lat, w, h, n):
55 x = num.linspace(-0.5*w, 0.5*w, n)
56 y = num.linspace(-0.5*h, 0.5*h, n)
57 slats, slons = od.ne_to_latlon(lat, lon, y[0], x)
58 nlats, nlons = od.ne_to_latlon(lat, lon, y[-1], x)
59 south = slats.min()
60 north = nlats.max()
62 wlats, wlons = od.ne_to_latlon(lat, lon, y, x[0])
63 elats, elons = od.ne_to_latlon(lat, lon, y, x[-1])
64 elons = num.where(elons < wlons, elons + 360., elons)
66 if elons.max() - elons.min() > 180 or wlons.max() - wlons.min() > 180.:
67 west = -180.
68 east = 180.
69 else:
70 west = wlons.min()
71 east = elons.max()
73 return topo.positive_region((west, east, south, north))
76class NoTopo(Exception):
77 pass
80class OutOfBounds(Exception):
81 pass
84class FloatTile(Object):
85 xmin = Float.T()
86 ymin = Float.T()
87 dx = Float.T()
88 dy = Float.T()
89 data = Array.T(shape=(None, None), dtype=float, serialize_as='table')
91 def __init__(self, xmin, ymin, dx, dy, data):
92 Object.__init__(self, init_props=False)
93 self.xmin = float(xmin)
94 self.ymin = float(ymin)
95 self.dx = float(dx)
96 self.dy = float(dy)
97 self.data = data
98 self._set_maxes()
100 def _set_maxes(self):
101 self.ny, self.nx = self.data.shape
102 self.xmax = self.xmin + (self.nx-1) * self.dx
103 self.ymax = self.ymin + (self.ny-1) * self.dy
105 def x(self):
106 return self.xmin + num.arange(self.nx) * self.dx
108 def y(self):
109 return self.ymin + num.arange(self.ny) * self.dy
111 def get(self, x, y):
112 ix = int(round((x - self.xmin) / self.dx))
113 iy = int(round((y - self.ymin) / self.dy))
114 if 0 <= ix < self.nx and 0 <= iy < self.ny:
115 return self.data[iy, ix]
116 else:
117 raise OutOfBounds()
120class City(Object):
121 def __init__(self, name, lat, lon, population=None, asciiname=None):
122 name = str(name)
123 lat = float(lat)
124 lon = float(lon)
125 if asciiname is None:
126 asciiname = name.encode('ascii', errors='replace')
128 if population is None:
129 population = 0
130 else:
131 population = int(population)
133 Object.__init__(self, name=name, lat=lat, lon=lon,
134 population=population, asciiname=asciiname)
136 name = Unicode.T()
137 lat = Float.T()
138 lon = Float.T()
139 population = Int.T()
140 asciiname = String.T()
143class Map(Object):
144 lat = Float.T(optional=True)
145 lon = Float.T(optional=True)
146 radius = Float.T(optional=True)
147 width = Float.T(default=20.)
148 height = Float.T(default=14.)
149 margins = List.T(Float.T())
150 illuminate = Bool.T(default=True)
151 skip_feature_factor = Float.T(default=0.02)
152 show_grid = Bool.T(default=False)
153 show_topo = Bool.T(default=True)
154 show_scale = Bool.T(default=False)
155 show_topo_scale = Bool.T(default=False)
156 show_center_mark = Bool.T(default=False)
157 show_rivers = Bool.T(default=True)
158 show_plates = Bool.T(default=False)
159 show_plate_velocities = Bool.T(default=False)
160 show_plate_names = Bool.T(default=False)
161 show_boundaries = Bool.T(default=False)
162 illuminate_factor_land = Float.T(default=0.5)
163 illuminate_factor_ocean = Float.T(default=0.25)
164 color_wet = Tuple.T(3, Int.T(), default=(216, 242, 254))
165 color_dry = Tuple.T(3, Int.T(), default=(172, 208, 165))
166 color_boundaries = Tuple.T(3, Int.T(), default=(1, 1, 1))
167 topo_resolution_min = Float.T(
168 default=40.,
169 help='minimum resolution of topography [dpi]')
170 topo_resolution_max = Float.T(
171 default=200.,
172 help='maximum resolution of topography [dpi]')
173 replace_topo_color_only = FloatTile.T(
174 optional=True,
175 help='replace topo color while keeping topographic shading')
176 topo_cpt_wet = String.T(default='light_sea')
177 topo_cpt_dry = String.T(default='light_land')
178 axes_layout = String.T(optional=True)
179 custom_cities = List.T(City.T())
180 gmt_config = Dict.T(String.T(), String.T())
181 comment = String.T(optional=True)
182 approx_ticks = Int.T(default=4)
184 def __init__(self, gmtversion='newest', **kwargs):
185 Object.__init__(self, **kwargs)
186 self._gmt = None
187 self._scaler = None
188 self._widget = None
189 self._corners = None
190 self._wesn = None
191 self._minarea = None
192 self._coastline_resolution = None
193 self._rivers = None
194 self._dems = None
195 self._have_topo_land = None
196 self._have_topo_ocean = None
197 self._jxyr = None
198 self._prep_topo_have = None
199 self._labels = []
200 self._area_labels = []
201 self._gmtversion = gmtversion
203 def save(self, outpath, resolution=75., oversample=2., size=None,
204 width=None, height=None, psconvert=False, crop_eps_mode=False):
206 '''
207 Save the image.
209 Save the image to ``outpath``. The format is determined by the filename
210 extension. Formats are handled as follows: ``'.eps'`` and ``'.ps'``
211 produce EPS and PS, respectively, directly with GMT. If the file name
212 ends with ``'.pdf'``, GMT output is fed through ``gmtpy-epstopdf`` to
213 create a PDF file. For any other filename extension, output is first
214 converted to PDF with ``gmtpy-epstopdf``, then with ``pdftocairo`` to
215 PNG with a resolution oversampled by the factor ``oversample`` and
216 finally the PNG is downsampled and converted to the target format with
217 ``convert``. The resolution of rasterized target image can be
218 controlled either by ``resolution`` in DPI or by specifying ``width``
219 or ``height`` or ``size``, where the latter fits the image into a
220 square with given side length. To save transparency use
221 ``psconvert=True``. To crop the output image with a rectangle to the
222 nearest non-white element set ``crop_eps_mode=True``.
223 '''
225 gmt = self.gmt
226 self.draw_labels()
227 self.draw_axes()
228 if self.show_topo and self.show_topo_scale:
229 self._draw_topo_scale()
231 gmt.save(outpath, resolution=resolution, oversample=oversample,
232 crop_eps_mode=crop_eps_mode,
233 size=size, width=width, height=height, psconvert=psconvert)
235 @property
236 def scaler(self):
237 if self._scaler is None:
238 self._setup_geometry()
240 return self._scaler
242 @property
243 def wesn(self):
244 if self._wesn is None:
245 self._setup_geometry()
247 return self._wesn
249 @property
250 def widget(self):
251 if self._widget is None:
252 self._setup()
254 return self._widget
256 @property
257 def layout(self):
258 if self._layout is None:
259 self._setup()
261 return self._layout
263 @property
264 def jxyr(self):
265 if self._jxyr is None:
266 self._setup()
268 return self._jxyr
270 @property
271 def pxyr(self):
272 if self._pxyr is None:
273 self._setup()
275 return self._pxyr
277 @property
278 def gmt(self):
279 if self._gmt is None:
280 self._setup()
282 if self._have_topo_ocean is None:
283 self._draw_background()
285 return self._gmt
287 def _setup(self):
288 if not self._widget:
289 self._setup_geometry()
291 self._setup_lod()
292 self._setup_gmt()
294 def _setup_geometry(self):
295 wpage, hpage = self.width, self.height
296 ml, mr, mt, mb = self._expand_margins()
297 wpage -= ml + mr
298 hpage -= mt + mb
300 wreg = self.radius * 2.0
301 hreg = self.radius * 2.0
302 if wpage >= hpage:
303 wreg *= wpage/hpage
304 else:
305 hreg *= hpage/wpage
307 self._wreg = wreg
308 self._hreg = hreg
310 self._corners = corners(self.lon, self.lat, wreg, hreg)
311 west, east, south, north = extent(self.lon, self.lat, wreg, hreg, 10)
313 x, y, z = ((west, east), (south, north), (-6000., 4500.))
315 xax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
316 yax = gmtpy.Ax(mode='min-max', approx_ticks=self.approx_ticks)
317 zax = gmtpy.Ax(mode='min-max', inc=1000., label='Height',
318 scaled_unit='km', scaled_unit_factor=0.001)
320 scaler = gmtpy.ScaleGuru(data_tuples=[(x, y, z)], axes=(xax, yax, zax))
322 par = scaler.get_params()
324 west = par['xmin']
325 east = par['xmax']
326 south = par['ymin']
327 north = par['ymax']
329 self._wesn = west, east, south, north
330 self._scaler = scaler
332 def _setup_lod(self):
333 w, e, s, n = self._wesn
334 if self.radius > 1500.*km:
335 coastline_resolution = 'i'
336 rivers = False
337 else:
338 coastline_resolution = 'f'
339 rivers = True
341 self._minarea = (self.skip_feature_factor * self.radius/km)**2
343 self._coastline_resolution = coastline_resolution
344 self._rivers = rivers
346 self._prep_topo_have = {}
347 self._dems = {}
349 cm2inch = gmtpy.cm/gmtpy.inch
351 dmin = 2.0 * self.radius * m2d / (self.topo_resolution_max *
352 (self.height * cm2inch))
353 dmax = 2.0 * self.radius * m2d / (self.topo_resolution_min *
354 (self.height * cm2inch))
356 for k in ['ocean', 'land']:
357 self._dems[k] = topo.select_dem_names(k, dmin, dmax, self._wesn)
358 if self._dems[k]:
359 logger.debug('using topography dataset %s for %s'
360 % (','.join(self._dems[k]), k))
362 def _expand_margins(self):
363 if len(self.margins) == 0 or len(self.margins) > 4:
364 ml = mr = mt = mb = 2.0
365 elif len(self.margins) == 1:
366 ml = mr = mt = mb = self.margins[0]
367 elif len(self.margins) == 2:
368 ml = mr = self.margins[0]
369 mt = mb = self.margins[1]
370 elif len(self.margins) == 4:
371 ml, mr, mt, mb = self.margins
373 return ml, mr, mt, mb
375 def _setup_gmt(self):
376 w, h = self.width, self.height
377 scaler = self._scaler
379 if gmtpy.is_gmt5(self._gmtversion):
380 gmtconf = dict(
381 MAP_TICK_PEN_PRIMARY='1.25p',
382 MAP_TICK_PEN_SECONDARY='1.25p',
383 MAP_TICK_LENGTH_PRIMARY='0.2c',
384 MAP_TICK_LENGTH_SECONDARY='0.6c',
385 FONT_ANNOT_PRIMARY='12p,1,black',
386 FONT_LABEL='12p,1,black',
387 PS_CHAR_ENCODING='ISOLatin1+',
388 MAP_FRAME_TYPE='fancy',
389 FORMAT_GEO_MAP='D',
390 PS_MEDIA='Custom_%ix%i' % (
391 w*gmtpy.cm,
392 h*gmtpy.cm),
393 PS_PAGE_ORIENTATION='portrait',
394 MAP_GRID_PEN_PRIMARY='thinnest,0/50/0',
395 MAP_ANNOT_OBLIQUE='6')
396 else:
397 gmtconf = dict(
398 TICK_PEN='1.25p',
399 TICK_LENGTH='0.2c',
400 ANNOT_FONT_PRIMARY='1',
401 ANNOT_FONT_SIZE_PRIMARY='12p',
402 LABEL_FONT='1',
403 LABEL_FONT_SIZE='12p',
404 CHAR_ENCODING='ISOLatin1+',
405 BASEMAP_TYPE='fancy',
406 PLOT_DEGREE_FORMAT='D',
407 PAPER_MEDIA='Custom_%ix%i' % (
408 w*gmtpy.cm,
409 h*gmtpy.cm),
410 GRID_PEN_PRIMARY='thinnest/0/50/0',
411 DOTS_PR_INCH='1200',
412 OBLIQUE_ANNOTATION='6')
414 gmtconf.update(
415 (k.upper(), v) for (k, v) in self.gmt_config.items())
417 gmt = gmtpy.GMT(config=gmtconf, version=self._gmtversion)
419 layout = gmt.default_layout()
421 layout.set_fixed_margins(*[x*cm for x in self._expand_margins()])
423 widget = layout.get_widget()
424 widget['P'] = widget['J']
425 widget['J'] = ('-JA%g/%g' % (self.lon, self.lat)) + '/%(width)gp'
426 scaler['R'] = '-R%g/%g/%g/%gr' % self._corners
428 # aspect = gmtpy.aspect_for_projection(
429 # gmt.installation['version'], *(widget.J() + scaler.R()))
431 aspect = self._map_aspect(jr=widget.J() + scaler.R())
432 widget.set_aspect(aspect)
434 self._gmt = gmt
435 self._layout = layout
436 self._widget = widget
437 self._jxyr = self._widget.JXY() + self._scaler.R()
438 self._pxyr = self._widget.PXY() + [
439 '-R%g/%g/%g/%g' % (0, widget.width(), 0, widget.height())]
440 self._have_drawn_axes = False
441 self._have_drawn_labels = False
443 def _draw_background(self):
444 self._have_topo_land = False
445 self._have_topo_ocean = False
446 if self.show_topo:
447 self._have_topo = self._draw_topo()
449 self._draw_basefeatures()
451 def _get_topo_tile(self, k):
452 t = None
453 demname = None
454 for dem in self._dems[k]:
455 t = topo.get(dem, self._wesn)
456 demname = dem
457 if t is not None:
458 break
460 if not t:
461 raise NoTopo()
463 return t, demname
465 def _prep_topo(self, k):
466 gmt = self._gmt
467 t, demname = self._get_topo_tile(k)
469 if demname not in self._prep_topo_have:
471 grdfile = gmt.tempfilename()
473 is_flat = num.all(t.data[0] == t.data)
475 gmtpy.savegrd(
476 t.x(), t.y(), t.data, filename=grdfile, naming='lonlat')
478 if self.illuminate and not is_flat:
479 if k == 'ocean':
480 factor = self.illuminate_factor_ocean
481 else:
482 factor = self.illuminate_factor_land
484 ilumfn = gmt.tempfilename()
485 gmt.grdgradient(
486 grdfile,
487 N='e%g' % factor,
488 A=-45,
489 G=ilumfn,
490 out_discard=True)
492 ilumargs = ['-I%s' % ilumfn]
493 else:
494 ilumargs = []
496 if self.replace_topo_color_only:
497 t2 = self.replace_topo_color_only
498 grdfile2 = gmt.tempfilename()
500 gmtpy.savegrd(
501 t2.x(), t2.y(), t2.data, filename=grdfile2,
502 naming='lonlat')
504 if gmt.is_gmt5():
505 gmt.grdsample(
506 grdfile2,
507 G=grdfile,
508 n='l',
509 I='%g/%g' % (t.dx, t.dy), # noqa
510 R=grdfile,
511 out_discard=True)
512 else:
513 gmt.grdsample(
514 grdfile2,
515 G=grdfile,
516 Q='l',
517 I='%g/%g' % (t.dx, t.dy), # noqa
518 R=grdfile,
519 out_discard=True)
521 gmt.grdmath(
522 grdfile, '0.0', 'AND', '=', grdfile2,
523 out_discard=True)
525 grdfile = grdfile2
527 self._prep_topo_have[demname] = grdfile, ilumargs
529 return self._prep_topo_have[demname]
531 def _draw_topo(self):
532 widget = self._widget
533 scaler = self._scaler
534 gmt = self._gmt
535 cres = self._coastline_resolution
536 minarea = self._minarea
538 JXY = widget.JXY()
539 R = scaler.R()
541 try:
542 grdfile, ilumargs = self._prep_topo('ocean')
543 gmt.pscoast(D=cres, S='c', A=minarea, *(JXY+R))
544 gmt.grdimage(grdfile, C=get_cpt_path(self.topo_cpt_wet),
545 *(ilumargs+JXY+R))
546 gmt.pscoast(Q=True, *(JXY+R))
547 self._have_topo_ocean = True
548 except NoTopo:
549 self._have_topo_ocean = False
551 try:
552 grdfile, ilumargs = self._prep_topo('land')
553 gmt.pscoast(D=cres, G='c', A=minarea, *(JXY+R))
554 gmt.grdimage(grdfile, C=get_cpt_path(self.topo_cpt_dry),
555 *(ilumargs+JXY+R))
556 gmt.pscoast(Q=True, *(JXY+R))
557 self._have_topo_land = True
558 except NoTopo:
559 self._have_topo_land = False
561 def _draw_topo_scale(self, label='Elevation [km]'):
562 dry = read_cpt(get_cpt_path(self.topo_cpt_dry))
563 wet = read_cpt(get_cpt_path(self.topo_cpt_wet))
564 combi = cpt_merge_wet_dry(wet, dry)
565 for level in combi.levels:
566 level.vmin /= km
567 level.vmax /= km
569 topo_cpt = self.gmt.tempfilename() + '.cpt'
570 write_cpt(combi, topo_cpt)
572 (w, h), (xo, yo) = self.widget.get_size()
573 self.gmt.psscale(
574 D='%gp/%gp/%gp/%gph' % (xo + 0.5*w, yo - 2.0*gmtpy.cm, w,
575 0.5*gmtpy.cm),
576 C=topo_cpt,
577 B='1:%s:' % label)
579 def _draw_basefeatures(self):
580 gmt = self._gmt
581 cres = self._coastline_resolution
582 rivers = self._rivers
583 minarea = self._minarea
585 color_wet = self.color_wet
586 color_dry = self.color_dry
588 if self.show_rivers and rivers:
589 rivers = ['-Ir/0.25p,%s' % gmtpy.color(self.color_wet)]
590 else:
591 rivers = []
593 fill = {}
594 if not self._have_topo_land:
595 fill['G'] = color_dry
597 if not self._have_topo_ocean:
598 fill['S'] = color_wet
600 if self.show_boundaries:
601 fill['N'] = '1/1p,%s,%s' % (
602 gmtpy.color(self.color_boundaries), 'solid')
604 gmt.pscoast(
605 D=cres,
606 W='thinnest,%s' % gmtpy.color(darken(gmtpy.color_tup(color_dry))),
607 A=minarea,
608 *(rivers+self._jxyr), **fill)
610 if self.show_plates:
611 self.draw_plates()
613 def _draw_axes(self):
614 gmt = self._gmt
615 scaler = self._scaler
616 widget = self._widget
618 if self.axes_layout is None:
619 if self.lat > 0.0:
620 axes_layout = 'WSen'
621 else:
622 axes_layout = 'WseN'
623 else:
624 axes_layout = self.axes_layout
626 scale_km = nice_value(self.radius/5.) / 1000.
628 if self.show_center_mark:
629 gmt.psxy(
630 in_rows=[[self.lon, self.lat]],
631 S='c20p', W='2p,black',
632 *self._jxyr)
634 if self.show_grid:
635 btmpl = ('%(xinc)gg%(xinc)g:%(xlabel)s:/'
636 '%(yinc)gg%(yinc)g:%(ylabel)s:')
637 else:
638 btmpl = '%(xinc)g:%(xlabel)s:/%(yinc)g:%(ylabel)s:'
640 if self.show_scale:
641 scale = 'x%gp/%gp/%g/%g/%gk' % (
642 6./7*widget.width(),
643 widget.height()/7.,
644 self.lon,
645 self.lat,
646 scale_km)
647 else:
648 scale = False
650 gmt.psbasemap(
651 B=(btmpl % scaler.get_params())+axes_layout,
652 L=scale,
653 *self._jxyr)
655 if self.comment:
656 font_size = self.gmt.label_font_size()
658 _, east, south, _ = self._wesn
659 if gmt.is_gmt5():
660 row = [
661 1, 0,
662 '%gp,%s,%s' % (font_size, 0, 'black'), 'BR',
663 self.comment]
665 farg = ['-F+f+j']
666 else:
667 row = [1, 0, font_size, 0, 0, 'BR', self.comment]
668 farg = []
670 gmt.pstext(
671 in_rows=[row],
672 N=True,
673 R=(0, 1, 0, 1),
674 D='%gp/%gp' % (-font_size*0.2, font_size*0.3),
675 *(widget.PXY() + farg))
677 def draw_axes(self):
678 if not self._have_drawn_axes:
679 self._draw_axes()
680 self._have_drawn_axes = True
682 def _have_coastlines(self):
683 gmt = self._gmt
684 cres = self._coastline_resolution
685 minarea = self._minarea
687 checkfile = gmt.tempfilename()
689 gmt.pscoast(
690 M=True,
691 D=cres,
692 W='thinnest,black',
693 A=minarea,
694 out_filename=checkfile,
695 *self._jxyr)
697 points = []
698 with open(checkfile, 'r') as f:
699 for line in f:
700 ls = line.strip()
701 if ls.startswith('#') or ls.startswith('>') or ls == '':
702 continue
703 plon, plat = [float(x) for x in ls.split()]
704 points.append((plat, plon))
706 points = num.array(points, dtype=float)
707 return num.any(points_in_region(points, self._wesn))
709 def have_coastlines(self):
710 self.gmt
711 return self._have_coastlines()
713 def project(self, lats, lons, jr=None):
714 onepoint = False
715 if isinstance(lats, float) and isinstance(lons, float):
716 lats = [lats]
717 lons = [lons]
718 onepoint = True
720 if jr is not None:
721 j, r = jr
722 gmt = gmtpy.GMT(version=self._gmtversion)
723 else:
724 j, _, _, r = self.jxyr
725 gmt = self.gmt
727 f = BytesIO()
728 gmt.mapproject(j, r, in_columns=(lons, lats), out_stream=f, D='p')
729 f.seek(0)
730 data = num.loadtxt(f, ndmin=2)
731 xs, ys = data.T
732 if onepoint:
733 xs = xs[0]
734 ys = ys[0]
735 return xs, ys
737 def _map_box(self, jr=None):
738 ll_lon, ll_lat, ur_lon, ur_lat = self._corners
740 xs_corner, ys_corner = self.project(
741 (ll_lat, ur_lat), (ll_lon, ur_lon), jr=jr)
743 w = xs_corner[1] - xs_corner[0]
744 h = ys_corner[1] - ys_corner[0]
746 return w, h
748 def _map_aspect(self, jr=None):
749 w, h = self._map_box(jr=jr)
750 return h/w
752 def _draw_labels(self):
753 points_taken = []
754 regions_taken = []
756 def no_points_in_rect(xs, ys, xmin, ymin, xmax, ymax):
757 xx = not num.any(la(la(xmin < xs, xs < xmax),
758 la(ymin < ys, ys < ymax)))
759 return xx
761 def roverlaps(a, b):
762 return (a[0] < b[2] and b[0] < a[2] and
763 a[1] < b[3] and b[1] < a[3])
765 w, h = self._map_box()
767 label_font_size = self.gmt.label_font_size()
769 if self._labels:
771 n = len(self._labels)
773 lons, lats, texts, sx, sy, colors, fonts, font_sizes, \
774 angles, styles = list(zip(*self._labels))
776 font_sizes = [
777 (font_size or label_font_size) for font_size in font_sizes]
779 sx = num.array(sx, dtype=float)
780 sy = num.array(sy, dtype=float)
782 xs, ys = self.project(lats, lons)
784 points_taken.append((xs, ys))
786 dxs = num.zeros(n)
787 dys = num.zeros(n)
789 for i in range(n):
790 dx, dy = gmtpy.text_box(
791 texts[i],
792 font=fonts[i],
793 font_size=font_sizes[i],
794 **styles[i])
796 dxs[i] = dx
797 dys[i] = dy
799 la = num.logical_and
800 anchors_ok = (
801 la(xs + sx + dxs < w, ys + sy + dys < h),
802 la(xs - sx - dxs > 0., ys - sy - dys > 0.),
803 la(xs + sx + dxs < w, ys - sy - dys > 0.),
804 la(xs - sx - dxs > 0., ys + sy + dys < h),
805 )
807 arects = [
808 (xs, ys, xs + sx + dxs, ys + sy + dys),
809 (xs - sx - dxs, ys - sy - dys, xs, ys),
810 (xs, ys - sy - dys, xs + sx + dxs, ys),
811 (xs - sx - dxs, ys, xs, ys + sy + dys)]
813 for i in range(n):
814 for ianch in range(4):
815 anchors_ok[ianch][i] &= no_points_in_rect(
816 xs, ys, *[xxx[i] for xxx in arects[ianch]])
818 anchor_choices = []
819 anchor_take = []
820 for i in range(n):
821 choices = [ianch for ianch in range(4)
822 if anchors_ok[ianch][i]]
823 anchor_choices.append(choices)
824 if choices:
825 anchor_take.append(choices[0])
826 else:
827 anchor_take.append(None)
829 def cost(anchor_take):
830 noverlaps = 0
831 for i in range(n):
832 for j in range(n):
833 if i != j:
834 i_take = anchor_take[i]
835 j_take = anchor_take[j]
836 if i_take is None or j_take is None:
837 continue
838 r_i = [xxx[i] for xxx in arects[i_take]]
839 r_j = [xxx[j] for xxx in arects[j_take]]
840 if roverlaps(r_i, r_j):
841 noverlaps += 1
843 return noverlaps
845 cur_cost = cost(anchor_take)
846 imax = 30
847 while cur_cost != 0 and imax > 0:
848 for i in range(n):
849 for t in anchor_choices[i]:
850 anchor_take_new = list(anchor_take)
851 anchor_take_new[i] = t
852 new_cost = cost(anchor_take_new)
853 if new_cost < cur_cost:
854 anchor_take = anchor_take_new
855 cur_cost = new_cost
857 imax -= 1
859 while cur_cost != 0:
860 for i in range(n):
861 anchor_take_new = list(anchor_take)
862 anchor_take_new[i] = None
863 new_cost = cost(anchor_take_new)
864 if new_cost < cur_cost:
865 anchor_take = anchor_take_new
866 cur_cost = new_cost
867 break
869 anchor_strs = ['BL', 'TR', 'TL', 'BR']
871 for i in range(n):
872 ianchor = anchor_take[i]
873 color = colors[i]
874 if color is None:
875 color = 'black'
877 if ianchor is not None:
878 regions_taken.append([xxx[i] for xxx in arects[ianchor]])
880 anchor = anchor_strs[ianchor]
882 yoff = [-sy[i], sy[i]][anchor[0] == 'B']
883 xoff = [-sx[i], sx[i]][anchor[1] == 'L']
884 if self.gmt.is_gmt5():
885 row = (
886 lons[i], lats[i],
887 '%i,%s,%s' % (font_sizes[i], fonts[i], color),
888 anchor,
889 texts[i])
891 farg = ['-F+f+j+a%g' % angles[i]]
892 else:
893 row = (
894 lons[i], lats[i],
895 font_sizes[i], angles[i], fonts[i], anchor,
896 texts[i])
897 farg = ['-G%s' % color]
899 self.gmt.pstext(
900 in_rows=[row],
901 D='%gp/%gp' % (xoff, yoff),
902 *(self.jxyr + farg),
903 **styles[i])
905 if self._area_labels:
907 for lons, lats, text, color, font, font_size, style in \
908 self._area_labels:
910 if font_size is None:
911 font_size = label_font_size
913 if color is None:
914 color = 'black'
916 if self.gmt.is_gmt5():
917 farg = ['-F+f+j']
918 else:
919 farg = ['-G%s' % color]
921 xs, ys = self.project(lats, lons)
922 dx, dy = gmtpy.text_box(
923 text, font=font, font_size=font_size, **style)
925 rects = [xs-0.5*dx, ys-0.5*dy, xs+0.5*dx, ys+0.5*dy]
927 locs_ok = num.ones(xs.size, dtype=bool)
929 for iloc in range(xs.size):
930 rcandi = [xxx[iloc] for xxx in rects]
932 locs_ok[iloc] = True
933 locs_ok[iloc] &= (
934 0 < rcandi[0] and rcandi[2] < w
935 and 0 < rcandi[1] and rcandi[3] < h)
937 overlap = False
938 for r in regions_taken:
939 if roverlaps(r, rcandi):
940 overlap = True
941 break
943 locs_ok[iloc] &= not overlap
945 for xs_taken, ys_taken in points_taken:
946 locs_ok[iloc] &= no_points_in_rect(
947 xs_taken, ys_taken, *rcandi)
949 if not locs_ok[iloc]:
950 break
952 rows = []
953 for iloc, (lon, lat) in enumerate(zip(lons, lats)):
954 if not locs_ok[iloc]:
955 continue
957 if self.gmt.is_gmt5():
958 row = (
959 lon, lat,
960 '%i,%s,%s' % (font_size, font, color),
961 'MC',
962 text)
964 else:
965 row = (
966 lon, lat,
967 font_size, 0, font, 'MC',
968 text)
970 rows.append(row)
972 regions_taken.append([xxx[iloc] for xxx in rects])
973 break
975 self.gmt.pstext(
976 in_rows=rows,
977 *(self.jxyr + farg),
978 **style)
980 def draw_labels(self):
981 self.gmt
982 if not self._have_drawn_labels:
983 self._draw_labels()
984 self._have_drawn_labels = True
986 def add_label(
987 self, lat, lon, text,
988 offset_x=5., offset_y=5.,
989 color=None,
990 font='1',
991 font_size=None,
992 angle=0,
993 style={}):
995 if 'G' in style:
996 style = style.copy()
997 color = style.pop('G')
999 self._labels.append(
1000 (lon, lat, text, offset_x, offset_y, color, font, font_size,
1001 angle, style))
1003 def add_area_label(
1004 self, lat, lon, text,
1005 color=None,
1006 font='3',
1007 font_size=None,
1008 style={}):
1010 self._area_labels.append(
1011 (lon, lat, text, color, font, font_size, style))
1013 def cities_in_region(self):
1014 from pyrocko.dataset import geonames
1015 cities = geonames.get_cities_region(region=self.wesn, minpop=0)
1016 cities.extend(self.custom_cities)
1017 cities.sort(key=lambda x: x.population)
1018 return cities
1020 def draw_cities(self,
1021 exact=None,
1022 include=[],
1023 exclude=[],
1024 nmax_soft=10,
1025 psxy_style=dict(S='s5p', G='black')):
1027 cities = self.cities_in_region()
1029 if exact is not None:
1030 cities = [c for c in cities if c.name in exact]
1031 minpop = None
1032 else:
1033 cities = [c for c in cities if c.name not in exclude]
1034 minpop = 10**3
1035 for minpop_new in [1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6, 3e6, 1e7]:
1036 cities_new = [
1037 c for c in cities
1038 if c.population > minpop_new or c.name in include]
1040 if len(cities_new) == 0 or (
1041 len(cities_new) < 3 and len(cities) < nmax_soft*2):
1042 break
1044 cities = cities_new
1045 minpop = minpop_new
1046 if len(cities) <= nmax_soft:
1047 break
1049 if cities:
1050 lats = [c.lat for c in cities]
1051 lons = [c.lon for c in cities]
1053 self.gmt.psxy(
1054 in_columns=(lons, lats),
1055 *self.jxyr, **psxy_style)
1057 for c in cities:
1058 try:
1059 text = c.name.encode('iso-8859-1').decode('iso-8859-1')
1060 except UnicodeEncodeError:
1061 text = c.asciiname
1063 self.add_label(c.lat, c.lon, text)
1065 self._cities_minpop = minpop
1067 def add_stations(self, stations, psxy_style=dict()):
1069 default_psxy_style = {
1070 'S': 't8p',
1071 'G': 'black'
1072 }
1073 default_psxy_style.update(psxy_style)
1075 lats, lons = zip(*[s.effective_latlon for s in stations])
1077 self.gmt.psxy(
1078 in_columns=(lons, lats),
1079 *self.jxyr, **default_psxy_style)
1081 for station in stations:
1082 self.add_label(
1083 station.effective_lat,
1084 station.effective_lon,
1085 '.'.join(x for x in (station.network, station.station) if x))
1087 def add_kite_scene(self, scene):
1088 tile = FloatTile(
1089 scene.frame.llLon,
1090 scene.frame.llLat,
1091 scene.frame.dLon,
1092 scene.frame.dLat,
1093 scene.displacement)
1095 return tile
1097 def add_gnss_campaign(self, campaign, psxy_style=None, offset_scale=None,
1098 labels=True, vertical=False, fontsize=10):
1100 stations = campaign.stations
1102 if offset_scale is None:
1103 offset_scale = num.zeros(campaign.nstations)
1104 for ista, sta in enumerate(stations):
1105 for comp in sta.components.values():
1106 offset_scale[ista] += comp.shift
1107 offset_scale = num.sqrt(offset_scale**2).max()
1109 size = math.sqrt(self.height**2 + self.width**2)
1110 scale = (size/10.) / offset_scale
1111 logger.debug('GNSS: Using offset scale %f, map scale %f',
1112 offset_scale, scale)
1114 lats, lons = zip(*[s.effective_latlon for s in stations])
1116 if self.gmt.is_gmt6():
1117 sign_factor = 1.
1118 arrow_head_placement = 'e'
1119 else:
1120 sign_factor = -1.
1121 arrow_head_placement = 'b'
1123 if vertical:
1124 rows = [[lons[ista], lats[ista],
1125 0., sign_factor * s.up.shift,
1126 (s.east.sigma + s.north.sigma) if s.east.sigma else 0.,
1127 s.up.sigma, 0.,
1128 s.code if labels else None]
1129 for ista, s in enumerate(stations)
1130 if s.up is not None]
1132 else:
1133 rows = [[lons[ista], lats[ista],
1134 sign_factor * s.east.shift, sign_factor * s.north.shift,
1135 s.east.sigma, s.north.sigma, s.correlation_ne,
1136 s.code if labels else None]
1137 for ista, s in enumerate(stations)
1138 if s.east is not None or s.north is not None]
1140 default_psxy_style = {
1141 'h': 0,
1142 'W': '2p,black',
1143 'A': '+p2p,black+{}+a40'.format(arrow_head_placement),
1144 'G': 'black',
1145 'L': True,
1146 'S': 'e%dc/0.95/%d' % (scale, fontsize),
1147 }
1149 if not labels:
1150 for row in rows:
1151 row.pop(-1)
1153 if psxy_style is not None:
1154 default_psxy_style.update(psxy_style)
1156 self.gmt.psvelo(
1157 in_rows=rows,
1158 *self.jxyr,
1159 **default_psxy_style)
1161 def draw_plates(self):
1162 from pyrocko.dataset import tectonics
1164 neast = 20
1165 nnorth = max(1, int(round(num.round(self._hreg/self._wreg * neast))))
1166 norths = num.linspace(-self._hreg*0.5, self._hreg*0.5, nnorth)
1167 easts = num.linspace(-self._wreg*0.5, self._wreg*0.5, neast)
1168 norths2 = num.repeat(norths, neast)
1169 easts2 = num.tile(easts, nnorth)
1170 lats, lons = od.ne_to_latlon(
1171 self.lat, self.lon, norths2, easts2)
1173 bird = tectonics.PeterBird2003()
1174 plates = bird.get_plates()
1176 color_plates = gmtpy.color('aluminium5')
1177 color_velocities = gmtpy.color('skyblue1')
1178 color_velocities_lab = gmtpy.color(darken(gmtpy.color_tup('skyblue1')))
1180 points = num.vstack((lats, lons)).T
1181 used = []
1182 for plate in plates:
1183 mask = plate.contains_points(points)
1184 if num.any(mask):
1185 used.append((plate, mask))
1187 if len(used) > 1:
1189 candi_fixed = {}
1191 label_data = []
1192 for plate, mask in used:
1194 mean_north = num.mean(norths2[mask])
1195 mean_east = num.mean(easts2[mask])
1196 iorder = num.argsort(num.sqrt(
1197 (norths2[mask] - mean_north)**2 +
1198 (easts2[mask] - mean_east)**2))
1200 lat_candis = lats[mask][iorder]
1201 lon_candis = lons[mask][iorder]
1203 candi_fixed[plate.name] = lat_candis.size
1205 label_data.append((
1206 lat_candis, lon_candis, plate, color_plates))
1208 boundaries = bird.get_boundaries()
1210 size = 1.
1212 psxy_kwargs = []
1214 for boundary in boundaries:
1215 if num.any(points_in_region(boundary.points, self._wesn)):
1216 for typ, part in boundary.split_types(
1217 [['SUB'],
1218 ['OSR', 'CRB'],
1219 ['OTF', 'CTF', 'OCB', 'CCB']]):
1221 lats, lons = part.T
1223 kwargs = {}
1225 kwargs['in_columns'] = (lons, lats)
1226 kwargs['W'] = '%gp,%s' % (size, color_plates)
1228 if typ[0] == 'SUB':
1229 if boundary.kind == '\\':
1230 kwargs['S'] = 'f%g/%gp+t+r' % (
1231 0.45*size, 3.*size)
1232 elif boundary.kind == '/':
1233 kwargs['S'] = 'f%g/%gp+t+l' % (
1234 0.45*size, 3.*size)
1236 kwargs['G'] = color_plates
1238 elif typ[0] in ['OSR', 'CRB']:
1239 kwargs_bg = {}
1240 kwargs_bg['in_columns'] = (lons, lats)
1241 kwargs_bg['W'] = '%gp,%s' % (
1242 size * 3, color_plates)
1243 psxy_kwargs.append(kwargs_bg)
1245 kwargs['W'] = '%gp,%s' % (size * 2, 'white')
1247 psxy_kwargs.append(kwargs)
1249 if boundary.kind == '\\':
1250 if boundary.plate_name2 in candi_fixed:
1251 candi_fixed[boundary.plate_name2] += \
1252 neast*nnorth
1254 elif boundary.kind == '/':
1255 if boundary.plate_name1 in candi_fixed:
1256 candi_fixed[boundary.plate_name1] += \
1257 neast*nnorth
1259 candi_fixed = [name for name in sorted(
1260 list(candi_fixed.keys()), key=lambda name: -candi_fixed[name])]
1262 candi_fixed.append(None)
1264 gsrm = tectonics.GSRM1()
1266 for name in candi_fixed:
1267 if name not in gsrm.plate_names() \
1268 and name not in gsrm.plate_alt_names():
1270 continue
1272 lats, lons, vnorth, veast, vnorth_err, veast_err, corr = \
1273 gsrm.get_velocities(name, region=self._wesn)
1275 fixed_plate_name = name
1277 if self.show_plate_velocities:
1278 self.gmt.psvelo(
1279 in_columns=(
1280 lons, lats, veast, vnorth, veast_err, vnorth_err,
1281 corr),
1282 W='0.25p,%s' % color_velocities,
1283 A='9p+e+g%s' % color_velocities,
1284 S='e0.2p/0.95/10',
1285 *self.jxyr)
1287 for _ in range(len(lons) // 50 + 1):
1288 ii = random.randint(0, len(lons)-1)
1289 v = math.sqrt(vnorth[ii]**2 + veast[ii]**2)
1290 self.add_label(
1291 lats[ii], lons[ii], '%.0f' % v,
1292 font_size=0.7*self.gmt.label_font_size(),
1293 style=dict(
1294 G=color_velocities_lab))
1296 break
1298 if self.show_plate_names:
1299 for (lat_candis, lon_candis, plate, color) in label_data:
1300 full_name = bird.full_name(plate.name)
1301 if plate.name == fixed_plate_name:
1302 full_name = '@_' + full_name + '@_'
1304 self.add_area_label(
1305 lat_candis, lon_candis,
1306 full_name,
1307 color=color,
1308 font='3')
1310 for kwargs in psxy_kwargs:
1311 self.gmt.psxy(*self.jxyr, **kwargs)
1314def rand(mi, ma):
1315 mi = float(mi)
1316 ma = float(ma)
1317 return random.random() * (ma-mi) + mi
1320def split_region(region):
1321 west, east, south, north = topo.positive_region(region)
1322 if east > 180:
1323 return [(west, 180., south, north),
1324 (-180., east-360., south, north)]
1325 else:
1326 return [region]
1329if __name__ == '__main__':
1330 from pyrocko import util
1331 util.setup_logging('pyrocko.automap', 'info')
1333 import sys
1334 if len(sys.argv) == 2:
1336 n = int(sys.argv[1])
1338 for i in range(n):
1339 m = Map(
1340 lat=rand(-60., 60.),
1341 lon=rand(-180., 180.),
1342 radius=math.exp(rand(math.log(500*km), math.log(3000*km))),
1343 width=30., height=30.,
1344 show_grid=True,
1345 show_topo=True,
1346 color_dry=(238, 236, 230),
1347 topo_cpt_wet='light_sea_uniform',
1348 topo_cpt_dry='light_land_uniform',
1349 illuminate=True,
1350 illuminate_factor_ocean=0.15,
1351 show_rivers=False,
1352 show_plates=True)
1354 m.draw_cities()
1355 print(m)
1356 m.save('map_%02i.pdf' % i)