Coverage for /usr/local/lib/python3.11/dist-packages/pyrocko/gui/snuffler/snufflings/map/__init__.py: 29%

236 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-10-06 06:59 +0000

1# https://pyrocko.org - GPLv3 

2# 

3# The Pyrocko Developers, 21st Century 

4# ---|P------/S----------~Lg---------- 

5 

6from pyrocko.gui.qt_compat import qc, qg 

7import os 

8from os import path as op 

9import sys 

10 

11import numpy as num 

12import logging 

13import posixpath 

14import shutil 

15import socket 

16 

17if sys.version_info >= (3, 0): # noqa 

18 from http.server import HTTPServer, SimpleHTTPRequestHandler 

19else: # noqa 

20 from BaseHTTPServer import HTTPServer 

21 from SimpleHTTPServer import SimpleHTTPRequestHandler 

22 

23from pyrocko.gui.snuffler.snuffling import Snuffling, Switch, Choice, \ 

24 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 

33 

34 

35g_counter = 0 

36 

37logger = logging.getLogger('pyrocko.snuffling.map') 

38 

39 

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 

43 

44 

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') 

54 

55 return port 

56 

57 

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) 

66 

67 

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 

75 

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]) 

87 

88 return xmleventmarker 

89 

90 

91class RootedHTTPServer(HTTPServer): 

92 

93 def __init__(self, base_path, *args, **kwargs): 

94 HTTPServer.__init__(self, *args, **kwargs) 

95 self.RequestHandlerClass.base_path = base_path 

96 

97 

98class RootedHTTPRequestHandler(SimpleHTTPRequestHandler): 

99 

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)) 

105 

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 

118 

119 

120class FileServingWorker(qc.QObject): 

121 

122 def __init__(self, dirname, *args, **kwargs): 

123 super(FileServingWorker, self).__init__(*args, **kwargs) 

124 self.dirname = dirname 

125 self.httpd = None 

126 

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() 

135 

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() 

142 

143 

144class SignalProxy(qc.QObject): 

145 content_to_serve = qc.pyqtSignal(int) 

146 

147 

148class MapMaker(Snuffling): 

149 

150 def help(self): 

151 return ''' 

152<html> 

153<body> 

154<h1>Map event and stations with OpenStreetMap or Google Maps</h1> 

155<p> 

156Invokes the standard browser if "Open in external browser" is selected. 

157Some browsers do not allow javascript to open and read the xml-file 

158containing the necessary information due to the "Same-Origin-Policy". 

159In that case you need to reset your standard browser. I.e.: Firefox on 

160Linux do: <tt>xdg-settings set default-web-browser firefox.desktop</tt> 

161</p> 

162<p> 

163Clicking one of the plate boundary lines shows a reference regarding that 

164plate boundary. 

165</p> 

166<p> 

167The plate boundary database is based on the work done by Peter Bird, who 

168kindly permitted usage. <br> 

169<p> 

170This snuffling can also be called from the command line, if it is stored in 

171the default pyrocko location under $HOME/.snufflings<br> 

172e.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 

180boundaries, Geochemistry Geophysics Geosystems, 4(3), 1027, 

181doi:10.1029/2001GC000252. 

182</i> 

183</p> 

184<br> 

185Also available at 

186<a href="http://peterbird.name/publications/2003_PB2002/2003_PB2002.htm"> 

187 http://www.peterbird.name</a> 

188<p> 

189Please note, that in the current implementation the orogens (cross-hatched 

190areas in 

191<a href="http://peterbird.name/publications/2003_PB2002/Figure_01.gif"> 

192figure 1</a>) 

193are not distinguished from plate boundaries. The orogens are meant to 

194mark areas where the plate model is known to be incomplete 

195(and/or inapplicable).<br> 

196This matter will be pointed out in future releases of this snuffling. 

197</body> 

198</html> 

199''' 

200 

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'])) 

208 

209 self.set_live_update(False) 

210 self.figcount = 0 

211 self.thread = qc.QThread() 

212 self.marker_tempdir = self.tempdir() 

213 

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 

221 

222 self.viewer_connected = False 

223 

224 def call(self): 

225 if self.port is None: 

226 self.port = get_port() 

227 

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 

232 

233 self.cleanup() 

234 

235 try: 

236 viewer = self.get_viewer() 

237 cli_mode = False 

238 except NoViewerSet: 

239 viewer = None 

240 cli_mode = True 

241 

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 

250 

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') 

261 

262 station_list.append(xml_station_marker) 

263 

264 active_station_list = StationMarkerList(stations=station_list) 

265 

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] 

278 

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) 

287 

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) 

292 

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' 

300 

301 url = 'http://localhost:' + str(self.port) + '/%s' % map_fn 

302 

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) 

312 

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 = [] 

324 

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())) 

332 

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)) 

363 

364 lats_all.extend(elats) 

365 lons_all.extend(elons) 

366 lats_all.extend(slats) 

367 lons_all.extend(slons) 

368 

369 lats_all = num.array(lats_all) 

370 lons_all = num.array(lons_all) 

371 

372 if len(lats_all) == 0: 

373 return 

374 

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) 

382 

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) 

397 

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]) 

401 

402 m.gmt.psmeca( 

403 in_rows=psmeca_input, S='m1.0', G='red', C='5p,0/0/0', *m.jxyr) 

404 

405 tmpdir = self.tempdir() 

406 

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') 

411 

412 def configure_cli_parser(self, parser): 

413 

414 parser.add_option( 

415 '--events', 

416 dest='events_filename', 

417 default=None, 

418 metavar='FILENAME', 

419 help='Read markers from FILENAME') 

420 

421 parser.add_option( 

422 '--markers', 

423 dest='markers_filename', 

424 default=None, 

425 metavar='FILENAME', 

426 help='Read markers from FILENAME') 

427 

428 parser.add_option( 

429 '--stations', 

430 dest='stations_filename', 

431 default=None, 

432 metavar='FILENAME', 

433 help='Read stations from FILENAME') 

434 

435 parser.add_option( 

436 '--provider', 

437 dest='map_provider', 

438 default='google', 

439 help='map provider [google | osm] (default=osm)') 

440 

441 def __del__(self): 

442 self.file_serving_worker.stop() 

443 self.thread.quit() 

444 self.thread.wait() 

445 

446 

447def __snufflings__(): 

448 return [MapMaker()] 

449 

450 

451if __name__ == '__main__': 

452 util.setup_logging('map.py', 'info') 

453 s = MapMaker() 

454 options, args, parser = s.setup_cli() 

455 s.markers = [] 

456 

457 if options.stations_filename: 

458 stations = model.load_stations(options.stations_filename) 

459 s.stations = stations 

460 else: 

461 s.stations = None 

462 

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) 

467 

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()