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