Skip to content

Commit

Permalink
Remove elog related to updating "sampling_rate" in resample.py
Browse files Browse the repository at this point in the history
  • Loading branch information
Cloud1e committed Oct 24, 2024
1 parent b459473 commit aa8a041
Show file tree
Hide file tree
Showing 2 changed files with 39 additions and 96 deletions.
96 changes: 20 additions & 76 deletions python/mspasspy/algorithms/resample.py
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,8 @@ def resample(self, mspass_object):
# We do this test at the top to avoid having returns testing for
# a dead datum in each of the if conditional blocks below
if isinstance(
mspass_object,
(TimeSeries, Seismogram, TimeSeriesEnsemble, SeismogramEnsemble),
mspass_object,
(TimeSeries, Seismogram, TimeSeriesEnsemble, SeismogramEnsemble),
):
if mspass_object.dead():
return mspass_object
Expand All @@ -221,58 +221,30 @@ def resample(self, mspass_object):

if isinstance(mspass_object, TimeSeries):
data_time_span = (
mspass_object.endtime() - mspass_object.t0 + mspass_object.dt
mspass_object.endtime() - mspass_object.t0 + mspass_object.dt
)
n_resampled = int(data_time_span * self.samprate)
rsdata = signal.resample(
mspass_object.data, n_resampled, window=self.window
)
mspass_object.set_npts(n_resampled)
mspass_object.dt = self.dt
# Check for "sampling_rate" attribute
if mspass_object.is_defined("sampling_rate"):
sampling_rate = mspass_object["sampling_rate"]
# Check if sampling_rate is consistent with 1/dt
if abs(sampling_rate - 1.0 / mspass_object.dt) > 1e-6:
# Record inconsistency in error log (elog)
message = "sampling_rate inconsistent with 1/dt; updating to 1/dt"
mspass_object.elog.log_error("resample",
message,
ErrorSeverity.Complaint)
# Update sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
else:
# Set sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
mspass_object["sampling_rate"] = self.samprate
# We have to go through this conversion to avoid TypeError exceptions
# i.e we can't just copy the entire vector rsdata to the data vector
dv = DoubleVector(rsdata)
mspass_object.data = dv
elif isinstance(mspass_object, Seismogram):
data_time_span = (
mspass_object.endtime() - mspass_object.t0 + mspass_object.dt
mspass_object.endtime() - mspass_object.t0 + mspass_object.dt
)
n_resampled = int(data_time_span * self.samprate)
rsdata = signal.resample(
mspass_object.data, n_resampled, window=self.window, axis=1
)
mspass_object.set_npts(n_resampled)
mspass_object.dt = self.dt
# Check for "sampling_rate" attribute
if mspass_object.is_defined("sampling_rate"):
sampling_rate = mspass_object["sampling_rate"]
# Check if sampling_rate is consistent with 1/dt
if abs(sampling_rate - 1.0 / mspass_object.dt) > 1e-6:
# Record inconsistency in error log (elog)
message = "sampling_rate inconsistent with 1/dt; updating to 1/dt"
mspass_object.elog.log_error("resample",
message,
ErrorSeverity.Complaint)
# Update sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
else:
# Set sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
mspass_object["sampling_rate"] = self.samprate
# We have to go through this conversion to avoid TypeError exceptions
# i.e we can't just copy the entire vector rsdata to the data vector
dm = dmatrix(rsdata)
Expand Down Expand Up @@ -376,8 +348,8 @@ def resample(self, mspass_object):
# We do this test at the top to avoid having returns testing for
# a dead datum in each of the if conditional blocks below
if isinstance(
mspass_object,
(TimeSeries, Seismogram, TimeSeriesEnsemble, SeismogramEnsemble),
mspass_object,
(TimeSeries, Seismogram, TimeSeriesEnsemble, SeismogramEnsemble),
):
if mspass_object.dead():
return mspass_object
Expand Down Expand Up @@ -406,21 +378,7 @@ def resample(self, mspass_object):
dsdata_npts = len(dsdata)
mspass_object.set_npts(dsdata_npts)
mspass_object.dt = self.dt
# Check for "sampling_rate" attribute
if mspass_object.is_defined("sampling_rate"):
sampling_rate = mspass_object["sampling_rate"]
# Check if sampling_rate is consistent with 1/dt
if abs(sampling_rate - 1.0 / mspass_object.dt) > 1e-6:
# Record inconsistency in error log (elog)
message = "sampling_rate inconsistent with 1/dt; updating to 1/dt"
mspass_object.elog.log_error("resample",
message,
ErrorSeverity.Complaint)
# Update sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
else:
# Set sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
mspass_object["sampling_rate"] = self.samprate
# We have to go through this conversion to avoid TypeError exceptions
# i.e we can't just copy the entire vector rsdata to the data vector
mspass_object.data = DoubleVector(dsdata)
Expand Down Expand Up @@ -450,21 +408,7 @@ def resample(self, mspass_object):
dsdata_npts = msize[1]
mspass_object.set_npts(dsdata_npts)
mspass_object.dt = self.dt
# Check for "sampling_rate" attribute
if mspass_object.is_defined("sampling_rate"):
sampling_rate = mspass_object["sampling_rate"]
# Check if sampling_rate is consistent with 1/dt
if abs(sampling_rate - 1.0 / mspass_object.dt) > 1e-6:
# Record inconsistency in error log (elog)
message = "sampling_rate inconsistent with 1/dt; updating to 1/dt"
mspass_object.elog.log_error("resample",
message,
ErrorSeverity.Complaint)
# Update sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
else:
# Set sampling_rate to 1/dt
mspass_object["sampling_rate"] = 1.0 / mspass_object.dt
mspass_object["sampling_rate"] = self.samprate
# We have to go through this conversion to avoid TypeError exceptions
# i.e we can't just copy the entire vector rsdata to the data vector
mspass_object.data = dmatrix(dsdata)
Expand All @@ -482,16 +426,16 @@ def resample(self, mspass_object):

@mspass_func_wrapper
def resample(
mspass_object,
decimator,
resampler,
verify_operators=True,
object_history=False,
alg_name="resample",
alg_id=None,
dryrun=False,
inplace_return=False,
function_return_key=None,
mspass_object,
decimator,
resampler,
verify_operators=True,
object_history=False,
alg_name="resample",
alg_id=None,
dryrun=False,
inplace_return=False,
function_return_key=None,
):
"""
Resample any valid data object to a common sample rate (sample interval).
Expand Down
39 changes: 19 additions & 20 deletions python/tests/algorithms/test_resample.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,7 @@ def test_resample():
upsampler = ScipyResampler(250.0)
tsup = upsampler.resample(ts)
assert np.isclose(tsup.dt, 0.004)
updateSamplingRateMessage = "sampling_rate inconsistent with 1/dt; updating to 1/dt"
assert np.isclose(tsup["sampling_rate"], 250.0) and tsup.elog.get_error_log()[0].algorithm == "resample" and tsup.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(tsup["sampling_rate"], 250.0)
# This computed npts is more robust. Otherwise changes in helper
# would break it
npup = int(ts_npts * 250.0 / 100.0)
Expand All @@ -50,7 +49,7 @@ def test_resample():

tsup = upsampler.resample(ts)
assert np.isclose(tsup.dt, 0.004)
assert ts.is_defined("sampling_rate") and np.isclose(tsup["sampling_rate"], 250.0) and (not tsup.elog.get_error_log() or tsup.elog.get_error_log()[0].message == updateSamplingRateMessage)
assert ts.is_defined("sampling_rate")

# now repeat for downsampling with resample algorithm
ts = TimeSeries(ts0)
Expand All @@ -60,15 +59,15 @@ def test_resample():
# Note plots of this output show the auto antialiasing works as
# advertised in scipy
assert np.isclose(tsds.dt, 0.2)
assert np.isclose(tsds["sampling_rate"], 5.0) and tsds.elog.get_error_log()[0].algorithm == "resample" and tsds.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(tsds["sampling_rate"], 5.0)
assert tsds.npts == int(ts_npts * 5.0 / 100.0)
# Repeat same downsampling with decimate
ts = TimeSeries(ts0)
decimator = ScipyDecimator(5.0)
assert ts.is_defined("sampling_rate")
tsds = decimator.resample(ts)
assert np.isclose(tsds.dt, 0.2)
assert np.isclose(tsds["sampling_rate"], 5.0) and tsds.elog.get_error_log()[0].algorithm == "resample" and tsds.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(tsds["sampling_rate"], 5.0)
# the documentation doesn't tell me why by the scipy decimate
# function seems to round npts up rather than use int
assert tsds.npts == int(ts_npts * 5.0 / 100.0) + 1
Expand All @@ -77,22 +76,22 @@ def test_resample():
assert seis.is_defined("sampling_rate")
seis = upsampler.resample(seis)
assert np.isclose(seis.dt, 0.004)
assert np.isclose(seis["sampling_rate"], 250.0) and seis.elog.get_error_log()[0].algorithm == "resample" and seis.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(seis["sampling_rate"], 250.0)
npup = int(seis0.npts * 250.0 / 20.0)
assert seis.npts == npup
seis = Seismogram(seis0)
assert seis.is_defined("sampling_rate")
seis = ds_resampler.resample(seis)
assert np.isclose(seis.dt, 0.2)
assert np.isclose(seis["sampling_rate"], 5.0) and seis.elog.get_error_log()[0].algorithm == "resample" and seis.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(seis["sampling_rate"], 5.0)
assert seis.npts == int(ts_npts * 5.0 / 20.0)
seis = Seismogram(seis0)
assert seis.is_defined("sampling_rate")
seis = decimator.resample(seis)
# again the round issue noted above
dec_npts = int(seis0.npts * 5.0 / 20.0) + 1
assert np.isclose(seis.dt, 0.2)
assert np.isclose(seis["sampling_rate"], 5.0) and seis.elog.get_error_log()[0].algorithm == "resample" and seis.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(seis["sampling_rate"], 5.0)
assert seis.npts == dec_npts

tse = get_live_timeseries_ensemble(5)
Expand All @@ -104,7 +103,7 @@ def test_resample():
for d in tse.member:
assert d.live
assert np.isclose(d.dt, 0.004)
assert np.isclose(d["sampling_rate"], 250.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 250.0)
assert d.npts == npup

tse = TimeSeriesEnsemble(tse0)
Expand All @@ -114,7 +113,7 @@ def test_resample():
for d in tse.member:
assert d.live
assert np.isclose(d.dt, 0.2)
assert np.isclose(d["sampling_rate"], 5.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 5.0)
assert d.npts == npup

tse = TimeSeriesEnsemble(tse0)
Expand All @@ -123,7 +122,7 @@ def test_resample():
for d in tse.member:
assert d.live
assert np.isclose(d.dt, 0.2)
assert np.isclose(d["sampling_rate"], 5.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 5.0)
assert d.npts == npup

seis_e = get_live_seismogram_ensemble(3)
Expand All @@ -134,7 +133,7 @@ def test_resample():
for d in seis_e.member:
assert d.live
assert np.isclose(d.dt, 0.004)
assert np.isclose(d["sampling_rate"], 250.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 250.0)
assert d.npts == npup

seis_e = SeismogramEnsemble(seis_e0)
Expand All @@ -144,7 +143,7 @@ def test_resample():
for d in seis_e.member:
assert d.live
assert np.isclose(d.dt, 0.2)
assert np.isclose(d["sampling_rate"], 5.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 5.0)
assert d.npts == npup

seis_e = SeismogramEnsemble(seis_e0)
Expand All @@ -154,7 +153,7 @@ def test_resample():
for d in seis_e.member:
assert d.live
assert np.isclose(d.dt, 0.2)
assert np.isclose(d["sampling_rate"], 5.0) and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert np.isclose(d["sampling_rate"], 5.0)
assert d.npts == npup
# Now test resample function. We define two operators
# for 40 sps target
Expand All @@ -163,15 +162,15 @@ def test_resample():
assert ts.is_defined("sampling_rate")
d = resample(ts, decimate40, resample40)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 104
# print(d.dt,d.live,d.npts)
assert ts0.is_defined("sampling_rate")
d = resample(ts0, decimate40, resample40)
# print(d.dt,d.live,d.npts)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 101
assert tse0.member[0].is_defined("sampling_rate")
Expand All @@ -180,7 +179,7 @@ def test_resample():
for d in tse0.member:
# print(d.dt,d.live,d.npts)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 510
assert tse.member[0].is_defined("sampling_rate")
Expand All @@ -189,7 +188,7 @@ def test_resample():
for d in tse0.member:
# print(d.dt,d.live,d.npts)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 510

Expand All @@ -205,7 +204,7 @@ def test_resample():
for d in seis_e.member:
# print(d.dt,d.live,d.npts)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 512
assert seis_e0.member[0].is_defined("sampling_rate")
Expand All @@ -214,7 +213,7 @@ def test_resample():
for d in seis_e0.member:
# print(d.dt,d.live,d.npts)
assert d.dt == 0.025
assert d["sampling_rate"] == 40.0 and d.elog.get_error_log()[0].algorithm == "resample" and d.elog.get_error_log()[0].message == updateSamplingRateMessage
assert d["sampling_rate"] == 40.0
assert d.live
assert d.npts == 510

Expand Down

0 comments on commit aa8a041

Please sign in to comment.