diff --git a/python/mspasspy/algorithms/resample.py b/python/mspasspy/algorithms/resample.py index 1f891c61..1386d074 100755 --- a/python/mspasspy/algorithms/resample.py +++ b/python/mspasspy/algorithms/resample.py @@ -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 @@ -221,7 +221,7 @@ 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( @@ -229,28 +229,14 @@ def resample(self, mspass_object): ) 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( @@ -258,21 +244,7 @@ def resample(self, mspass_object): ) 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) @@ -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 @@ -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) @@ -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) @@ -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). diff --git a/python/tests/algorithms/test_resample.py b/python/tests/algorithms/test_resample.py index ec4520c3..0922e24e 100755 --- a/python/tests/algorithms/test_resample.py +++ b/python/tests/algorithms/test_resample.py @@ -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) @@ -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) @@ -60,7 +59,7 @@ 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) @@ -68,7 +67,7 @@ def test_resample(): 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 @@ -77,14 +76,14 @@ 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") @@ -92,7 +91,7 @@ def test_resample(): # 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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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) @@ -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 @@ -163,7 +162,7 @@ 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) @@ -171,7 +170,7 @@ def test_resample(): 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") @@ -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") @@ -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 @@ -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") @@ -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