From b8297fc603efe030bb76b0e01e2ab6f4cae81424 Mon Sep 17 00:00:00 2001 From: Yiyu Ni Date: Wed, 6 Nov 2024 15:21:36 -0800 Subject: [PATCH] feat: update to noisepy-io v0.2.0 (#327) * enforce single sampling rate for one station Signed-off-by: Yiyu Ni * update tutorial Signed-off-by: Yiyu Ni * update tutorial Signed-off-by: Yiyu Ni * update io ver Signed-off-by: Yiyu Ni * update tests for single_freq Signed-off-by: Yiyu Ni * update get_started Signed-off-by: Yiyu Ni * update monitoring Signed-off-by: Yiyu Ni --------- Signed-off-by: Yiyu Ni --- .gitignore | 3 + README.md | 2 +- pyproject.toml | 2 +- src/noisepy/seis/correlate.py | 25 ++++- tests/test_cross_correlation.py | 33 ++++--- tutorials/cloud/config.yaml | 2 +- tutorials/config.yml | 2 +- tutorials/get_started.ipynb | 34 +++---- tutorials/monitoring/monitoring_demo.ipynb | 104 ++++++++++----------- tutorials/tutorial_compositestore.ipynb | 32 +++---- tutorials/tutorial_ncedc.ipynb | 38 ++++---- tutorials/tutorial_pnwstore.ipynb | 59 ++++++------ tutorials/tutorial_scedc.ipynb | 27 +++--- 13 files changed, 197 insertions(+), 166 deletions(-) diff --git a/.gitignore b/.gitignore index b73e73e5..e4acf05b 100644 --- a/.gitignore +++ b/.gitignore @@ -180,6 +180,9 @@ tutorials/ncedc_data/ tutorials/cli/tmpdata tutorials/pnw_data tutorials/composite_data +tutorials/monitoring/CCF_ASDF/ +tutorials/monitoring/Monitoring_output.csv +tutorials/monitoring/monito_config.yml # Jupyterbook _build/ diff --git a/README.md b/README.md index ed69d8cd..bd853245 100644 --- a/README.md +++ b/README.md @@ -95,7 +95,7 @@ Taxonomy of the NoisePy variables. * ``step`` is the window that get skipped when sliding windows in seconds * ``smooth_N`` number of points for smoothing the time or frequency domain discrete arrays. * ``maxlag`` maximum length in seconds saved in files in each side of the correlation (save on storage) -* ``substack,substack_len`` boolean, window length over which to substack the correlation (to save storage or do monitoring), it has to be a multiple of ``cc_len``. +* ``substack, substack_windows`` boolean, number of window over which to substack the correlation (to save storage or do monitoring). * ``time_chunk, nchunk`` refers to the time unit that defined a single job. for instace, ``cc_len`` is the correlation length (e.g., 1 hour, 30 min), the overall duration of the experiment is the total length (1 month, 1 year, ...). The time chunk could be 1 day: the code would loop through each cc_len window in a for loop. But each day will be sent as a thread. ## Acknowledgements diff --git a/pyproject.toml b/pyproject.toml index 97f39863..34832af1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,7 +50,7 @@ dependencies = [ "PyYAML==6.0", "pydantic-yaml==1.0", "psutil>=5.9.5,<6.0.0", - "noisepy-seis-io>=0.1.16", + "noisepy-seis-io>=0.2.0", "scipy==1.12.0" ] diff --git a/src/noisepy/seis/correlate.py b/src/noisepy/seis/correlate.py index 530ca594..dc099451 100644 --- a/src/noisepy/seis/correlate.py +++ b/src/noisepy/seis/correlate.py @@ -502,20 +502,35 @@ def _filter_channel_data( logging.warning(f"No data available with sampling rate >= {sampling_rate}") return [] if single_freq: - closest_freq = min( - frequencies, - key=lambda f: max(f - sampling_rate, 0), - ) - logger.info(f"Picked {closest_freq} as the closest sampling rate to {sampling_rate}. ") + closest_freq = _get_closest_freq(frequencies, sampling_rate) + logger.info(f"Picked {closest_freq} as the closest sampling_rate to {sampling_rate}. ") filtered_tuples = list(filter(lambda tup: tup[1].sampling_rate == closest_freq, tuples)) logger.info(f"Filtered to {len(filtered_tuples)}/{len(tuples)} channels with sampling rate == {closest_freq}") else: filtered_tuples = list(filter(lambda tup: tup[1].sampling_rate >= sampling_rate, tuples)) + # for each station, pick the closest >= to sampling_rate + tmp = list( + map( + lambda s: [t for t in filtered_tuples if t[0].station == s], + set([t[0].station for t in filtered_tuples]), + ) + ) + filtered_tuples = sum(list(map(lambda t: _filt_single_station(t, sampling_rate), tmp)), []) logger.info(f"Filtered to {len(filtered_tuples)}/{len(tuples)} channels with sampling rate >= {sampling_rate}") return filtered_tuples +def _get_closest_freq(frequencies, sampling_rate: int): + return min(frequencies, key=lambda f: max(f - sampling_rate, 0)) + + +def _filt_single_station(tuples: List[Tuple[Channel, ChannelData]], sampling_rate: int): + frequencies = set(t[1].sampling_rate for t in tuples) + closest_freq = _get_closest_freq(frequencies, sampling_rate) + return [t for t in tuples if t[1].sampling_rate == closest_freq] + + def check_memory(params: ConfigParameters, nsta: int) -> int: # maximum memory allowed # TODO: Is this needed? Should it be configurable? diff --git a/tests/test_cross_correlation.py b/tests/test_cross_correlation.py index d1b926b8..917bfd90 100644 --- a/tests/test_cross_correlation.py +++ b/tests/test_cross_correlation.py @@ -25,21 +25,32 @@ def test_read_channels(): CLOSEST_FREQ = 60 sampling_rate = 40 - freqs = [10, 39, CLOSEST_FREQ, 100] - ch_data = [] - for f in freqs: + ch_data1 = [] + N = 5 + for f in [10, 39, CLOSEST_FREQ, 100]: cd = ChannelData.empty() cd.sampling_rate = f - ch_data.append(cd) - N = 5 - tuples = [(Channel("foo", Station("CI", "bar")), cd) for cd in ch_data] * N + ch_data1.append(cd) + + cd = ChannelData.empty() + cd.sampling_rate = 100 + ch_data2 = [cd] + + tuples = [(Channel("foo", Station("CI", "FOO")), cd) for cd in ch_data1] * N + tuples += [(Channel("bar", Station("CI", "BAR")), cd) for cd in ch_data2] * N - # we should pick the closest frequency that is >= to the target freq, 60 in this case + # we should pick the closest frequency that is >= to the target sampling_rate + # 60 Hz in this case, for both stations + # CI.FOO returns 60 Hz + # CI.BAR returns nothing filtered = _filter_channel_data(tuples, sampling_rate, single_freq=True) assert N == len(filtered) assert [t[1].sampling_rate for t in filtered] == [CLOSEST_FREQ] * N - # we should get all data at >= 40 Hz (60 and 100) + # we should pick the closest frequency that is >= to the target sampling_rate + # but might be different for each station + # CI.FOO returns 60 Hz + # CI.BAR returns 100 Hz filtered = _filter_channel_data(tuples, sampling_rate, single_freq=False) assert N * 2 == len(filtered) assert all([t[1].sampling_rate >= sampling_rate for t in filtered]) @@ -79,11 +90,11 @@ def get_inventory(self, timespan: DateTimeRange, station: Station) -> obspy.Inve @pytest.mark.parametrize("rm_resp", [RmResp.NO, RmResp.INV]) # add tests for SPECTRUM, RESP and POLES_ZEROS @pytest.mark.parametrize("cc_method", [CCMethod.XCORR, CCMethod.COHERENCY, CCMethod.DECONV]) @pytest.mark.parametrize("substack", [True, False]) -@pytest.mark.parametrize("substack_len", [1, 2]) +@pytest.mark.parametrize("substack_windows", [1, 2]) @pytest.mark.parametrize("inc_hours", [0, 24]) @pytest.mark.parametrize("dpath", ["./data/cc", "./data/acc"]) def test_cross_correlation( - rm_resp: RmResp, cc_method: CCMethod, substack: bool, substack_len: int, inc_hours: int, dpath: str + rm_resp: RmResp, cc_method: CCMethod, substack: bool, substack_windows: int, inc_hours: int, dpath: str ): config = ConfigParameters() config.sampling_rate = 1.0 @@ -92,7 +103,7 @@ def test_cross_correlation( config.inc_hours = inc_hours if substack: config.substack = substack - config.substack_len = substack_len * config.cc_len + config.substack_windows = substack_windows path = os.path.join(os.path.dirname(__file__), dpath) raw_store = SCEDCS3DataStore(path, MockCatalogMock()) diff --git a/tutorials/cloud/config.yaml b/tutorials/cloud/config.yaml index 8126edd3..ea87dcfa 100644 --- a/tutorials/cloud/config.yaml +++ b/tutorials/cloud/config.yaml @@ -35,6 +35,6 @@ stationxml: false step: 450.0 storage_options: {} substack: false -substack_len: 1800 +substack_windows: 1 time_norm: no xcorr_only: true diff --git a/tutorials/config.yml b/tutorials/config.yml index d996b3db..7a858421 100644 --- a/tutorials/config.yml +++ b/tutorials/config.yml @@ -15,7 +15,7 @@ cc_method: xcorr smooth_N: 10 smoothspect_N: 10 substack: true -substack_len: 3600 +substack_windows: 1 maxlag: 200 inc_hours: 12 max_over_std: 10 diff --git a/tutorials/get_started.ipynb b/tutorials/get_started.ipynb index 3ec88670..5d3220fd 100644 --- a/tutorials/get_started.ipynb +++ b/tutorials/get_started.ipynb @@ -64,7 +64,7 @@ "from noisepy.seis.io import plotting_modules\n", "from noisepy.seis.io.asdfstore import ASDFRawDataStore, ASDFCCStore, ASDFStackStore\n", "from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, FreqNorm, RmResp, TimeNorm\n", - "from dateutil.parser import isoparse\n", + "from datetime import datetime, timezone\n", "from datetimerange import DateTimeRange\n", "import os\n", "import shutil\n", @@ -95,8 +95,8 @@ "source": [ "config = ConfigParameters() # default config parameters which can be customized\n", "config.inc_hours = 12\n", - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", - "config.cc_len= 3600 # (float) basic unit of data length for fft (sec)\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.cc_len = 3600 # (float) basic unit of data length for fft (sec)\n", " # criteria for data selection\n", "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", @@ -111,32 +111,32 @@ "config.lomax = -115 # max longitude\n", "\n", "# pre-processing parameters\n", - "config.step= 1800.0 # (float) overlapping between each cc_len (sec)\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "config.step = 1800.0 # (float) overlapping between each cc_len (sec)\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "config.freqmin = 0.05\n", "config.freqmax = 2.0\n", "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "# TEMPORAL and SPECTRAL NORMALISATION\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 10 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "config.time_norm = TimeNorm.NO # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", " # TODO: change time_norm option from \"no\" to \"None\"\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", "\n", "# OUTPUTS:\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = 1 # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", "\n", - "config.maxlag= 200 # lags of cross-correlation to save (sec)\n", + "config.maxlag = 200 # lags of cross-correlation to save (sec)\n", "config.substack = True" ] }, @@ -173,9 +173,9 @@ "source": [ "config.networks = [\"*\"]\n", "config.stations = [\"A*\"]\n", - "config.channels = [\"BHE\",\"BHN\",\"BHZ\"]\n", - "config.start_date = isoparse(\"2019-02-01T00:00:00Z\")\n", - "config.end_date = isoparse(\"2019-02-02T00:00:00Z\")\n", + "config.channels = [\"BH?\"]\n", + "config.start_date = datetime(2019, 2, 1, tzinfo=timezone.utc)\n", + "config.end_date = datetime(2019, 2, 2, tzinfo=timezone.utc)\n", "timerange = DateTimeRange(config.start_date, config.end_date)\n", "\n", "# Download data locally. Enters raw data path, channel types, stations, config, and fdsn server.\n", @@ -218,7 +218,7 @@ "file = os.path.join(raw_data_path, \"2019_02_01_00_00_00T2019_02_01_12_00_00.h5\")\n", "raw_store = ASDFRawDataStore(raw_data_path) # Store for reading raw data\n", "timespans = raw_store.get_timespans()\n", - "plotting_modules.plot_waveform(raw_store, timespans[0], 'CI','ADO',0.01,0.4) # this function takes for input: filename, network, station, freqmin, freqmax for a bandpass filter" + "plotting_modules.plot_waveform(raw_store, timespans[0], 'CI','ADO', 0.01, 0.4) # this function takes for input: filename, network, station, freqmin, freqmax for a bandpass filter" ] }, { diff --git a/tutorials/monitoring/monitoring_demo.ipynb b/tutorials/monitoring/monitoring_demo.ipynb index d86d0f78..617239d0 100644 --- a/tutorials/monitoring/monitoring_demo.ipynb +++ b/tutorials/monitoring/monitoring_demo.ipynb @@ -111,16 +111,15 @@ "config.start_date = start_date\n", "config.end_date = end_date\n", "\n", - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", - "config.cc_len= 600 # (int) basic unit of data length for fft (sec)\n", - " # criteria for data selection\n", - "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", - "\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.cc_len = 600 # (int) basic unit of data length for fft (sec)\n", + " \n", + "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", "config.acorr_only = True # only perform auto-correlation or not\n", "config.xcorr_only = False # only perform cross-correlation or not\n", "\n", - "\n", + "# criteria for data selection\n", "config.lamin = 31 # min latitude\n", "config.lamax = 42 # max latitude\n", "config.lomin = -124 # min longitude\n", @@ -129,32 +128,32 @@ "config.stations = [\"LJR\"] # station names, e.g. [\"LJR\",\"DLA\",\"LAF\"]\n", "config.channels = [\"BH?\", \"HH?\"]\n", "\n", - " # pre-processing parameters\n", - "config.step= 600 # (int) overlapping between each cc_len (sec)\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "# pre-processing parameters\n", + "config.step = 600 # (int) overlapping between each cc_len (sec)\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "config.freqmin = 0.05\n", "config.freqmax = 8.0\n", "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "# TEMPORAL and SPECTRAL NORMALISATION\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 1 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "config.time_norm = TimeNorm.ONE_BIT # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", " # TODO: change time_norm option from \"no\" to \"None\"\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", - "config.stack_method=StackMethod.ALL\n", + "config.stack_method = StackMethod.ALL\n", + "\n", "# OUTPUTS:\n", - "num_stack=24\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = num_stack * config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = 24 # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", "\n", "config.maxlag= 60 # lags of cross-correlation to save (sec)" ] @@ -187,9 +186,10 @@ "outputs": [], "source": [ "# S3 data store\n", - "timerange=DateTimeRange(config.start_date, config.end_date)\n", + "timerange = DateTimeRange(config.start_date, config.end_date)\n", "catalog = XMLStationChannelCatalog(S3_STATION_XML, storage_options=S3_STORAGE_OPTIONS) # Station catalog\n", - "raw_store = SCEDCS3DataStore(S3_DATA, catalog, channel_filter(config.networks, config.stations, config.channels), \\\n", + "raw_store = SCEDCS3DataStore(S3_DATA, catalog, \n", + " channel_filter(config.networks, config.stations, config.channels),\n", " timerange, storage_options=S3_STORAGE_OPTIONS) # Store for reading raw data from S3 bucket" ] }, @@ -201,7 +201,7 @@ "source": [ "# CC store \n", "cc_data_path = os.path.join(path, \"CCF_ASDF\")\n", - "cc_store=ASDFCCStore(cc_data_path)\n", + "cc_store = ASDFCCStore(cc_data_path)\n", "\n", "# For this tutorial make sure the previous run is empty\n", "os.system(f\"rm -rf {cc_data_path}\")\n", @@ -209,7 +209,7 @@ "cross_correlate(raw_store, config, cc_store)\n", "\n", "# Save config parameters\n", - "xcorr_config_fn='xcorr_config.yml'\n", + "xcorr_config_fn = 'xcorr_config.yml'\n", "config.save_yaml(xcorr_config_fn)" ] }, @@ -272,8 +272,8 @@ "print('pairs: ',pairs_all)\n", "print('stations: ', stations)\n", "\n", - "src=\"CI.LJR\"\n", - "rec=\"CI.LJR\"\n", + "src = \"CI.LJR\"\n", + "rec = \"CI.LJR\"\n", "timespans = cc_store.get_timespans(src,rec)\n", "print(timespans)\n", "plot_substack_cc(cc_store, timespans[0], 0.1, 1, config.maxlag, False)" @@ -342,7 +342,7 @@ "config_monito.cvel = 2.6 # Rayleigh wave velocities over the freqency bands\n", "config_monito.atten_tbeg = 2.0\n", "config_monito.atten_tend = 10.0\n", - "config_monito.intb_interval_base=0.01 # interval base of intrinsic absorption parameter for a grid-searching process\n" + "config_monito.intb_interval_base = 0.01 # interval base of intrinsic absorption parameter for a grid-searching process\n" ] }, { @@ -450,7 +450,7 @@ }, "outputs": [], "source": [ - "nccomp=config.ncomp**2\n", + "nccomp = config.ncomp**2\n", "cc_array = np.zeros((nccomp*num_chunk * num_segmts, npts_one_segmt), dtype=np.float32)\n", "cc_time = np.zeros( nccomp*num_chunk * num_segmts, dtype=np.float32)\n", "cc_ngood = np.zeros( nccomp*num_chunk * num_segmts, dtype=np.int16)\n", @@ -562,9 +562,9 @@ "outputs": [], "source": [ "# in single-station cross-component case : enz_system = [ \"EN\", \"EZ\", \"NZ\"]\n", - "freq1=config_monito.freq[0]\n", - "freq2=config_monito.freq[1]\n", - "dt=1/config.sampling_rate\n", + "freq1 = config_monito.freq[0]\n", + "freq2 = config_monito.freq[1]\n", + "dt = 1/config.sampling_rate\n", "\n", "## Choose a targeted component\n", "comp = \"EN\"\n", @@ -578,7 +578,7 @@ "# bandpass filter the data.\n", "tcur = np.zeros(shape=(len(indx),npts_one_segmt))\n", "for i in range(len(indx)):\n", - " tcur[i,:]=bandpass(cc_array[indx[i]], freq1, freq2, int(1 / dt), corners=4, zerophase=True)\n", + " tcur[i,:]=bandpass(cc_array[indx[i]], freq1, freq2, int(1 / dt), corners=4, zerophase=True)\n", "\n", "# output stacked data\n", "(\n", @@ -592,7 +592,7 @@ ") = noise_module.stacking(tcur, cc_time[indx], cc_ngood[indx], config)\n", "\n", "# Plot\n", - "fig,ax=plt.subplots(3,1)\n", + "fig, ax = plt.subplots(3,1)\n", "ax[0].imshow(tcur,extent=[-config.maxlag,config.maxlag,iseg,0],aspect='auto',vmin=-0.001,vmax=0.001)\n", "ax[0].set_title(comp)\n", "ax[0].set_xlabel('Lag-time (sec)')\n", @@ -654,8 +654,8 @@ "pcor_cc = np.zeros(shape=(nwin), dtype=np.float32)\n", "ncor_cc = np.zeros(shape=(nwin), dtype=np.float32)\n", "for i in range(nwin):\n", - " pcor_cc[i] = np.corrcoef(tref[pwin_indx], tcur[i, pwin_indx])[0, 1]\n", - " ncor_cc[i] = np.corrcoef(tref[nwin_indx], tcur[i, nwin_indx])[0, 1]\n", + " pcor_cc[i] = np.corrcoef(tref[pwin_indx], tcur[i, pwin_indx])[0, 1]\n", + " ncor_cc[i] = np.corrcoef(tref[nwin_indx], tcur[i, nwin_indx])[0, 1]\n", "\n", "# Plot\n", "plt.figure(figsize=(8,2));plt.grid(True)\n", @@ -698,7 +698,7 @@ "dvv_stretch = np.zeros(shape=(nwin, 4), dtype=np.float32)\n", "\n", "# define the parameters for stretching\n", - "para=dict()\n", + "para = dict()\n", "para[\"t\"] = np.arange(-config.maxlag,config.maxlag+dt,dt)\n", "para[\"freq\"] = [freq1, freq2]\n", "para[\"twin\"] = [config_monito.coda_tbeg, config_monito.coda_tend]\n", @@ -783,7 +783,7 @@ "metadata": {}, "outputs": [], "source": [ - "num_indx=86400//config.substack_len\n", + "num_indx = 86400//config.substack_len\n", "\n", "# select the data from the three cross-component correlations\n", "# that have sufficiently good windows\n", @@ -792,8 +792,8 @@ "indx0 = np.where( (cc_comp.lower() == comp[0].lower()) & cc_ngood==1)[0]\n", "indx1 = np.where( (cc_comp.lower() == comp[1].lower()) & cc_ngood==1)[0]\n", "indx2 = np.where( (cc_comp.lower() == comp[2].lower()) & cc_ngood==1)[0]\n", - "indx=np.intersect1d(np.intersect1d(indx0, indx1-(1*num_indx)),indx2-(3*num_indx))\n", - "nwin=len(indx)\n", + "indx = np.intersect1d(np.intersect1d(indx0, indx1-(1*num_indx)),indx2-(3*num_indx))\n", + "nwin = len(indx)\n", "print(\"Update the window number for the used components :\", nwin)\n", "\n", "# print(comp[0].lower(),indx0)\n", @@ -801,7 +801,7 @@ "# print(comp[2].lower(),indx2)\n", "#print(\"Final common indx between multiple components: \\n\",indx)\n", "\n", - "indx_all=np.zeros(shape=(3,nwin),dtype=np.integer)\n", + "indx_all = np.zeros(shape=(3,nwin),dtype=np.integer)\n", "indx_all[0]=indx\n", "indx_all[1]=indx+(1*num_indx)\n", "indx_all[2]=indx+(3*num_indx)\n" @@ -823,14 +823,14 @@ "metadata": {}, "outputs": [], "source": [ - "print( src_sta, rec_sta, nwin)\n", - "nccomp=len(enz_system)\n", + "print(src_sta, rec_sta, nwin)\n", + "nccomp = len(enz_system)\n", "\n", "# initializing arrays\n", - "all_dvv= np.zeros(shape=(nccomp, nwin), dtype=np.float32)\n", - "all_err= np.zeros(shape=(nccomp, nwin), dtype=np.float32)\n", - "results_dvv= np.zeros(shape=(nwin), dtype=np.float32)\n", - "results_err= np.zeros(shape=(nwin), dtype=np.float32)\n", + "all_dvv = np.zeros(shape=(nccomp, nwin), dtype=np.float32)\n", + "all_err = np.zeros(shape=(nccomp, nwin), dtype=np.float32)\n", + "results_dvv = np.zeros(shape=(nwin), dtype=np.float32)\n", + "results_err = np.zeros(shape=(nwin), dtype=np.float32)\n", " \n", "for icomp in range(0,nccomp):\n", " comp = enz_system[icomp] \n", @@ -840,7 +840,7 @@ " # bandpass filter the data.\n", " tcur = np.zeros(shape=(len(indx),npts_one_segmt))\n", " for i in range(len(indx)):\n", - " tcur[i,:]=bandpass(cc_array[indx[i]], freq1, freq2, int(1 / dt), corners=4, zerophase=True)\n", + " tcur[i,:]=bandpass(cc_array[indx[i]], freq1, freq2, int(1 / dt), corners=4, zerophase=True)\n", " \n", " # output stacked data\n", " (\n", @@ -876,12 +876,12 @@ " all_dvv[icomp]=(dvv_stretch[:, 0]+dvv_stretch[:, 2])/2.0\n", " all_err[icomp]=np.sqrt(dvv_stretch[:, 1]**2+dvv_stretch[:, 3]**2)\n", " print('component: ',comp,' completed. ')\n", - " results_dvv+=all_dvv[icomp]\n", - " results_err+=all_err[icomp]**2\n", - "results_dvv=results_dvv/nccomp\n", - "results_err=np.sqrt(results_err)\n", + " results_dvv += all_dvv[icomp]\n", + " results_err += all_err[icomp]**2\n", + "results_dvv = results_dvv/nccomp\n", + "results_err = np.sqrt(results_err)\n", "print(all_dvv.shape, results_dvv.shape, )\n", - "nwin=len(results_dvv)\n", + "nwin = len(results_dvv)\n", " " ] }, @@ -891,9 +891,9 @@ "metadata": {}, "outputs": [], "source": [ - "fig,ax=plt.subplots(2,2,figsize=(12,4))\n", + "fig, ax = plt.subplots(2, 2, figsize=(12, 4))\n", "ax[0,0].plot(all_dvv.T)\n", - "ax[0,0].legend(enz_system, loc='upper right', bbox_to_anchor=(1.15,1))\n", + "ax[0,0].legend(enz_system, loc='upper right', bbox_to_anchor=(1.15, 1))\n", "ax[0,0].grid(True)\n", "ax[0,0].set_title('Estimated dv/v')\n", "ax[0,0].set_xlabel('Window number')\n", diff --git a/tutorials/tutorial_compositestore.ipynb b/tutorials/tutorial_compositestore.ipynb index 8a417819..c104e50b 100644 --- a/tutorials/tutorial_compositestore.ipynb +++ b/tutorials/tutorial_compositestore.ipynb @@ -169,48 +169,48 @@ " # being a day of data.\n", "\n", "# pre-processing parameters\n", - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", "config.single_freq = False\n", - "config.cc_len= 3600 # (float) basic unit of data length for fft (sec)\n", - "config.step= 1800.0 # (float) overlapping between each cc_len (sec)\n", + "config.cc_len = 3600 # (float) basic unit of data length for fft (sec)\n", + "config.step = 1800.0 # (float) overlapping between each cc_len (sec)\n", "\n", "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", " # If True, the stationXML file is assumed to be provided.\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "\n", "############## NOISE PRE-PROCESSING ##################\n", - "config.freqmin,config.freqmax = 0.05,2.0 # broad band filtering of the data before cross correlation\n", + "config.freqmin, config.freqmax = 0.05, 2.0 # broad band filtering of the data before cross correlation\n", "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "################### SPECTRAL NORMALIZATION ############\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 10 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "#################### TEMPORAL NORMALIZATION ##########\n", "config.time_norm = TimeNorm.ONE_BIT # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", "############ cross correlation ##############\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", "\n", "# OUTPUTS:\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = 1 # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", " # if substack=False, the cross correlation will be stacked over the inc_hour window\n", "\n", "### For monitoring applications ####\n", "## we recommend substacking ever 2-4 cross correlations and storing the substacks\n", "# e.g. \n", "# config.substack = True \n", - "# config.substack_len = 4* config.cc_len\n", + "# config.substack_windows = 4\n", "\n", - "config.maxlag= 200 # lags of cross-correlation to save (sec)\n", + "config.maxlag = 200 # lags of cross-correlation to save (sec)\n", "\n", "# the network list that are used here\n", "config.networks = ['CI', 'NC']" @@ -260,11 +260,11 @@ " storage_options=S3_STORAGE_OPTIONS)\n", "\n", "scedc_store = SCEDCS3DataStore(SCEDC_DATA, scedc_catalog, \n", - " channel_filter([\"CI\"], [\"RPV\", \"SVD\", \"BBR\"], [\"BH?\", \"HH?\"]), \n", + " channel_filter([\"CI\"], [\"VES\", \"SVD\", \"BBR\"], [\"BH?\", \"HH?\"]), \n", " timerange, storage_options=S3_STORAGE_OPTIONS)\n", "\n", "ncedc_store = NCEDCS3DataStore(NCEDC_DATA, ncedc_catalog, \n", - " channel_filter([\"NC\"], [\"KCT\", \"KRP\", \"KHMB\"], [\"BH?\", \"HH?\"]), \n", + " channel_filter([\"NC\"], [\"BBGB\", \"AFD\", \"GDXB\"], [\"BH?\", \"HH?\"]), \n", " timerange, storage_options=S3_STORAGE_OPTIONS)\n", "\n", "raw_store = CompositeRawStore({\"CI\": scedc_store, \n", diff --git a/tutorials/tutorial_ncedc.ipynb b/tutorials/tutorial_ncedc.ipynb index 4e966f3b..2b6e64b4 100644 --- a/tutorials/tutorial_ncedc.ipynb +++ b/tutorials/tutorial_ncedc.ipynb @@ -60,12 +60,12 @@ "source": [ "%load_ext autoreload\n", "%autoreload 2\n", - "from noisepy.seis import cross_correlate, stack_cross_correlations, __version__ # noisepy core functions\n", - "from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore # Object to store ASDF data within noisepy\n", - "from noisepy.seis.io.s3store import NCEDCS3DataStore # Object to query SCEDC data from on S3\n", + "from noisepy.seis import cross_correlate, stack_cross_correlations, __version__ # noisepy core functions\n", + "from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore # Object to store ASDF data within noisepy\n", + "from noisepy.seis.io.s3store import NCEDCS3DataStore # Object to query SCEDC data from on S3\n", "from noisepy.seis.io.channel_filter_store import channel_filter\n", - "from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, FreqNorm, RmResp, StackMethod, TimeNorm # Main configuration object\n", - "from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog # Required stationXML handling object\n", + "from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, FreqNorm, RmResp, StackMethod, TimeNorm\n", + "from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog # Required stationXML handling object\n", "import os\n", "from datetime import datetime, timezone\n", "from datetimerange import DateTimeRange\n", @@ -181,46 +181,46 @@ " # being a day of data.\n", "\n", "# pre-processing parameters\n", - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", - "config.cc_len= 3600 # (float) basic unit of data length for fft (sec)\n", - "config.step= 1800.0 # (float) overlapping between each cc_len (sec)\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.cc_len = 3600 # (float) basic unit of data length for fft (sec)\n", + "config.step = 1800.0 # (float) overlapping between each cc_len (sec)\n", "\n", "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", " # If True, the stationXML file is assumed to be provided.\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "\n", "############## NOISE PRE-PROCESSING ##################\n", - "config.freqmin,config.freqmax = 0.05,2.0 # broad band filtering of the data before cross correlation\n", - "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", + "config.freqmin, config.freqmax = 0.05, 2.0 # broad band filtering of the data before cross correlation\n", + "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "################### SPECTRAL NORMALIZATION ############\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 10 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "#################### TEMPORAL NORMALIZATION ##########\n", "config.time_norm = TimeNorm.ONE_BIT # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", "\n", "############ cross correlation ##############\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", "\n", "# OUTPUTS:\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = 1 # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", " # if substack=False, the cross correlation will be stacked over the inc_hour window\n", "\n", "### For monitoring applications ####\n", "## we recommend substacking ever 2-4 cross correlations and storing the substacks\n", "# e.g. \n", "# config.substack = True \n", - "# config.substack_len = 4* config.cc_len\n", + "# config.substack_windows = 4\n", "\n", "config.maxlag= 200 # lags of cross-correlation to save (sec)" ] diff --git a/tutorials/tutorial_pnwstore.ipynb b/tutorials/tutorial_pnwstore.ipynb index a26472f8..1c698878 100644 --- a/tutorials/tutorial_pnwstore.ipynb +++ b/tutorials/tutorial_pnwstore.ipynb @@ -68,13 +68,13 @@ }, "outputs": [], "source": [ - "from noisepy.seis import cross_correlate, stack_cross_correlations # noisepy core functions\n", + "from noisepy.seis import cross_correlate, stack_cross_correlations # noisepy core functions\n", "from noisepy.seis.io import plotting_modules\n", - "from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore # Object to store ASDF data within noisepy\n", + "from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore # Object to store ASDF data within noisepy\n", "from noisepy.seis.io.channel_filter_store import channel_filter\n", "from noisepy.seis.io.pnwstore import PNWDataStore\n", - "from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, Channel, ChannelData, ChannelType, FreqNorm, RmResp, Station, TimeNorm # Main configuration object\n", - "from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog # Required stationXML handling object\n", + "from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, Channel, ChannelData, ChannelType, FreqNorm, RmResp, Station, TimeNorm\n", + "from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog # Required stationXML handling object\n", "import os\n", "from datetime import datetime\n", "from datetimerange import DateTimeRange\n", @@ -92,9 +92,12 @@ "metadata": {}, "outputs": [], "source": [ - "STATION_XML = \"/1-fnp/pnwstore1/p-wd11/PNWStationXML/\" # storage of stationXML\n", - "DATA = \"/1-fnp/pnwstore1/p-wd00/PNW2020/__/\" # __ indicates any two chars (network code)\n", - "DB_PATH=\"/1-fnp/pnwstore1/p-wd00/PNW2020/timeseries.sqlite\"\n", + "# storage of stationXML\n", + "STATION_XML = \"/1-fnp/pnwstore1/p-wd11/PNWStationXML/\" \n", + "\n", + "# __ indicates any two chars (network code). However, it will slow down the query\n", + "DATA = \"/1-fnp/pnwstore1/p-wd00/PNW2020/__/\"\n", + "DB_PATH =\"/1-fnp/pnwstore1/p-wd00/PNW2020/timeseries.sqlite\"\n", "# timeframe for analysis\n", "start = datetime(2020, 1, 2)\n", "end = datetime(2020, 1, 4)\n", @@ -160,8 +163,8 @@ }, "outputs": [], "source": [ - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", - "config.cc_len= 3600 # (float) basic unit of data length for fft (sec)\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.cc_len = 3600 # (float) basic unit of data length for fft (sec)\n", " # criteria for data selection\n", "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", @@ -170,33 +173,33 @@ "\n", "# config.inc_hours = 24 # if the data is first \n", "\n", - " # pre-processing parameters\n", - "config.step= 1800.0 # (float) overlapping between each cc_len (sec)\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "# pre-processing parameters\n", + "config.step = 1800.0 # (float) overlapping between each cc_len (sec)\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "config.freqmin = 0.05\n", "config.freqmax = 2.0\n", "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "# TEMPORAL and SPECTRAL NORMALISATION\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitenning or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 10 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "config.time_norm = TimeNorm.NO # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", " # TODO: change time_norm option from \"no\" to \"None\"\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", "\n", "# OUTPUTS:\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = 1 # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", "\n", - "config.maxlag= 200 # lags of cross-correlation to save (sec)\n", + "config.maxlag = 200 # lags of cross-correlation to save (sec)\n", "config.substack = True" ] }, @@ -226,16 +229,14 @@ "metadata": {}, "outputs": [], "source": [ - "# cross network, cross channel type\n", - "# UW.BBO..HH\n", - "# UW.BABR..HH\n", - "# UW.SHUK..BH\n", - "# CC.PANH..BH\n", - "stations = \"BBO,BABR,SHUK,PANH\".split(\",\") # filter to these stations\n", + "config.networks = [\"UW\", \"UO\", \"PB\", \"CC\"]\n", + "config.stations = [\"BBO\", \"BABR\", \"SHUK\", \"PANH\"]\n", + "config.channels = [\"BH?\", \"HH?\"]\n", + "\n", "catalog = XMLStationChannelCatalog(STATION_XML, path_format=\"{network}/{network}.{name}.xml\")\n", - "raw_store = PNWDataStore(DATA, catalog, DB_PATH, \n", - " channel_filter([\"UW\", \"UO\", \"PB\", \"CC\"], stations, \n", - " [\"BH?\", \"HH?\"]), date_range=timerange) # Store for reading raw data from S3 bucket\n", + "raw_store = PNWDataStore(DATA, DB_PATH, catalog,\n", + " channel_filter(config.networks, config.stations, config.channels), \n", + " date_range=timerange) # Store for reading raw data from S3 bucket\n", "cc_store = ASDFCCStore(cc_data_path) # Store for writing CC data" ] }, diff --git a/tutorials/tutorial_scedc.ipynb b/tutorials/tutorial_scedc.ipynb index 4deebcfa..85e5996f 100644 --- a/tutorials/tutorial_scedc.ipynb +++ b/tutorials/tutorial_scedc.ipynb @@ -182,47 +182,48 @@ " # being a day of data.\n", "\n", "# pre-processing parameters\n", - "config.sampling_rate= 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", - "config.cc_len= 3600 # (float) basic unit of data length for fft (sec)\n", - "config.step= 1800.0 # (float) overlapping between each cc_len (sec)\n", + "config.sampling_rate = 20 # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)\n", + "config.single_freq = False\n", + "config.cc_len = 3600 # (float) basic unit of data length for fft (sec)\n", + "config.step = 1800.0 # (float) overlapping between each cc_len (sec)\n", "\n", "config.ncomp = 3 # 1 or 3 component data (needed to decide whether do rotation)\n", "\n", - "config.stationxml= False # station.XML file used to remove instrument response for SAC/miniseed data\n", + "config.stationxml = False # station.XML file used to remove instrument response for SAC/miniseed data\n", " # If True, the stationXML file is assumed to be provided.\n", - "config.rm_resp= RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", + "config.rm_resp = RmResp.INV # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',\n", "\n", "############## NOISE PRE-PROCESSING ##################\n", "config.freqmin,config.freqmax = 0.05,2.0 # broad band filtering of the data before cross correlation\n", "config.max_over_std = 10 # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them\n", "\n", "################### SPECTRAL NORMALIZATION ############\n", - "config.freq_norm= FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", + "config.freq_norm = FreqNorm.RMA # choose between \"rma\" for a soft whitening or \"no\" for no whitening. Pure whitening is not implemented correctly at this point.\n", "config.smoothspect_N = 10 # moving window length to smooth spectrum amplitude (points)\n", " # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)\n", "\n", "#################### TEMPORAL NORMALIZATION ##########\n", "config.time_norm = TimeNorm.ONE_BIT # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,\n", - "config.smooth_N= 10 # moving window length for time domain normalization if selected (points)\n", + "config.smooth_N = 10 # moving window length for time domain normalization if selected (points)\n", "\n", "############ cross correlation ##############\n", - "config.cc_method= CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", + "config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;\n", " # FOR \"COHERENCY\" PLEASE set freq_norm to \"rma\", time_norm to \"no\" and cc_method to \"xcorr\"\n", "\n", "# OUTPUTS:\n", "config.substack = True # True = smaller stacks within the time chunk. False: it will stack over inc_hours\n", - "config.substack_len = config.cc_len # how long to stack over (for monitoring purpose): need to be multiples of cc_len\n", - " # if substack=True, substack_len=2*cc_len, then you pre-stack every 2 correlation windows.\n", - " # for instance: substack=True, substack_len=cc_len means that you keep ALL of the correlations\n", + "config.substack_windows = config.cc_len # how long to stack over (for monitoring purpose)\n", + " # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.\n", + " # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations\n", " # if substack=False, the cross correlation will be stacked over the inc_hour window\n", "\n", "### For monitoring applications ####\n", "## we recommend substacking ever 2-4 cross correlations and storing the substacks\n", "# e.g. \n", "# config.substack = True \n", - "# config.substack_len = 4* config.cc_len\n", + "# config.substack_windows = 4\n", "\n", - "config.maxlag= 200 # lags of cross-correlation to save (sec)" + "config.maxlag = 200 # lags of cross-correlation to save (sec)" ] }, {