1# https://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from __future__ import absolute_import, print_function, division
8# import copy
9import logging
10import operator
11import calendar
12import numpy as num
14from pyrocko.guts import \
15 Object, StringChoice, String, List, Float
17from pyrocko import table, model # , automap
18from pyrocko.client import fdsn, catalog
19from pyrocko.gui.qt_compat import qw
20from pyrocko.plot import beachball
21# from pyrocko.himesh import HiMesh
23# from pyrocko.gui.vtk_util import TrimeshPipe
25from .. import common
27from .table import TableElement, TableState
29logger = logging.getLogger('pyrocko.gui.sparrow.elements.catalog')
31guts_prefix = 'sparrow'
34attribute_names = [
35 'time', 'lat', 'lon', 'northing', 'easting', 'depth', 'magnitude']
37attribute_dtypes = [
38 'f16', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8']
40name_to_icol = dict(
41 (name, icol) for (icol, name) in enumerate(attribute_names))
43event_dtype = num.dtype(list(zip(attribute_names, attribute_dtypes)))
45t_time = num.float64
48def binned_statistic(values, ibins, function):
49 order = num.argsort(ibins)
50 values_sorted = values[order]
51 ibins_sorted = ibins[order]
52 parts = num.concatenate((
53 [0],
54 num.where(num.diff(ibins_sorted) != 0)[0] + 1,
55 [ibins.size]))
57 results = []
58 ibins_result = []
59 for ilow, ihigh in zip(parts[:-1], parts[1:]):
60 values_part = values_sorted[ilow:ihigh]
61 results.append(function(values_part))
62 ibins_result.append(ibins_sorted[ilow])
64 return num.array(ibins_result, dtype=num.int64), num.array(results)
67def load_text(
68 filepath,
69 column_names=('time', 'lat', 'lon', 'depth', 'magnitude'),
70 time_format='seconds'):
72 with open(filepath, 'r') as f:
73 if column_names == 'from_header':
74 line = f.readline()
75 column_names = line.split()
77 name_to_icol_in = dict(
78 (name, icol) for (icol, name) in enumerate(column_names)
79 if name in attribute_names)
81 data_in = num.loadtxt(filepath, skiprows=1)
83 nevents = data_in.shape[0]
84 c5 = num.zeros((nevents, 5))
85 c5[:, 0] = data_in[:, name_to_icol_in['lat']]
86 c5[:, 1] = data_in[:, name_to_icol_in['lon']]
87 c5[:, 2] = 0.0
88 c5[:, 3] = 0.0
89 c5[:, 4] = data_in[:, name_to_icol_in['depth']] * 1000.
91 tab = table.Table()
92 loc_rec = table.LocationRecipe()
93 tab.add_recipe(loc_rec)
95 tab.add_col(loc_rec.c5_header, c5)
96 for k, unit in [
97 ('time', 's'),
98 ('magnitude', None)]:
100 values = data_in[:, name_to_icol_in[k]]
102 if k == 'time' and time_format == 'year':
103 values = decimal_year_to_time(values)
105 tab.add_col(table.Header(k, unit), values)
107 return tab
110def decimal_year_to_time(year):
111 iyear_start = num.floor(year).astype(num.int64)
112 iyear_end = iyear_start + 1
114 iyear_min = num.min(iyear_start)
115 iyear_max = num.max(iyear_end)
117 iyear_to_time = num.zeros(iyear_max - iyear_min + 1, dtype=t_time)
118 for iyear in range(iyear_min, iyear_max+1):
119 iyear_to_time[iyear-iyear_min] = calendar.timegm(
120 (iyear, 1, 1, 0, 0, 0))
122 tyear_start = iyear_to_time[iyear_start - iyear_min]
123 tyear_end = iyear_to_time[iyear_end - iyear_min]
125 t = tyear_start + (year - iyear_start) * (tyear_end - tyear_start)
127 return t
130def oa_to_array(objects, attribute):
131 return num.fromiter(
132 map(operator.attrgetter(attribute), objects),
133 num.float64,
134 len(objects))
137def eventtags_to_array(events, tab):
139 ks = set()
140 event_tags = []
141 for ev in events:
142 tags = ev.tags_as_dict()
143 ks.update(tags.keys())
144 event_tags.append(tags)
146 for k in sorted(ks):
147 column = [tags.get(k, None) for tags in event_tags]
148 if all(isinstance(v, int) for v in column):
149 dtype = int
150 elif all(isinstance(v, float) for v in column):
151 dtype = float
152 else:
153 dtype = num.string_
154 column = [v or '' for v in column]
156 arr = num.array(column, dtype=dtype)
158 tab.add_col(table.Header('tag_%s' % k), arr)
161def eventextras_to_array(events, tab):
162 n_extras = num.array([len(ev.extras) for ev in events])
163 evts_all_extras = num.arange(len(events))[n_extras != 0]
165 if evts_all_extras.shape[0] == 0:
166 return tab
168 ev0 = events[evts_all_extras[0]]
170 if num.unique(n_extras).shape[0] > 1:
171 msg = 'Not all events have equal number of extras.'
172 if num.unique(n_extras).shape[0] == 2 and num.min(n_extras) == 0:
173 logger.warn(msg + ' Zero length lists are filled with NaNs.')
174 else:
175 raise IndexError(
176 msg + ' Several non-zero shapes detected. Please check.')
178 for key, val in ev0.extras.items():
179 dtype = num.string_
180 values = num.array(['' for x in range(n_extras.shape[0])], dtype=dtype)
182 if type(val) is float:
183 dtype, values = num.float64, num.ones_like(n_extras) * num.nan
185 elif type(val) is int:
186 dtype = num.int64
187 values = num.ones_like(n_extras, dtype=num.int64) * num.nan
189 header = 'extra_%s' % key
190 values[evts_all_extras] = num.array([
191 events[iev].extras[key] for iev in evts_all_extras], dtype=dtype)
193 tab.add_col(table.Header(header), values)
196def events_to_table(events):
197 c5 = num.zeros((len(events), 5))
198 m6 = None
200 for i, ev in enumerate(events):
201 c5[i, :] = (
202 ev.lat, ev.lon, ev.north_shift, ev.east_shift, ev.depth)
204 if ev.moment_tensor:
205 if m6 is None:
206 m6 = num.zeros((len(events), 6))
207 m6[:, 0:3] = 1.0
209 m6[i, :] = beachball.deco_part(ev.moment_tensor, 'deviatoric').m6()
211 tab = table.Table()
213 ev_rec = table.EventRecipe()
214 tab.add_recipe(ev_rec)
215 tab.add_col(ev_rec.get_header('c5'), c5)
217 for k in ['time', 'magnitude']:
218 tab.add_col(ev_rec.get_header(k), oa_to_array(events, k))
220 if events:
221 eventtags_to_array(events, tab)
222 eventextras_to_array(events, tab)
224 if m6 is not None:
225 mt_rec = table.MomentTensorRecipe()
226 tab.add_recipe(mt_rec)
227 tab.add_col(mt_rec.m6_header, m6)
229 return tab
232class LoadingChoice(StringChoice):
233 choices = [choice.upper() for choice in [
234 'file',
235 'fdsn']]
238class FDSNSiteChoice(StringChoice):
239 choices = [key.upper() for key in fdsn.g_site_abbr.keys()]
242class CatalogSelection(Object):
243 pass
246class MemoryCatalogSelection(CatalogSelection):
248 def __init__(self, events=None):
249 if events is None:
250 events = []
252 self._events = events
254 def get_table(self):
255 return events_to_table(self._events)
258class FileCatalogSelection(CatalogSelection):
259 paths = List.T(String.T())
261 def get_table(self):
262 from pyrocko.io import quakeml
264 events = []
265 for path in self.paths:
266 fn_ext = path.split('.')[-1].lower()
267 if fn_ext in ['xml', 'qml', 'quakeml']:
268 qml = quakeml.QuakeML.load_xml(filename=path)
269 events.extend(qml.get_pyrocko_events())
271 if fn_ext in ['dat', 'csv']:
272 assert len(self.paths) == 1
273 tab = load_text(
274 path, column_names='from_header', time_format='year')
276 return tab
278 else:
279 events.extend(model.load_events(path))
281 return events_to_table(events)
284g_catalogs = {
285 'Geofon': catalog.Geofon(),
286 'USGS/NEIC US': catalog.USGS('us'),
287 'Global-CMT': catalog.GlobalCMT(),
288 'Saxony (Uni-Leipzig)': catalog.Saxony()
289}
291g_fdsn_has_events = ['ISC', 'SCEDC', 'NCEDC', 'IRIS', 'GEONET']
293g_sites = sorted(g_catalogs.keys())
294g_sites.extend(g_fdsn_has_events)
297class OnlineCatalogSelection(CatalogSelection):
298 site = String.T()
299 tmin = Float.T()
300 tmax = Float.T()
301 magnitude_min = Float.T(optional=True)
303 def get_table(self):
304 logger.info('Getting events from "%s" catalog.' % self.site)
306 cat = g_catalogs.get(self.site, None)
307 try:
308 if cat:
309 kwargs = {}
310 if self.magnitude_min is not None:
311 kwargs['magmin'] = self.magnitude_min
313 events = cat.get_events(
314 time_range=(self.tmin, self.tmax), **kwargs)
316 else:
317 kwargs = {}
318 if self.magnitude_min is not None:
319 kwargs['minmagnitude'] = self.magnitude_min
321 request = fdsn.event(
322 starttime=self.tmin,
323 endtime=self.tmax,
324 site=self.site, **kwargs)
326 from pyrocko.io import quakeml
327 qml = quakeml.QuakeML.load_xml(request)
328 events = qml.get_pyrocko_events()
330 logger.info('Got %i event%s from "%s" catalog.' % (
331 len(events), '' if len(events) == 1 else 's', self.site))
333 except Exception as e:
334 logger.error(
335 'Getting events from "%s" catalog failed: %s' % (
336 self.site, str(e)))
338 events = []
340 return events_to_table(events)
343class CatalogState(TableState):
344 selection = CatalogSelection.T(optional=True)
345 size_parameter = String.T(default='magnitude', optional=True)
346 color_parameter = String.T(default='depth', optional=True)
348 @classmethod
349 def get_name(self):
350 return 'Catalog'
352 def create(self):
353 element = CatalogElement()
354 return element
357class CatalogElement(TableElement):
359 def __init__(self, *args, **kwargs):
360 TableElement.__init__(self, *args, **kwargs)
361 self._selection_view = None
362 # self._himesh = HiMesh(order=6)
364 # cpt_data = [
365 # (0.0, 0.0, 0.0, 0.0),
366 # (1.0, 0.9, 0.9, 0.2)]
367 #
368 # self.cpt_mesh = automap.CPT(
369 # levels=[
370 # automap.CPTLevel(
371 # vmin=a[0],
372 # vmax=b[0],
373 # color_min=[255*x for x in a[1:]],
374 # color_max=[255*x for x in b[1:]])
375 # for (a, b) in zip(cpt_data[:-1], cpt_data[1:])])
377 def get_name(self):
378 return 'Catalog'
380 def get_size_parameter_extra_entries(self):
381 return ['magnitude']
383 def get_color_parameter_extra_entries(self):
384 return ['depth']
386 def bind_state(self, state):
387 TableElement.bind_state(self, state)
388 self.register_state_listener3(self.update, state, 'selection')
390 def update(self, *args):
391 state = self._state
392 # ifaces = None
393 if self._selection_view is not state.selection:
394 self.set_table(state.selection.get_table())
395 self._selection_view = state.selection
397 # ifaces = self._himesh.points_to_faces(self._table.get_col('xyz'))
399 TableElement.update(self, *args)
401 # if ifaces is not None:
402 # ifaces_x, sizes = binned_statistic(
403 # ifaces, ifaces, lambda part: part.shape[0])
404 #
405 # vertices = self._himesh.get_vertices()
406 # # vertices *= 0.95
407 # faces = self._himesh.get_faces()
408 #
409 # values = num.zeros(faces.shape[0])
410 # values[ifaces_x] = num.log(1+sizes)
411 #
412 # self._mesh = TrimeshPipe(vertices, faces, values=values)
413 # cpt = copy.deepcopy(self.cpt_mesh)
414 # cpt.scale(num.min(values), num.max(values))
415 # self._mesh.set_cpt(cpt)
416 # self._mesh.set_opacity(0.5)
418 # self._parent.add_actor(self._mesh.actor)
420 def open_file_dialog(self):
421 caption = 'Select one or more files to open'
423 fns, _ = qw.QFileDialog.getOpenFileNames(
424 self._parent, caption, options=common.qfiledialog_options)
426 self._state.selection = FileCatalogSelection(
427 paths=[str(fn) for fn in fns])
429 def open_catalog_load_dialog(self):
430 dialog = qw.QDialog(self._parent)
431 dialog.setWindowTitle('Get events from online catalog.')
433 layout = qw.QHBoxLayout(dialog)
435 layout.addWidget(qw.QLabel('Site'))
437 cb = qw.QComboBox()
438 for i, s in enumerate(g_sites):
439 cb.insertItem(i, s)
441 layout.addWidget(cb)
443 pb = qw.QPushButton('Cancel')
444 pb.clicked.connect(dialog.reject)
445 layout.addWidget(pb)
447 pb = qw.QPushButton('Ok')
448 pb.clicked.connect(dialog.accept)
449 layout.addWidget(pb)
451 dialog.exec_()
453 site = str(cb.currentText())
455 vstate = self._parent.state
457 if dialog.result() == qw.QDialog.Accepted:
458 self._state.selection = OnlineCatalogSelection(
459 site=site,
460 tmin=vstate.tmin,
461 tmax=vstate.tmax)
463 def _get_table_widgets_start(self):
464 return 1 # used as y arg in addWidget calls
466 def _get_controls(self):
467 if not self._controls:
468 frame = TableElement._get_controls(self) # sets self._controls
469 layout = frame.layout()
471 lab = qw.QLabel('Load from')
472 pb_file = qw.QPushButton('File')
474 layout.addWidget(lab, 0, 0)
475 layout.addWidget(pb_file, 0, 1)
477 pb_file.clicked.connect(self.open_file_dialog)
479 pb_file = qw.QPushButton('Online Catalog')
481 layout.addWidget(lab, 0, 0)
482 layout.addWidget(pb_file, 0, 2)
484 pb_file.clicked.connect(self.open_catalog_load_dialog)
486 return self._controls
489__all__ = [
490 'CatalogSelection',
491 'FileCatalogSelection',
492 'MemoryCatalogSelection',
493 'CatalogElement',
494 'CatalogState',
495]