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