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