Commit bee0e4c2 authored by Kris Vanneste's avatar Kris Vanneste
Browse files

Empty streams in clear_all method of MainWindow.

Added load_from_asdf method to MainWindow.
Added 'download_meth' argument to event_chosen method of MainWindow.
More robust determination of event label in update_event_display method of MainWindow.
Small improvements in _on_download_rob_data method of MainWindow.
Added 'force_reprocessing' argument to _on_process_data method of MainWindow.
Updated _update_mean_parameter_display method of MainWindow to use minus/plus errors for fc, r, SD and Q.
Calculate RMSE if missing in _on_file_load_belshake method of MainWindow.
Added calc_rmse method to MainWindow.
parent 089d61ae
......@@ -462,6 +462,8 @@ class MainWindow(QtWidgets.QMainWindow):
self.ui.seishub_server.setText("")
self.update_hypocentral_velocities()
self.time_windows = {'P': {}, 'S': {}}
self.streams = {}
self.proc_streams = {}
self.update_menu()
def _get_options_submenu(self, submenu_name):
......@@ -554,7 +556,7 @@ class MainWindow(QtWidgets.QMainWindow):
action.setChecked(self.current_state['restitution_method'] == action_model)
if self.current_state["event"]:
self._on_process_data()
self._on_process_data(force_reprocessing=True)
def _on_phase_comp_sel_change(self, i):
"""
......@@ -619,9 +621,79 @@ class MainWindow(QtWidgets.QMainWindow):
window.raise_()
self.hide()
def load_from_asdf(self, asdf_file, tag):
"""
Load event and waveforms from ASDF file
"""
ds = robspy.io.asdf.open_asdf(asdf_file)
## Extract event and associated picks
ev = robspy.io.asdf.extract_event(ds, keep_open=True)
if not ev:
QtWidgets.QMessageBox.critical(self, "Error",
"Selected event not found")
return
if len(ev.picks) == 0:
QtWidgets.QMessageBox.critical(self, "Error",
"Selected event has no associated picks")
## Determine available network-station codes
nw_stat_codes = []
for pick in ev.picks:
nw_stat_code = pick.waveform_id.network_code + '.' + pick.waveform_id.station_code
if not nw_stat_code in nw_stat_codes:
nw_stat_codes.append(nw_stat_code)
## Initiate progress dialog
progress_dialog = QtWidgets.QProgressDialog(
"Downloading waveform data...", "Cancel", 0,
len(nw_stat_codes))
progress_dialog.setWindowModality(QtCore.Qt.WindowModal)
progress_dialog.forceShow()
## Extract waveform data
self.streams = {}
self.proc_streams = {}
for _i, nw_stat_code in enumerate(nw_stat_codes):
if progress_dialog.wasCanceled():
break
progress_dialog.setValue(_i)
sg = None
if robspy.io.asdf.has_tag(ds, nw_stat_code, tag, keep_open=True):
sg = robspy.io.asdf.extract_seismogram(ds, tag, nw_stat_code,
attach_response=True, keep_open=True)
if sg:
print(sg)
network, station = nw_stat_code.split('.')
self.streams[(network, station)] = sg
else:
print("Could not download waveform data for %s.%s"
% (network, station))
## Remove picks for stations without any data
pick_idxs = []
for p, pick in enumerate(ev.picks):
if (pick.waveform_id.network_code == network
and pick.waveform_id.station_code == station):
pick_idxs.append(p)
for idx in sorted(pick_idxs)[::-1]:
ev.picks.pop(idx)
# Finish the progress dialog.
ds._close()
progress_dialog.setValue(len(nw_stat_codes))
self.current_state['event'] = ev
self.event_chosen(ev, download_meth=None)
self._on_process_data()
def _load_rob_event(self):
"""
Read event and phase picks from ROB database
Read event and phase picks from ROB or Belshake database
"""
event_id = self.ui.seishub_server.text()
try:
......@@ -640,15 +712,15 @@ class MainWindow(QtWidgets.QMainWindow):
if not ev:
QtWidgets.QMessageBox.critical(self, "Error",
"Selected event not found")
"Selected event not found")
return
if len(ev.picks) == 0:
QtWidgets.QMessageBox.critical(self, "Error",
"Selected event has no associated picks")
"Selected event has no associated picks")
self.current_state['event'] = ev
self.event_chosen(ev)
self.event_chosen(ev, '_on_download_rob_data')
def fetch_event_from_db(self, event_id):
"""
......@@ -679,7 +751,6 @@ class MainWindow(QtWidgets.QMainWindow):
return ev
# TODO: store focal mechanism in database
def fetch_belshake_event(self, event_id):
"""
"""
......@@ -704,7 +775,7 @@ class MainWindow(QtWidgets.QMainWindow):
return ev
def event_chosen(self, event):
def event_chosen(self, event, download_meth='_on_download_data'):
"""
This is called when a new event has been selected.
......@@ -724,7 +795,8 @@ class MainWindow(QtWidgets.QMainWindow):
self.update_hypocentral_velocities()
# Download the data upon loading.
self._on_download_rob_data()
if download_meth:
getattr(self, download_meth)()
self.populate_result_table()
self.update_menu()
......@@ -737,7 +809,8 @@ class MainWindow(QtWidgets.QMainWindow):
event = self.current_state["event"]
if event:
event_label = "%s - %s"
event_label %= (event.resource_id, event.event_descriptions[0])
event_name = event.event_descriptions[0] if len(event.event_descriptions) else ''
event_label %= (event.resource_id, event_name)
lon_label = "%.4f" % event.origins[0].longitude
lat_label = "%.4f" % event.origins[0].latitude
depth_label = "%.4f" % (event.origins[0].depth / 1000.)
......@@ -843,10 +916,10 @@ class MainWindow(QtWidgets.QMainWindow):
"""
from waveforms import download_rob_seismogram
if "event" not in self.current_state:
event = self.current_state.get("event")
if not event:
return
event = self.current_state["event"]
event_ID = int(event.resource_id.id.split("/")[-1])
print(event)
......@@ -871,6 +944,7 @@ class MainWindow(QtWidgets.QMainWindow):
belshakedb = BelshakeDB()
self.streams = {}
self.proc_streams = {}
for _i, ((network, station), station_picks) in enumerate(sorted(picks.items())):
if progress_dialog.wasCanceled():
break
......@@ -922,13 +996,13 @@ class MainWindow(QtWidgets.QMainWindow):
self._on_process_data()
def _on_process_data(self):
def _on_process_data(self, force_reprocessing=False):
"""
(Re)process waveforms and store trimmed traces in Pick objects
of selected event
"""
import warnings
from waveforms import (process_waveforms, get_phase_pick)
from waveforms import (process_waveforms, post_process_waveforms, get_phase_pick)
## Initiate progress dialog
progress_dialog = QtWidgets.QProgressDialog(
......@@ -941,6 +1015,7 @@ class MainWindow(QtWidgets.QMainWindow):
phase_component_sel = self.current_state['phase_component_sel']
buffer_len = self.ui.buffer_seconds.value()
for _i, ((network, station), sg) in enumerate(sorted(self.streams.items())):
if progress_dialog.wasCanceled():
break
......@@ -952,8 +1027,12 @@ class MainWindow(QtWidgets.QMainWindow):
with warnings.catch_warnings():
warnings.simplefilter("ignore")
restitution_method = self.current_state['restitution_method']
process_waveforms(sg, event, station_picks, buffer_len,
phase_component_sel, restitution_method)
proc_sg = self.proc_streams.get((network, station))
if proc_sg is None or force_reprocessing:
proc_sg = process_waveforms(sg, event, restitution_method)
self.proc_streams[(network, station)] = proc_sg
post_process_waveforms(proc_sg, sg, event, station_picks, buffer_len,
phase_component_sel)
# Finish the progress dialog.
progress_dialog.setValue(len(self.streams))
......@@ -1246,6 +1325,7 @@ class MainWindow(QtWidgets.QMainWindow):
pick = self.current_state["pick"]
imin = int(((pick.time - trace.stats.starttime) + tmin) / trace.stats.delta)
imax = int(((pick.time - trace.stats.starttime) + tmax) / trace.stats.delta)
print(imin, imax)
freq, spec, jackknife_errors = calculate_spectrum(trace, imin, imax,
multi_tapered=multi_tapered)
......@@ -1505,9 +1585,10 @@ class MainWindow(QtWidgets.QMainWindow):
# Update the values.
self._update_spectral_parameter_display()
## Compute RMSE
print(omega_0, fc, Q, travel_time)
ax.rmse = calc_spectral_misfit(ax.spectrum, ax.frequencies,
theoretical_spectrum, freqs, ax.fmin, ax.fmax)
print('RMSE: %f' % ax.rmse)
print('RMSE: %f' % ax.rmse, ax.fmin, ax.fmax)
def _draw_fc(self):
"""
......@@ -1994,9 +2075,9 @@ class MainWindow(QtWidgets.QMainWindow):
and (np.isclose(final_result[wt]['MW'], self.final_result[wt]['MW'])
and np.isclose(final_result[wt]['MW_sigma'], self.final_result[wt]['MW_sigma'])
and np.isclose(final_result[wt]['r'], self.final_result[wt]['r'])
and np.isclose(final_result[wt]['r_sigma'], self.final_result[wt]['r_sigma'])
and np.isclose(final_result[wt]['r_min'], self.final_result[wt]['r_min'])
and np.isclose(final_result[wt]['SD'], self.final_result[wt]['SD'])
and np.isclose(final_result[wt]['SD_sigma'], self.final_result[wt]['SD_sigma']))):
and np.isclose(final_result[wt]['SD_min'], self.final_result[wt]['SD_max']))):
final_result[wt]['is_modified'] = self.final_result[wt]["is_modified"]
self.final_result.update(final_result)
......@@ -2020,24 +2101,31 @@ class MainWindow(QtWidgets.QMainWindow):
MW = phase_result["MW"]
MW_sigma = phase_result.get("MW_sigma", 0)
r = phase_result["r"]
r_sigma = phase_result.get("r_sigma", 0)
#r_sigma = phase_result.get("r_sigma", 0)
r_min = phase_result.get("r_min", np.nan)
r_max = phase_result.get("r_max", np.nan)
SD = phase_result["SD"] * 1E-5
SD_sigma = phase_result.get("SD_sigma", 0) * 1E-5
#SD_sigma = phase_result.get("SD_sigma", 0) * 1E-5
SD_min = phase_result.get("SD_min", np.nan) * 1E-5
SD_max = phase_result.get("SD_max", np.nan) * 1E-5
if phase_result and not '+' in phase:
Qmin = phase_result["Qmin"]
Qmax = phase_result["Qmax"]
Q = phase_result["Q"]
Qmin = phase_result.get("Qmin", np.nan)
Qmax = phase_result.get("Qmax", np.nan)
fc = phase_result["fc"]
fc_sigma = phase_result["fc_sigma"]
#fc_sigma = phase_result["fc_sigma"]
fc_min = phase_result.get("fc_min", np.nan)
fc_max = phase_result.get("fc_max", np.nan)
txt1 += (u'[%s (n=%d)] Mw: %.2f ± %.2f || M₀: %.3e Nm ||'
' Q: %.0f - %.0f || fc: %.1f ± %.1f Hz ||'
' r: %.2f ± %.2f km || SD: %.0f ± %.0f bar\n')
txt1 %= (phase, n, MW, MW_sigma, M0, Qmin, Qmax, fc, fc_sigma,
r, r_sigma, SD, SD_sigma)
' Q: %.0f [%.0f] %.0f || fc: %.1f [%.1f] %.1f Hz ||'
' r: %.2f [%.2f] %.2f km || SD: %.0f [%.0f] %.0f bar\n')
txt1 %= (phase, n, MW, MW_sigma, M0, Qmin, Q, Qmax, fc_min, fc, fc_max,
r_min, r, r_max, SD_min, SD, SD_max)
elif '+' in phase:
txt2 = (u'[%s (n=%d)] Mw: %.2f ± %.2f || M₀: %.3e Nm ||'
'r: %.2f ± %.2f km || SD: %.0f ± %.0f bar')
txt2 %= (phase, n, MW, MW_sigma, M0, r, r_sigma, SD, SD_sigma)
'r: %.2f [%.2f] %.2f km || SD: %.0f [%.0f] %.0f bar')
txt2 %= (phase, n, MW, MW_sigma, M0, r_min, r, r_max, SD_min, SD, SD_max)
self.ui.mean_parameter_line_1.setText(txt1)
self.ui.mean_parameter_line_2.setText(txt2)
......@@ -2113,10 +2201,10 @@ class MainWindow(QtWidgets.QMainWindow):
gmt = self.current_state["spectral_gmt"]
try:
omega_0, f_c, Q, _, _, _ = fit_spectrum(spec, freq, travel_time,
omega_0, f_c, Q, fix_Q=True,
fmin=fmin, fmax=fmax,
sigma=None, model=model,
gmt=gmt)
omega_0, f_c, Q, fix_Q=True,
fmin=fmin, fmax=fmax,
sigma=None, model=model,
gmt=gmt)
except:
print("Error fitting spectrum")
continue
......@@ -2404,7 +2492,10 @@ class MainWindow(QtWidgets.QMainWindow):
:return:
None, :prop:`results` and :prop:`final_result` are modified in place
"""
from waveforms import get_phase_pick
# TODO: selection box if multiple methods / ...
# TODO: show dialog if no solutions are found
# TODO: calculate RMSE if missing
from waveforms import (get_phase_pick, calculate_spectrum)
event = self.current_state["event"]
......@@ -2414,8 +2505,9 @@ class MainWindow(QtWidgets.QMainWindow):
from hazard.belshakelib import BelshakeDB
db = BelshakeDB()
event_ID = int(event.resource_id.id.split("/")[1])
method = "Spectral fit"
event_ID = int(event.resource_id.id.split("/")[2])
method = "sourcespec"
#method = "Spectral fit"
# TODO: empty self.results first?
# or only overwrite existing results for nw_stat_code, comp_code and phase
......@@ -2430,7 +2522,7 @@ class MainWindow(QtWidgets.QMainWindow):
for key in ("M0", "MW", "source_radius", "stress_drop",
"fmin", "fmax", "tmin", "tmax"):
result[key] = mw_params[key]
result["RMSE"] = mw_params["rmse_fit"]
result["RMSE"] = mw_params["rmse_fit"] or np.nan
result["A/M"] = {True: "M", False: "A"}[mw_params["manual"]]
result["use_for_SD"] = mw_params["use_for_SD"]
result["is_modified"] = False
......@@ -2443,13 +2535,15 @@ class MainWindow(QtWidgets.QMainWindow):
for trace in pick.data:
if trace.id[-1] == comp_code:
channel = result["channel"] = trace.id
start_time = trace.stats.starttime
dt = trace.stats.delta
tmin = mw_params["tmin"] - pick.time
tmax = mw_params["tmax"] - pick.time
idx = self.append_result(result, overwrite=overwrite)
if idx is not None:
self.time_windows[phase][channel] = [tmin, tmax]
## Calculate RMSE if missing
if np.isnan(result["RMSE"]):
result["RMSE"] = self._calc_rmse(result, trace, pick)
break
if not phase in phase_component_sel:
......@@ -2476,3 +2570,47 @@ class MainWindow(QtWidgets.QMainWindow):
self._calculate_final_values()
self._update_mean_parameter_display()
db.close()
def _calc_rmse(self, result, trace, pick):
"""
Recalculate RMSE for given result
:param result:
dict, result dictionary
:param trace:
instance of :class:`obspy.Stream`
:param pick:
instance of :class:`robspy.PhasePick`
:return:
float, RMSE
"""
from waveforms import calculate_spectrum
multi_tapered = self.current_state['spectral_calculation'] == 'multi-tapered'
omega_0 = result["omega_0"]
fc = result["corner_frequency"]
Q = result["quality_factor"]
start_time = trace.stats.starttime
dt = trace.stats.delta
tmin = result["tmin"]
tmax = result["tmax"]
imin = int(round((tmin - start_time) / dt))
imax = imin + int(round((tmax - tmin) / dt))
print(imin, imax)
freq, spec, jackknife_errors = calculate_spectrum(trace, imin, imax,
multi_tapered=multi_tapered)
event = self.current_state["event"]
origin = event.origins[0]
travel_time = pick.time - origin.time
fitted_freqs, fitted_spectrum = self.calc_theoretical_spectrum(
omega_0, fc, Q, travel_time)
fmin, fmax = result["fmin"], result["fmax"]
#print(omega_0, fc, Q, travel_time)
rmse = calc_spectral_misfit(spec, freq, fitted_spectrum, fitted_freqs,
fmin, fmax)
return rmse
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment