1#!/usr/bin/env python
2# http://pyrocko.org - GPLv3
3#
4# The Pyrocko Developers, 21st Century
5# ---|P------/S----------~Lg----------
7import sys
8import logging
9import numpy as num
10from pyrocko import util, model, orthodrome as od, gmtpy
11from pyrocko import moment_tensor as pmt
12from pyrocko.plot import automap
14from optparse import OptionParser
16km = 1000.
17logger = logging.getLogger('pyrocko.apps.automap')
19program_name = 'automap'
21usage = '''
22usage: %s [options] [--] <lat> <lon> <radius_km> <output.(pdf|png)>
23 %s [--download-etopo1] [--download-srtmgl3]
24'''.strip() % (program_name, program_name)
26description = '''Create a simple map with topography.'''
29def latlon_arrays(locs):
30 return num.array(
31 [(x.effective_lat, x.effective_lon) for x in locs]).T
34def main(args=None):
35 if args is None:
36 args = sys.argv[1:]
38 parser = OptionParser(
39 usage=usage,
40 description=description)
42 parser.add_option(
43 '--width',
44 dest='width',
45 type='float',
46 default=20.0,
47 metavar='FLOAT',
48 help='set width of output image [cm] (%default)')
50 parser.add_option(
51 '--height',
52 dest='height',
53 type='float',
54 default=15.0,
55 metavar='FLOAT',
56 help='set height of output image [cm] (%default)')
58 parser.add_option(
59 '--topo-resolution-min',
60 dest='topo_resolution_min',
61 type='float',
62 default=40.0,
63 metavar='FLOAT',
64 help='minimum resolution of topography [dpi] (%default)')
66 parser.add_option(
67 '--topo-resolution-max',
68 dest='topo_resolution_max',
69 type='float',
70 default=200.0,
71 metavar='FLOAT',
72 help='maximum resolution of topography [dpi] (%default)')
74 parser.add_option(
75 '--no-grid',
76 dest='show_grid',
77 default=True,
78 action='store_false',
79 help="don't show grid lines")
81 parser.add_option(
82 '--no-topo',
83 dest='show_topo',
84 default=True,
85 action='store_false',
86 help="don't show topography")
88 parser.add_option(
89 '--no-cities',
90 dest='show_cities',
91 default=True,
92 action='store_false',
93 help="don't show cities")
95 parser.add_option(
96 '--no-illuminate',
97 dest='illuminate',
98 default=True,
99 action='store_false',
100 help='deactivate artificial illumination of topography')
102 parser.add_option(
103 '--illuminate-factor-land',
104 dest='illuminate_factor_land',
105 type='float',
106 metavar='FLOAT',
107 help='set factor for artificial illumination of land (0.5)')
109 parser.add_option(
110 '--illuminate-factor-ocean',
111 dest='illuminate_factor_ocean',
112 type='float',
113 metavar='FLOAT',
114 help='set factor for artificial illumination of ocean (0.25)')
116 parser.add_option(
117 '--theme',
118 choices=['topo', 'soft'],
119 default='topo',
120 help='select color theme, available: topo, soft (topo)"')
122 parser.add_option(
123 '--download-etopo1',
124 dest='download_etopo1',
125 action='store_true',
126 help='download full ETOPO1 topography dataset')
128 parser.add_option(
129 '--download-srtmgl3',
130 dest='download_srtmgl3',
131 action='store_true',
132 help='download full SRTMGL3 topography dataset')
134 parser.add_option(
135 '--make-decimated-topo',
136 dest='make_decimated',
137 action='store_true',
138 help='pre-make all decimated topography datasets')
140 parser.add_option(
141 '--stations',
142 dest='stations_fn',
143 metavar='FILENAME',
144 help='load station coordinates from FILENAME')
146 parser.add_option(
147 '--events',
148 dest='events_fn',
149 metavar='FILENAME',
150 help='load event coordinates from FILENAME')
152 parser.add_option(
153 '--debug',
154 dest='debug',
155 action='store_true',
156 default=False,
157 help='print debugging information to stderr')
159 (options, args) = parser.parse_args(args)
161 if options.debug:
162 util.setup_logging(program_name, 'debug')
163 else:
164 util.setup_logging(program_name, 'info')
166 if options.download_etopo1:
167 import pyrocko.datasets.topo.etopo1
168 pyrocko.datasets.topo.etopo1.download()
170 if options.download_srtmgl3:
171 import pyrocko.datasets.topo.srtmgl3
172 pyrocko.datasets.topo.srtmgl3.download()
174 if options.make_decimated:
175 import pyrocko.datasets.topo
176 pyrocko.datasets.topo.make_all_missing_decimated()
178 if (options.download_etopo1 or options.download_srtmgl3 or
179 options.make_decimated) and len(args) == 0:
181 sys.exit(0)
183 if options.theme == 'soft':
184 color_kwargs = {
185 'illuminate_factor_land': options.illuminate_factor_land or 0.2,
186 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.15,
187 'color_wet': (216, 242, 254),
188 'color_dry': (238, 236, 230),
189 'topo_cpt_wet': 'light_sea_uniform',
190 'topo_cpt_dry': 'light_land_uniform'}
191 elif options.theme == 'topo':
192 color_kwargs = {
193 'illuminate_factor_land': options.illuminate_factor_land or 0.5,
194 'illuminate_factor_ocean': options.illuminate_factor_ocean or 0.25}
196 events = []
197 if options.events_fn:
198 events = model.load_events(options.events_fn)
200 stations = []
202 if options.stations_fn:
203 stations = model.load_stations(options.stations_fn)
205 if not (len(args) == 4 or (
206 len(args) == 1 and (events or stations))):
208 parser.print_help()
209 sys.exit(1)
211 if len(args) == 4:
212 try:
213 lat = float(args[0])
214 lon = float(args[1])
215 radius = float(args[2])*km
216 except Exception:
217 parser.print_help()
218 sys.exit(1)
219 else:
220 lats, lons = latlon_arrays(stations+events)
221 lat, lon = map(float, od.geographic_midpoint(lats, lons))
222 radius = float(
223 num.max(od.distance_accurate50m_numpy(lat, lon, lats, lons)))
224 radius *= 1.1
226 m = automap.Map(
227 width=options.width,
228 height=options.height,
229 lat=lat,
230 lon=lon,
231 radius=radius,
232 topo_resolution_max=options.topo_resolution_max,
233 topo_resolution_min=options.topo_resolution_min,
234 show_topo=options.show_topo,
235 show_grid=options.show_grid,
236 illuminate=options.illuminate,
237 **color_kwargs)
239 logger.debug('map configuration:\n%s' % str(m))
241 if options.show_cities:
242 m.draw_cities()
244 if stations:
245 lats = [s.lat for s in stations]
246 lons = [s.lon for s in stations]
248 m.gmt.psxy(
249 in_columns=(lons, lats),
250 S='t8p',
251 G='black',
252 *m.jxyr)
254 for s in stations:
255 m.add_label(s.lat, s.lon, '%s' % '.'.join(
256 x for x in s.nsl() if x))
258 if events:
259 beachball_symbol = 'mt'
260 beachball_size = 20.0
261 for ev in events:
262 if ev.moment_tensor is None:
263 m.gmt.psxy(
264 in_rows=[[ev.lon, ev.lat]],
265 S='c12p',
266 G=gmtpy.color('scarletred2'),
267 W='1p,black',
268 *m.jxyr)
270 else:
271 devi = ev.moment_tensor.deviatoric()
272 mt = devi.m_up_south_east()
273 mt = mt / ev.moment_tensor.scalar_moment() \
274 * pmt.magnitude_to_moment(5.0)
275 m6 = pmt.to6(mt)
276 data = (ev.lon, ev.lat, 10) + tuple(m6) + (1, 0, 0)
278 if m.gmt.is_gmt5():
279 kwargs = dict(
280 M=True,
281 S='%s%g' % (
282 beachball_symbol[0], beachball_size / gmtpy.cm))
283 else:
284 kwargs = dict(
285 S='%s%g' % (
286 beachball_symbol[0], beachball_size*2 / gmtpy.cm))
288 m.gmt.psmeca(
289 in_rows=[data],
290 G=gmtpy.color('chocolate1'),
291 E='white',
292 W='1p,%s' % gmtpy.color('chocolate3'),
293 *m.jxyr,
294 **kwargs)
296 m.save(args[-1])
299if __name__ == '__main__':
300 main()