1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from __future__ import absolute_import, print_function, division
8from pyrocko.gui.qt_compat import qc, qg
9import os
10from os import path as op
11import sys
13import numpy as num
14import logging
15import posixpath
16import shutil
17import socket
19if sys.version_info >= (3, 0): # noqa
20 from http.server import HTTPServer, SimpleHTTPRequestHandler
21else: # noqa
22 from BaseHTTPServer import HTTPServer
23 from SimpleHTTPServer import SimpleHTTPRequestHandler
25from pyrocko.gui.snuffler.snuffling import Snuffling, Switch, Choice, \
26 NoViewerSet
27from pyrocko.guts import dump_xml
28from pyrocko import util, model, orthodrome as ortho
29from pyrocko.gui import util as gui_util
30from pyrocko import moment_tensor
31from pyrocko.plot.automap import Map
32from pyrocko.util import unquote
33from .xmlMarker import XMLEventMarker, EventMarkerList, XMLStationMarker
34from .xmlMarker import StationMarkerList, MarkerLists
37g_counter = 0
39logger = logging.getLogger('pyrocko.snuffling.map')
42def port_in_use(port):
43 with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
44 return s.connect_ex(('localhost', port)) == 0
47def get_port():
48 port = 9998
49 while port < 11000:
50 if port_in_use(port):
51 port += 1
52 else:
53 break
54 else:
55 raise Exception('No free port found')
57 return port
60def get_magnitude(event):
61 if event.magnitude:
62 mag = event.magnitude
63 elif event.moment_tensor:
64 mag = event.moment_tensor.moment_magnitude()
65 else:
66 mag = 0.
67 return float(mag)
70def convert_event_marker(marker):
71 ev = marker.get_event()
72 depth = None
73 if ev is not None:
74 depth = ev.depth
75 else:
76 return None
78 if depth is None:
79 depth = 0.0
80 ev_name = ev.name if ev.name else '(Event)'
81 xmleventmarker = XMLEventMarker(
82 eventname=ev_name,
83 longitude=float(ev.effective_lon),
84 latitude=float(ev.effective_lat),
85 origintime=util.time_to_str(ev.time),
86 depth=float(depth),
87 magnitude=float(get_magnitude(ev)),
88 active=['no', 'yes'][marker.active])
90 return xmleventmarker
93class RootedHTTPServer(HTTPServer):
95 def __init__(self, base_path, *args, **kwargs):
96 HTTPServer.__init__(self, *args, **kwargs)
97 self.RequestHandlerClass.base_path = base_path
100class RootedHTTPRequestHandler(SimpleHTTPRequestHandler):
102 def log_message(self, format, *args):
103 logger.debug("%s - - [%s] %s\n" % (
104 self.client_address[0],
105 self.log_date_time_string(),
106 format % args))
108 def translate_path(self, path):
109 path = posixpath.normpath(unquote(path))
110 words = path.split('/')
111 words = filter(None, words)
112 path = self.base_path
113 for word in words:
114 drive, word = os.path.splitdrive(word)
115 head, word = os.path.split(word)
116 if word in (os.curdir, os.pardir):
117 continue
118 path = os.path.join(path, word)
119 return path
122class FileServingWorker(qc.QObject):
124 def __init__(self, dirname, *args, **kwargs):
125 super(FileServingWorker, self).__init__(*args, **kwargs)
126 self.dirname = dirname
127 self.httpd = None
129 @qc.pyqtSlot(int)
130 def run(self, port):
131 server_address = ('', port)
132 self.httpd = RootedHTTPServer(
133 self.dirname, server_address, RootedHTTPRequestHandler)
134 sa = self.httpd.socket.getsockname()
135 logger.debug('Serving on port %s' % sa[1])
136 self.httpd.serve_forever()
138 @qc.pyqtSlot()
139 def stop(self):
140 if self.httpd:
141 logger.debug('Shutdown server')
142 self.httpd.shutdown()
143 self.httpd.socket.close()
146class SignalProxy(qc.QObject):
147 content_to_serve = qc.pyqtSignal(int)
150class MapMaker(Snuffling):
151 '''
152 <html>
153 <body>
154 <h1>Map event and stations with OpenStreetMap or Google Maps</h1>
155 <p>
156 Invokes the standard browser if "Open in external browser" is selected.
157 Some browsers do not allow javascript to open and read the xml-file
158 containing the necessary information due to the "Same-Origin-Policy".
159 In that case you need to reset your standard browser. I.e.: Firefox on
160 Linux do: <tt>xdg-settings set default-web-browser firefox.desktop</tt>
161 </p>
162 <p>
163 Clicking one of the plate boundary lines shows a reference regarding that
164 plate boundary.
165 </p>
166 <p>
167 The plate boundary database is based on the work done by Peter Bird, who
168 kindly permitted usage. <br>
169 <p>
170 This snuffling can also be called from the command line, if it is stored in
171 the default pyrocko location under $HOME/.snufflings<br>
172 e.g.:
173 <p>
174 <code>
175python $HOME/.snufflings/map/snuffling.py --stations=stations.pf
176--events=events_test.pf
177</code>
178 <h2>References</h2>
179 <i>50. Bird, P. (2003) An updated digital model of plate
180 boundaries, Geochemistry Geophysics Geosystems, 4(3), 1027,
181 doi:10.1029/2001GC000252.
182 </i>
183 </p>
184 <br>
185 Also available at
186 <a href="http://peterbird.name/publications/2003_PB2002/2003_PB2002.htm">
187 http://www.peterbird.name</a>
188 <p>
189 Please note, that in the current implementation the orogens (cross-hatched
190 areas in
191 <a href="http://peterbird.name/publications/2003_PB2002/Figure_01.gif">
192 figure 1</a>)
193 are not distinguished from plate boundaries. The orogens are meant to
194 mark areas where the plate model is known to be incomplete
195 (and/or inapplicable).<br>
196 This matter will be pointed out in future releases of this snuffling.
197 </body>
198 </html>
199 '''
201 def setup(self):
202 self.set_name('Map')
203 self.add_parameter(Switch('Only active event', 'only_active', False))
204 self.add_parameter(Switch('Open in external browser',
205 'open_external', False))
206 self.add_parameter(Choice('Provider', 'map_kind', 'OpenStreetMap',
207 ['OpenStreetMap', 'GMT', 'Google Maps']))
209 self.set_live_update(False)
210 self.figcount = 0
211 self.thread = qc.QThread()
212 self.marker_tempdir = self.tempdir()
214 self.file_serving_worker = FileServingWorker(self.marker_tempdir)
215 self.file_serving_worker.moveToThread(self.thread)
216 self.data_proxy = SignalProxy()
217 self.data_proxy.content_to_serve.connect(
218 self.file_serving_worker.run)
219 self.thread.start()
220 self.port = None
222 self.viewer_connected = False
224 def call(self):
225 if self.port is None:
226 self.port = get_port()
228 if not self.viewer_connected:
229 self.get_viewer().about_to_close.connect(
230 self.file_serving_worker.stop)
231 self.viewer_connected = True
233 self.cleanup()
235 try:
236 viewer = self.get_viewer()
237 cli_mode = False
238 except NoViewerSet:
239 viewer = None
240 cli_mode = True
242 if not cli_mode:
243 if self.only_active:
244 _, active_stations = \
245 self.get_active_event_and_stations()
246 else:
247 active_stations = viewer.stations.values()
248 elif cli_mode:
249 active_stations = self.stations
251 station_list = []
252 if active_stations:
253 for stat in active_stations:
254 is_blacklisted = util.match_nslc(viewer.blacklist, stat.nsl())
255 if (viewer and not is_blacklisted) or cli_mode:
256 xml_station_marker = XMLStationMarker(
257 nsl='.'.join(stat.nsl()),
258 longitude=float(stat.effective_lon),
259 latitude=float(stat.effective_lat),
260 active='yes')
262 station_list.append(xml_station_marker)
264 active_station_list = StationMarkerList(stations=station_list)
266 if self.only_active:
267 markers = [viewer.get_active_event_marker()]
268 else:
269 if cli_mode:
270 markers = self.markers
271 else:
272 markers = self.get_selected_markers()
273 if len(markers) == 0:
274 tmin, tmax = self.get_selected_time_range(fallback=True)
275 markers = [m for m in viewer.get_markers()
276 if isinstance(m, gui_util.EventMarker) and
277 m.tmin >= tmin and m.tmax <= tmax]
279 ev_marker_list = []
280 for m in markers:
281 if not isinstance(m, gui_util.EventMarker):
282 continue
283 xmleventmarker = convert_event_marker(m)
284 if xmleventmarker is None:
285 continue
286 ev_marker_list.append(xmleventmarker)
288 event_list = EventMarkerList(events=ev_marker_list)
289 event_station_list = MarkerLists(
290 station_marker_list=active_station_list,
291 event_marker_list=event_list)
293 event_station_list.validate()
294 if self.map_kind != 'GMT':
295 tempdir = self.marker_tempdir
296 if self.map_kind == 'Google Maps':
297 map_fn = 'map_googlemaps.html'
298 elif self.map_kind == 'OpenStreetMap':
299 map_fn = 'map_osm.html'
301 url = 'http://localhost:' + str(self.port) + '/%s' % map_fn
303 files = ['loadxmldoc.js', 'map_util.js', 'plates.kml', map_fn]
304 snuffling_dir = op.dirname(op.abspath(__file__))
305 for entry in files:
306 shutil.copy(os.path.join(snuffling_dir, entry),
307 os.path.join(tempdir, entry))
308 logger.debug('copied data to %s' % tempdir)
309 markers_fn = os.path.join(self.marker_tempdir, 'markers.xml')
310 self.data_proxy.content_to_serve.emit(self.port)
311 dump_xml(event_station_list, filename=markers_fn)
313 if self.open_external:
314 qg.QDesktopServices.openUrl(qc.QUrl(url))
315 else:
316 global g_counter
317 g_counter += 1
318 self.web_frame(
319 url,
320 name='Map %i (%s)' % (g_counter, self.map_kind))
321 else:
322 lats_all = []
323 lons_all = []
325 slats = []
326 slons = []
327 slabels = []
328 for s in active_stations:
329 slats.append(s.effective_lat)
330 slons.append(s.effective_lon)
331 slabels.append('.'.join(s.nsl()))
333 elats = []
334 elons = []
335 elats = []
336 elons = []
337 psmeca_input = []
338 markers = self.get_selected_markers()
339 for m in markers:
340 if isinstance(m, gui_util.EventMarker):
341 e = m.get_event()
342 elats.append(e.effective_lat)
343 elons.append(e.effective_lon)
344 if e.moment_tensor is not None:
345 mt = e.moment_tensor.m6()
346 psmeca_input.append(
347 (e.effective_lon, e.effective_lat,
348 e.depth/1000., mt[0], mt[1],
349 mt[2], mt[3], mt[4], mt[5],
350 1., e.effective_lon, e.effective_lat, e.name))
351 else:
352 if e.magnitude is None:
353 moment = -1.
354 else:
355 moment = moment_tensor.magnitude_to_moment(
356 e.magnitude)
357 psmeca_input.append(
358 (e.effective_lon, e.effective_lat,
359 e.depth/1000.,
360 moment/3., moment/3., moment/3.,
361 0., 0., 0., 1.,
362 e.effective_lon, e.effective_lat, e.name))
364 lats_all.extend(elats)
365 lons_all.extend(elons)
366 lats_all.extend(slats)
367 lons_all.extend(slons)
369 lats_all = num.array(lats_all)
370 lons_all = num.array(lons_all)
372 if len(lats_all) == 0:
373 return
375 center_lat, center_lon = ortho.geographic_midpoint(
376 lats_all, lons_all)
377 ntotal = len(lats_all)
378 clats = num.ones(ntotal) * center_lat
379 clons = num.ones(ntotal) * center_lon
380 dists = ortho.distance_accurate50m_numpy(
381 clats, clons, lats_all, lons_all)
383 maxd = num.max(dists) or 0.
384 m = Map(
385 lat=center_lat, lon=center_lon,
386 radius=max(10000., maxd) * 1.1,
387 width=35, height=25,
388 show_grid=True,
389 show_topo=True,
390 color_dry=(238, 236, 230),
391 topo_cpt_wet='light_sea_uniform',
392 topo_cpt_dry='light_land_uniform',
393 illuminate=True,
394 illuminate_factor_ocean=0.15,
395 show_rivers=False,
396 show_plates=False)
398 m.gmt.psxy(in_columns=(slons, slats), S='t15p', G='black', *m.jxyr)
399 for i in range(len(active_stations)):
400 m.add_label(slats[i], slons[i], slabels[i])
402 m.gmt.psmeca(
403 in_rows=psmeca_input, S='m1.0', G='red', C='5p,0/0/0', *m.jxyr)
405 tmpdir = self.tempdir()
407 self.outfn = os.path.join(tmpdir, '%i.png' % self.figcount)
408 m.save(self.outfn)
409 f = self.pixmap_frame(self.outfn) # noqa
410 # f = self.svg_frame(self.outfn, name='XX')
412 def configure_cli_parser(self, parser):
414 parser.add_option(
415 '--events',
416 dest='events_filename',
417 default=None,
418 metavar='FILENAME',
419 help='Read markers from FILENAME')
421 parser.add_option(
422 '--markers',
423 dest='markers_filename',
424 default=None,
425 metavar='FILENAME',
426 help='Read markers from FILENAME')
428 parser.add_option(
429 '--stations',
430 dest='stations_filename',
431 default=None,
432 metavar='FILENAME',
433 help='Read stations from FILENAME')
435 parser.add_option(
436 '--provider',
437 dest='map_provider',
438 default='google',
439 help='map provider [google | osm] (default=osm)')
441 def __del__(self):
442 self.file_serving_worker.stop()
443 self.thread.quit()
444 self.thread.wait()
447def __snufflings__():
448 return [MapMaker()]
451if __name__ == '__main__':
452 util.setup_logging('map.py', 'info')
453 s = MapMaker()
454 options, args, parser = s.setup_cli()
455 s.markers = []
457 if options.stations_filename:
458 stations = model.load_stations(options.stations_filename)
459 s.stations = stations
460 else:
461 s.stations = None
463 if options.events_filename:
464 events = model.load_events(filename=options.events_filename)
465 markers = [gui_util.EventMarker(e) for e in events]
466 s.markers.extend(markers)
468 if options.markers_filename:
469 markers = gui_util.load_markers(options.markers_filename)
470 s.markers.extend(markers)
471 s.open_external = True
472 mapmap = {'google': 'Google Maps', 'osm': 'OpenStreetMap'}
473 s.map_kind = mapmap[options.map_provider]
474 s.call()