Python PPSD.add方法代码示例

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_ppsd_w_iris_against_obspy_results(self):
        """
        Test against results obtained after merging of #1108.
        """
        # Read in ANMO data for one day
        st = read(os.path.join(self.path, "IUANMO.seed"))

        # Read in metadata in various different formats
        paz = {
            "gain": 86298.5,
            "zeros": [0, 0],
            "poles": [-59.4313, -22.7121 + 27.1065j, -22.7121 + 27.1065j, -0.0048004, -0.073199],
            "sensitivity": 3.3554 * 10 ** 9,
        }
        resp = os.path.join(self.path, "IUANMO.resp")
        parser = Parser(os.path.join(self.path, "IUANMO.dataless"))
        inv = read_inventory(os.path.join(self.path, "IUANMO.xml"))

        # load expected results, for both only PAZ and full response
        filename_paz = os.path.join(self.path, "IUANMO_ppsd_paz.npz")
        results_paz = PPSD.load_npz(filename_paz, metadata=None)
        filename_full = os.path.join(self.path, "IUANMO_ppsd_fullresponse.npz")
        results_full = PPSD.load_npz(filename_full, metadata=None)

        # Calculate the PPSDs and test against expected results
        # first: only PAZ
        ppsd = PPSD(st[0].stats, paz)
        ppsd.add(st)
        # commented code to generate the test data:
        # ## np.savez(filename_paz,
        # ##          **dict([(k, getattr(ppsd, k))
        # ##                  for k in PPSD.NPZ_STORE_KEYS]))
        for key in PPSD.NPZ_STORE_KEYS_ARRAY_TYPES:
            np.testing.assert_allclose(getattr(ppsd, key), getattr(results_paz, key), rtol=1e-5)
        for key in PPSD.NPZ_STORE_KEYS_LIST_TYPES:
            for got, expected in zip(getattr(ppsd, key), getattr(results_paz, key)):
                np.testing.assert_allclose(got, expected, rtol=1e-5)
        for key in PPSD.NPZ_STORE_KEYS_SIMPLE_TYPES:
            if key in ["obspy_version", "numpy_version", "matplotlib_version"]:
                continue
            self.assertEqual(getattr(ppsd, key), getattr(results_paz, key))
        # second: various methods for full response
        for metadata in [parser, inv, resp]:
            ppsd = PPSD(st[0].stats, metadata)
            ppsd.add(st)
            # commented code to generate the test data:
            # ## np.savez(filename_full,
            # ##          **dict([(k, getattr(ppsd, k))
            # ##                  for k in PPSD.NPZ_STORE_KEYS]))
            for key in PPSD.NPZ_STORE_KEYS_ARRAY_TYPES:
                np.testing.assert_allclose(getattr(ppsd, key), getattr(results_full, key), rtol=1e-5)
            for key in PPSD.NPZ_STORE_KEYS_LIST_TYPES:
                for got, expected in zip(getattr(ppsd, key), getattr(results_full, key)):
                    np.testing.assert_allclose(got, expected, rtol=1e-5)
            for key in PPSD.NPZ_STORE_KEYS_SIMPLE_TYPES:
                if key in ["obspy_version", "numpy_version", "matplotlib_version"]:
                    continue
                self.assertEqual(getattr(ppsd, key), getattr(results_full, key))

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
def _get_ppsd():
    """
    Returns ready computed ppsd for testing purposes.
    """
    tr, paz = _get_sample_data()
    st = Stream([tr])
    ppsd = PPSD(tr.stats, paz, db_bins=(-200, -50, 0.5))
    ppsd.add(st)
    return ppsd

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_PPSD_w_IRIS(self):
        # Bands to be used this is the upper and lower frequency band pairs
        fres = zip([0.1, 0.05], [0.2, 0.1])

        file_dataANMO = os.path.join(self.path, 'IUANMO.seed')
        # Read in ANMO data for one day
        st = read(file_dataANMO)

        # Use a canned ANMO response which will stay static
        paz = {'gain': 86298.5, 'zeros': [0, 0],
               'poles': [-59.4313, -22.7121 + 27.1065j, -22.7121 + 27.1065j,
                         -0.0048004, -0.073199], 'sensitivity': 3.3554*10**9}

        # Make an empty PPSD and add the data
        ppsd = PPSD(st[0].stats, paz)
        ppsd.add(st)
        ppsd.calculate_histogram()

        # Get the 50th percentile from the PPSD
        (per, perval) = ppsd.get_percentile(percentile=50)

        # Read in the results obtained from a Mustang flat file
        file_dataIRIS = os.path.join(self.path, 'IRISpdfExample')
        freq, power, hits = np.genfromtxt(file_dataIRIS, comments='#',
                                          delimiter=',', unpack=True)

        # For each frequency pair we want to compare the mean of the bands
        for fre in fres:
            pervalGoodOBSPY = []

            # Get the values for the bands from the PPSD
            perinv = 1 / per
            mask = (fre[0] < perinv) & (perinv < fre[1])
            pervalGoodOBSPY = perval[mask]

            # Now we sort out all of the data from the IRIS flat file
            mask = (fre[0] < freq) & (freq < fre[1])
            triples = list(zip(freq[mask], hits[mask], power[mask]))
            # We now have all of the frequency values of interest
            # We will get the distinct frequency values
            freqdistinct = sorted(list(set(freq[mask])), reverse=True)
            percenlist = []
            # We will loop through the frequency values and compute a
            # 50th percentile
            for curfreq in freqdistinct:
                tempvalslist = []
                for triple in triples:
                    if np.isclose(curfreq, triple[0], atol=1e-3, rtol=0.0):
                        tempvalslist += [int(triple[2])] * int(triple[1])
                percenlist.append(np.percentile(tempvalslist, 50))
            # Here is the actual test
            np.testing.assert_allclose(np.mean(pervalGoodOBSPY),
                                       np.mean(percenlist), rtol=0.0, atol=1.0)

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_PPSD_w_IRIS_against_obspy_results(self):
        """
        Test against results obtained after merging of #1108.
        """
        # Read in ANMO data for one day
        st = read(os.path.join(self.path, 'IUANMO.seed'))

        # Read in metadata in various different formats
        paz = {'gain': 86298.5, 'zeros': [0, 0],
               'poles': [-59.4313, -22.7121 + 27.1065j, -22.7121 + 27.1065j,
                         -0.0048004, -0.073199], 'sensitivity': 3.3554*10**9}
        resp = os.path.join(self.path, 'IUANMO.resp')
        parser = Parser(os.path.join(self.path, 'IUANMO.dataless'))
        inv = read_inventory(os.path.join(self.path, 'IUANMO.xml'))

        # load expected results, for both only PAZ and full response
        filename_paz = os.path.join(self.path, 'IUANMO_ppsd_paz.npz')
        results_paz = np.load(filename_paz)
        filename_full = os.path.join(self.path, 'IUANMO_ppsd_fullresponse.npz')
        results_full = np.load(filename_full)
        arrays_to_check = ['_times_data', '_times_processed', '_times_gaps',
                           '_spec_octaves', 'per_octaves', 'per_octaves_right',
                           'per_octaves_left', 'period_bin_centers',
                           'spec_bins', 'period_bins']
        arrays_to_check = [native_str(key) for key in arrays_to_check]

        # Calculate the PPSDs and test against expected results
        # first: only PAZ
        ppsd = PPSD(st[0].stats, paz)
        ppsd.add(st)
        # commented code to generate the test data:
        # ## np.savez(filename_paz,
        # ##          **dict([(k, getattr(ppsd, k)) for k in arrays_to_check]))
        for key in arrays_to_check:
            self.assertTrue(np.allclose(
                getattr(ppsd, key), results_paz[key], rtol=1e-5))
        # second: various methods for full response
        # (also test various means of initialization, basically testing the
        #  decorator that maps the deprecated keywords)
        for metadata in [parser, inv, resp]:
            ppsd = PPSD(st[0].stats, paz=metadata)
            ppsd = PPSD(st[0].stats, parser=metadata)
            ppsd = PPSD(st[0].stats, metadata)
            ppsd.add(st)
            # commented code to generate the test data:
            # ## np.savez(filename_full,
            # ##          **dict([(k, getattr(ppsd, k))
            # ##                  for k in arrays_to_check]))
            for key in arrays_to_check:
                self.assertTrue(np.allclose(
                    getattr(ppsd, key), results_full[key], rtol=1e-5))

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_exclude_last_sample(self):
        start = UTCDateTime("2017-01-01T00:00:00")
        header = {
            "starttime": start,
            "network": "GR",
            "station": "FUR",
            "channel": "BHZ"
        }
        # 49 segments of 30 minutes to allow 30 minutes overlap in next day
        tr = Trace(data=np.arange(30 * 60 * 4, dtype=np.int32), header=header)

        ppsd = PPSD(tr.stats, read_inventory())
        ppsd.add(tr)

        self.assertEqual(3, len(ppsd._times_processed))
        self.assertEqual(3600, ppsd.len)
        for i, time in enumerate(ppsd._times_processed):
            current = start.ns + (i * 30 * 60) * 1e9
            self.assertTrue(time == current)

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
 def test_issue1216(self):
     tr, paz = _get_sample_data()
     st = Stream([tr])
     ppsd = PPSD(tr.stats, paz, db_bins=(-200, -50, 0.5))
     ppsd.add(st)
     # After adding data internal representation of hist stack is None
     self.assertIsNone(ppsd._current_hist_stack)
     # Accessing the current_histogram property calculates the default stack
     self.assertIsNotNone(ppsd.current_histogram)
     self.assertIsNotNone(ppsd._current_hist_stack)
     # Adding the same data again does not invalidate the internal stack
     # but raises "UserWarning: Already covered time spans detected"
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter('always', UserWarning)
         ppsd.add(st)
         self.assertEqual(len(w), 4)
         for w_ in w:
             self.assertTrue(str(w_.message).startswith(
                 "Already covered time spans detected"))
     self.assertIsNotNone(ppsd._current_hist_stack)
     # Adding new data invalidates the internal stack
     tr.stats.starttime += 3600
     st2 = Stream([tr])
     # raises "UserWarning: Already covered time spans detected"
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter('always', UserWarning)
         ppsd.add(st2)
         self.assertEqual(len(w), 2)
         for w_ in w:
             self.assertTrue(str(w_.message).startswith(
                 "Already covered time spans detected"))
     self.assertIsNone(ppsd._current_hist_stack)
     # Accessing current_histogram again calculates the stack
     self.assertIsNotNone(ppsd.current_histogram)
     self.assertIsNotNone(ppsd._current_hist_stack)

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
 def test_wrong_trace_id_message(self):
     """
     Test that we get the expected warning message on waveform/metadata
     mismatch.
     """
     tr, _paz = _get_sample_data()
     inv = read_inventory(os.path.join(self.path, 'IUANMO.xml'))
     st = Stream([tr])
     ppsd = PPSD(tr.stats, inv)
     # metadata doesn't fit the trace ID specified via stats
     # should show a warning..
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter('always')
         ret = ppsd.add(st)
         # the trace is sliced into four segments, so we get the warning
         # message four times..
         self.assertEqual(len(w), 4)
         for w_ in w:
             self.assertTrue(str(w_.message).startswith(
                 "Error getting response from provided metadata"))
     # should not add the data to the ppsd
     self.assertFalse(ret)

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_ppsd_w_iris(self):
        # Bands to be used this is the upper and lower frequency band pairs
        fres = zip([0.1, 0.05], [0.2, 0.1])

        file_data_anmo = os.path.join(self.path, 'IUANMO.seed')
        # Read in ANMO data for one day
        st = read(file_data_anmo)

        # Use a canned ANMO response which will stay static
        paz = {'gain': 86298.5, 'zeros': [0, 0],
               'poles': [-59.4313, -22.7121 + 27.1065j, -22.7121 + 27.1065j,
                         -0.0048004, -0.073199],
               'sensitivity': 3.3554 * 10 ** 9}

        # Make an empty PPSD and add the data
        # use highest frequency given by IRIS Mustang noise-pdf web service
        # (0.475683 Hz == 2.10224036 s) as center of first bin, so that we
        # end up with the same bins.
        ppsd = PPSD(st[0].stats, paz, period_limits=(2.10224036, 1400))
        ppsd.add(st)
        ppsd.calculate_histogram()

        # Get the 50th percentile from the PPSD
        (per, perval) = ppsd.get_percentile(percentile=50)
        perinv = 1 / per

        # Read in the results obtained from a Mustang flat file
        file_data_iris = os.path.join(self.path, 'IRISpdfExample')
        data = np.genfromtxt(
            file_data_iris, comments='#', delimiter=',',
            dtype=[(native_str("freq"), np.float64),
                   (native_str("power"), np.int32),
                   (native_str("hits"), np.int32)])
        freq = data["freq"]
        power = data["power"]
        hits = data["hits"]
        # cut data to same period range as in the ppsd we computed
        # (Mustang returns more long periods, probably due to some zero padding
        # or longer nfft in psd)
        num_periods = len(ppsd.period_bin_centers)
        freqdistinct = np.array(sorted(set(freq), reverse=True)[:num_periods])
        # just make sure that we compare the same periods in the following
        # (as we access both frequency arrays by indices from now on)
        np.testing.assert_allclose(freqdistinct, 1 / ppsd.period_bin_centers,
                                   rtol=1e-4)

        # For each frequency pair we want to compare the mean of the bands
        for fre in fres:
            # determine which bins we want to compare
            mask = (fre[0] < perinv) & (perinv < fre[1])

            # Get the values for the bands from the PPSD
            per_val_good_obspy = perval[mask]

            percenlist = []
            # Now we sort out all of the data from the IRIS flat file
            # We will loop through the frequency values and compute a
            # 50th percentile
            for curfreq in freqdistinct[mask]:
                mask_ = curfreq == freq
                tempvalslist = np.repeat(power[mask_], hits[mask_])
                percenlist.append(np.percentile(tempvalslist, 50))
            # Here is the actual test
            np.testing.assert_allclose(np.mean(per_val_good_obspy),
                                       np.mean(percenlist), rtol=0.0, atol=1.2)

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_PPSD(self):
        """
        Test PPSD routine with some real data. Data was downsampled to 100Hz
        so the ppsd is a bit distorted which does not matter for the purpose
        of testing.
        """
        # load test file
        file_data = os.path.join(
            self.path, 'BW.KW1._.EHZ.D.2011.090_downsampled.asc.gz')
        file_histogram = os.path.join(
            self.path,
            'BW.KW1._.EHZ.D.2011.090_downsampled__ppsd_hist_stack.npy')
        file_binning = os.path.join(
            self.path, 'BW.KW1._.EHZ.D.2011.090_downsampled__ppsd_mixed.npz')
        # parameters for the test
        data = np.loadtxt(file_data)
        stats = {'_format': 'MSEED',
                 'calib': 1.0,
                 'channel': 'EHZ',
                 'delta': 0.01,
                 'endtime': UTCDateTime(2011, 3, 31, 2, 36, 0, 180000),
                 'location': '',
                 'mseed': {'dataquality': 'D', 'record_length': 512,
                           'encoding': 'STEIM2', 'byteorder': '>'},
                 'network': 'BW',
                 'npts': 936001,
                 'sampling_rate': 100.0,
                 'starttime': UTCDateTime(2011, 3, 31, 0, 0, 0, 180000),
                 'station': 'KW1'}
        tr = Trace(data, stats)
        st = Stream([tr])
        paz = {'gain': 60077000.0,
               'poles': [(-0.037004 + 0.037016j), (-0.037004 - 0.037016j),
                         (-251.33 + 0j), (-131.04 - 467.29j),
                         (-131.04 + 467.29j)],
               'sensitivity': 2516778400.0,
               'zeros': [0j, 0j]}
        ppsd = PPSD(tr.stats, paz)
        ppsd.add(st)
        # read results and compare
        result_hist = np.load(file_histogram)
        self.assertEqual(len(ppsd.times), 4)
        self.assertEqual(ppsd.nfft, 65536)
        self.assertEqual(ppsd.nlap, 49152)
        np.testing.assert_array_equal(ppsd.hist_stack, result_hist)
        # add the same data a second time (which should do nothing at all) and
        # test again - but it will raise UserWarnings, which we omit for now
        with warnings.catch_warnings(record=True):
            warnings.simplefilter('ignore', UserWarning)
            ppsd.add(st)
            np.testing.assert_array_equal(ppsd.hist_stack, result_hist)
        # test the binning arrays
        binning = np.load(file_binning)
        np.testing.assert_array_equal(ppsd.spec_bins, binning['spec_bins'])
        np.testing.assert_array_equal(ppsd.period_bins, binning['period_bins'])

        # test saving and loading of the PPSD (using a temporary file)
        with NamedTemporaryFile() as tf:
            filename = tf.name
            # test saving and loading an uncompressed file
            ppsd.save(filename, compress=False)
            ppsd_loaded = PPSD.load(filename)
            self.assertEqual(len(ppsd_loaded.times), 4)
            self.assertEqual(ppsd_loaded.nfft, 65536)
            self.assertEqual(ppsd_loaded.nlap, 49152)
            np.testing.assert_array_equal(ppsd_loaded.hist_stack, result_hist)
            np.testing.assert_array_equal(ppsd_loaded.spec_bins,
                                          binning['spec_bins'])
            np.testing.assert_array_equal(ppsd_loaded.period_bins,
                                          binning['period_bins'])
            # test saving and loading a compressed file
            ppsd.save(filename, compress=True)
            ppsd_loaded = PPSD.load(filename)
            self.assertEqual(len(ppsd_loaded.times), 4)
            self.assertEqual(ppsd_loaded.nfft, 65536)
            self.assertEqual(ppsd_loaded.nlap, 49152)
            np.testing.assert_array_equal(ppsd_loaded.hist_stack, result_hist)
            np.testing.assert_array_equal(ppsd_loaded.spec_bins,
                                          binning['spec_bins'])
            np.testing.assert_array_equal(ppsd_loaded.period_bins,
                                          binning['period_bins'])

# 需要导入模块: from obspy.signal.spectral_estimation import PPSD [as 别名]
# 或者: from obspy.signal.spectral_estimation.PPSD import add [as 别名]
    def test_ppsd_w_iris_against_obspy_results(self):
        """
        Test against results obtained after merging of #1108.
        """
        # Read in ANMO data for one day
        st = read(os.path.join(self.path, 'IUANMO.seed'))

        # Read in metadata in various different formats
        paz = {'gain': 86298.5, 'zeros': [0, 0],
               'poles': [-59.4313, -22.7121 + 27.1065j, -22.7121 + 27.1065j,
                         -0.0048004, -0.073199], 'sensitivity': 3.3554*10**9}
        resp = os.path.join(self.path, 'IUANMO.resp')
        parser = Parser(os.path.join(self.path, 'IUANMO.dataless'))
        inv = read_inventory(os.path.join(self.path, 'IUANMO.xml'))

        # load expected results, for both only PAZ and full response
        filename_paz = os.path.join(self.path, 'IUANMO_ppsd_paz.npz')
        results_paz = PPSD.load_npz(filename_paz, metadata=None)
        filename_full = os.path.join(self.path,
                                     'IUANMO_ppsd_fullresponse.npz')
        results_full = PPSD.load_npz(filename_full, metadata=None)

        # Calculate the PPSDs and test against expected results
        # first: only PAZ
        ppsd = PPSD(st[0].stats, paz)
        ppsd.add(st)
        # commented code to generate the test data:
        # ## np.savez(filename_paz,
        # ##          **dict([(k, getattr(ppsd, k))
        # ##                  for k in PPSD.NPZ_STORE_KEYS]))
        for key in PPSD.NPZ_STORE_KEYS_ARRAY_TYPES:
            np.testing.assert_allclose(
                getattr(ppsd, key), getattr(results_paz, key), rtol=1e-5)
        for key in PPSD.NPZ_STORE_KEYS_LIST_TYPES:
            for got, expected in zip(getattr(ppsd, key),
                                     getattr(results_paz, key)):
                np.testing.assert_allclose(got, expected, rtol=1e-5)
        for key in PPSD.NPZ_STORE_KEYS_SIMPLE_TYPES:
            if key in ["obspy_version", "numpy_version", "matplotlib_version"]:
                continue
            self.assertEqual(getattr(ppsd, key), getattr(results_paz, key))
        # second: various methods for full response
        # (also test various means of initialization, basically testing the
        #  decorator that maps the deprecated keywords)
        for metadata in [parser, inv, resp]:
            with warnings.catch_warnings(record=True) as w:
                warnings.simplefilter('always')
                ppsd = PPSD(st[0].stats, paz=metadata)
            self.assertEqual(len(w), 1)
            self.assertIs(w[0].category, ObsPyDeprecationWarning)

            with warnings.catch_warnings(record=True) as w:
                warnings.simplefilter('always')
                ppsd = PPSD(st[0].stats, parser=metadata)
            self.assertEqual(len(w), 1)
            self.assertIs(w[0].category, ObsPyDeprecationWarning)

            ppsd = PPSD(st[0].stats, metadata)
            ppsd.add(st)
            # commented code to generate the test data:
            # ## np.savez(filename_full,
            # ##          **dict([(k, getattr(ppsd, k))
            # ##                  for k in PPSD.NPZ_STORE_KEYS]))
            for key in PPSD.NPZ_STORE_KEYS_ARRAY_TYPES:
                np.testing.assert_allclose(
                    getattr(ppsd, key), getattr(results_full, key), rtol=1e-5)
            for key in PPSD.NPZ_STORE_KEYS_LIST_TYPES:
                for got, expected in zip(getattr(ppsd, key),
                                         getattr(results_full, key)):
                    np.testing.assert_allclose(got, expected, rtol=1e-5)
            for key in PPSD.NPZ_STORE_KEYS_SIMPLE_TYPES:
                if key in ["obspy_version", "numpy_version",
                           "matplotlib_version"]:
                    continue
                self.assertEqual(getattr(ppsd, key),
                                 getattr(results_full, key))