Composite Datastore Tutorial

Composite Datastore Tutorial#

Noisepy is a python software package to process ambient seismic noise cross correlations. This tutorial aims to introduce the use of noisepy for a toy problem on a composite data store. It can be ran locally or on the cloud.

The data is stored on AWS S3:

SCEDC Data Set: https://scedc.caltech.edu/data/getstarted-pds.html
NCEDC Data Set: https://ncedc.org/db/cloud/getstarted-pds.html

First, we install the noisepy-seis package

# Uncomment and run this line if the environment doesn't have noisepy already installed:
# ! pip install noisepy-seis 

Warning: NoisePy uses obspy as a core Python module to manipulate seismic data. If you use Google Colab, restart the runtime now for proper installation of obspy on Colab.

Import necessary modules#

Then we import the basic modules

%load_ext autoreload
%autoreload 2
from noisepy.seis import cross_correlate, stack_cross_correlations, __version__    # noisepy core functions
from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore                  # Object to store ASDF data within noisepy
from noisepy.seis.io.compositerawstore import CompositeRawStore
from noisepy.seis.io.s3store import SCEDCS3DataStore, NCEDCS3DataStore
from noisepy.seis.io.channel_filter_store import channel_filter
from noisepy.seis.io.datatypes import Channel, CCMethod, ConfigParameters, FreqNorm, RmResp, StackMethod, TimeNorm        # Main configuration object
from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog                # Required stationXML handling object
import os
import obspy
import shutil
from datetime import datetime, timezone
from datetimerange import DateTimeRange

from noisepy.seis.io.plotting_modules import plot_all_moveout

print(f"Using NoisePy version {__version__}")

path = "./composite_data" 

os.makedirs(path, exist_ok=True)
cc_data_path = os.path.join(path, "CCF")
stack_data_path = os.path.join(path, "STACK")
S3_STORAGE_OPTIONS = {"s3": {"anon": True}}

/opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/noisepy/seis/io/utils.py:13: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

Using NoisePy version 0.1.dev1

We will work with a single day worth of data on SCEDC and NCEDC. The continuous data is organized with a single day and channel per miniseed. For this example, you can choose any year since 2002. We will just cross correlate a single day.

# SCEDC S3 bucket common URL characters for that day.
SCEDC_DATA = "s3://scedc-pds/continuous_waveforms/"
NCEDC_DATA = "s3://ncedc-pds/continuous_waveforms/NC/"

SCEDC_STATION_XML = "s3://scedc-pds/FDSNstationXML/CI/"  
NCEDC_STATION_XML = "s3://ncedc-pds/FDSNstationXML/NC/"

# timeframe for analysis
start = datetime(2012, 1, 1, tzinfo=timezone.utc)
end = datetime(2012, 1, 3, tzinfo=timezone.utc)
timerange = DateTimeRange(start, end)
print(timerange)

2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000

Ambient Noise Project Configuration#

We prepare the configuration of the workflow by declaring and storing parameters into the ConfigParameters() object and/or editing the config.yml file.

# Initialize ambient noise workflow configuration
config = ConfigParameters() # default config parameters which can be customized

Customize the job parameters below:

config.start_date = start
config.end_date = end
config.acorr_only = False # only perform auto-correlation or not
config.xcorr_only = True # only perform cross-correlation or not

config.inc_hours = 24 # INC_HOURS is used in hours (integer) as the 
        #chunk of time that the paralelliztion will work.
        # data will be loaded in memory, so reduce memory with smaller 
        # inc_hours if there are over 400+ stations.
        # At regional scale for SCEDC, we can afford 20Hz data and inc_hour 
        # being a day of data.

# pre-processing parameters
config.sampling_rate = 20  # (int) Sampling rate in Hz of desired processing (it can be different than the data sampling rate)
config.single_freq = False
config.cc_len = 3600  # (float) basic unit of data length for fft (sec)
config.step = 1800.0  # (float) overlapping between each cc_len (sec)

config.ncomp = 3  # 1 or 3 component data (needed to decide whether do rotation)

config.stationxml = False  # station.XML file used to remove instrument response for SAC/miniseed data
      # If True, the stationXML file is assumed to be provided.
config.rm_resp = RmResp.INV  # select 'no' to not remove response and use 'inv' if you use the stationXML,'spectrum',

############## NOISE PRE-PROCESSING ##################
config.freqmin, config.freqmax = 0.05, 2.0  # broad band filtering of the data before cross correlation
config.max_over_std  = 10  # threshold to remove window of bad signals: set it to 10*9 if prefer not to remove them

################### SPECTRAL NORMALIZATION ############
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.
config.smoothspect_N = 10  # moving window length to smooth spectrum amplitude (points)
    # here, choose smoothspect_N for the case of a strict whitening (e.g., phase_only)

#################### TEMPORAL NORMALIZATION ##########
config.time_norm = TimeNorm.ONE_BIT # 'no' for no normalization, or 'rma', 'one_bit' for normalization in time domain,
config.smooth_N = 10  # moving window length for time domain normalization if selected (points)

############ cross correlation ##############
config.cc_method = CCMethod.XCORR # 'xcorr' for pure cross correlation OR 'deconv' for deconvolution;
    # FOR "COHERENCY" PLEASE set freq_norm to "rma", time_norm to "no" and cc_method to "xcorr"

# OUTPUTS:
config.substack = True  # True = smaller stacks within the time chunk. False: it will stack over inc_hours
config.substack_windows = 1  # how long to stack over (for monitoring purpose)
    # if substack=True, substack_windows=2, then you pre-stack every 2 correlation windows.
    # for instance: substack=True, substack_windows=1 means that you keep ALL of the correlations
    # if substack=False, the cross correlation will be stacked over the inc_hour window

### For monitoring applications ####
## we recommend substacking ever 2-4 cross correlations and storing the substacks
# e.g. 
# config.substack = True 
# config.substack_windows = 4

config.maxlag = 200  # lags of cross-correlation to save (sec)

# the network list that are used here
config.networks = ['CI', 'NC']

# For this tutorial make sure the previous run is empty
#os.system(f"rm -rf {cc_data_path}")
if os.path.exists(cc_data_path):
    shutil.rmtree(cc_data_path)

Step 1: Cross-correlation#

In this instance, we read directly the data from the scedc bucket into the cross correlation code. We are not attempting to recreate a data store. Therefore we go straight to step 1 of the cross correlations.

We first declare the data and cross correlation stores

scedc_catalog = XMLStationChannelCatalog(SCEDC_STATION_XML, 
                                         storage_options=S3_STORAGE_OPTIONS)
ncedc_catalog = XMLStationChannelCatalog(NCEDC_STATION_XML, "{network}.{name}.xml", 
                                         storage_options=S3_STORAGE_OPTIONS)

scedc_store = SCEDCS3DataStore(SCEDC_DATA, scedc_catalog,  
                               channel_filter(["CI"], ["VES", "SVD", "BBR"], ["BH?", "HH?"]), 
                               timerange, storage_options=S3_STORAGE_OPTIONS)

ncedc_store = NCEDCS3DataStore(NCEDC_DATA, ncedc_catalog, 
                               channel_filter(["NC"], ["BBGB", "AFD", "GDXB"], ["BH?", "HH?"]), 
                               timerange, storage_options=S3_STORAGE_OPTIONS)

raw_store = CompositeRawStore({"CI": scedc_store, 
                               "NC": ncedc_store}) # Composite store for reading data from both SCEDC and NCEDC
cc_store = ASDFCCStore(cc_data_path)               # Store for writing CC data

get the time range of the data in the data store inventory

span = raw_store.get_timespans()
print(span)

[2012-01-01T00:00:00+0000 - 2012-01-02T00:00:00+0000, 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000]

Get the channels available during a given time spane

channel_list=raw_store.get_channels(span[1])
print(channel_list)

2025-01-14 01:07:36,164 140435400354688 INFO utils.log_raw(): TIMING:  0.558 secs for Listing 1731 files from s3://scedc-pds/continuous_waveforms/2012/2012_002/

2025-01-14 01:07:36,199 140435400354688 INFO utils.log_raw(): TIMING:  0.034 secs for Init: 1 timespans and 18 channels

2025-01-14 01:07:36,332 140433818650304 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://scedc-pds/FDSNstationXML/CI/CI_BBR.xml

2025-01-14 01:07:36,544 140433780901568 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://scedc-pds/FDSNstationXML/CI/CI_SVD.xml

2025-01-14 01:07:36,557 140433799775936 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://scedc-pds/FDSNstationXML/CI/CI_VES.xml

2025-01-14 01:07:40,066 140435400354688 INFO s3store.get_channels(): Getting 18 channels for 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:07:40,278 140435400354688 INFO utils.log_raw(): TIMING:  0.187 secs for Listing 804 files from s3://ncedc-pds/continuous_waveforms/NC/2012/2012.002/

2025-01-14 01:07:40,294 140435400354688 INFO utils.log_raw(): TIMING:  0.017 secs for Init: 1 timespans and 9 channels

2025-01-14 01:07:40,355 140433799775936 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://ncedc-pds/FDSNstationXML/NC/NC.GDXB.xml

2025-01-14 01:07:40,378 140433590060736 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://ncedc-pds/FDSNstationXML/NC/NC.AFD.xml

2025-01-14 01:07:40,383 140433571186368 INFO channelcatalog._get_inventory_from_file(): Reading StationXML file s3://ncedc-pds/FDSNstationXML/NC/NC.BBGB.xml

2025-01-14 01:07:40,713 140435400354688 INFO s3store.get_channels(): Getting 9 channels for 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

[CI.BBR.BHE, CI.BBR.BHN, CI.BBR.BHZ, CI.BBR.HHE, CI.BBR.HHN, CI.BBR.HHZ, CI.SVD.BHE, CI.SVD.BHN, CI.SVD.BHZ, CI.SVD.HHE, CI.SVD.HHN, CI.SVD.HHZ, CI.VES.BHE, CI.VES.BHN, CI.VES.BHZ, CI.VES.HHE, CI.VES.HHN, CI.VES.HHZ, NC.AFD.HHE, NC.AFD.HHN, NC.AFD.HHZ, NC.BBGB.HHE, NC.BBGB.HHN, NC.BBGB.HHZ, NC.GDXB.HHE, NC.GDXB.HHN, NC.GDXB.HHZ]

Perform the cross correlation#

The data will be pulled from SCEDC & NCEDC, cross correlated, and stored locally if this notebook is ran locally. If you are re-calculating, we recommend to clear the old cc_store.

# we also define a channel pair filter that limits cross-correlation to stations pairs closer than 600 km
def pair_filter(src: Channel, rec: Channel) -> bool:
    latS = src.station.lat
    lonS = src.station.lon
    latR = rec.station.lat
    lonR = rec.station.lon
    dist, _, _ = obspy.geodetics.base.gps2dist_azimuth(latS, lonS, latR, lonR)
    dist /= 1e3 # to km
    if dist <= 600:
        return True
    else:
        return False

cross_correlate(raw_store, config, cc_store, pair_filter=pair_filter)

2025-01-14 01:07:40,766 140435400354688 INFO correlate.cross_correlate(): Starting Cross-Correlation with 4 cores

2025-01-14 01:07:41,349 140435400354688 INFO utils.log_raw(): TIMING:  0.582 secs for Listing 1731 files from s3://scedc-pds/continuous_waveforms/2012/2012_001/

2025-01-14 01:07:41,382 140435400354688 INFO utils.log_raw(): TIMING:  0.033 secs for Init: 2 timespans and 36 channels

2025-01-14 01:07:41,383 140435400354688 INFO s3store.get_channels(): Getting 18 channels for 2012-01-01T00:00:00+0000 - 2012-01-02T00:00:00+0000

2025-01-14 01:07:41,534 140435400354688 INFO utils.log_raw(): TIMING:  0.142 secs for Listing 804 files from s3://ncedc-pds/continuous_waveforms/NC/2012/2012.001/

2025-01-14 01:07:41,551 140435400354688 INFO utils.log_raw(): TIMING:  0.017 secs for Init: 2 timespans and 18 channels

2025-01-14 01:07:41,552 140435400354688 INFO s3store.get_channels(): Getting 9 channels for 2012-01-01T00:00:00+0000 - 2012-01-02T00:00:00+0000

2025-01-14 01:07:41,555 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.788 secs for get 27 channels

2025-01-14 01:07:41,564 140435400354688 INFO correlate.cc_timespan(): Checking for stations already done: 17 pairs

2025-01-14 01:07:41,566 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.010 secs for check for 6 stations already done (warm up cache)

2025-01-14 01:07:41,570 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.004 secs for check for stations already done

2025-01-14 01:07:41,571 140435400354688 INFO correlate.cc_timespan(): Still need to process: 6/6 stations, 27/27 channels, 17/17 pairs for 2012-01-01T00:00:00+0000 - 2012-01-02T00:00:00+0000

2025-01-14 01:07:43,948 140435400354688 INFO correlate._filter_channel_data(): Filtered to 18/27 channels with sampling rate >= 20.0

2025-01-14 01:07:43,952 140435400354688 INFO utils.log_raw(): TIMING CC Main:  2.382 secs for Read channel data: 18 channels

2025-01-14 01:07:45,816 140433195796160 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:46,014 140433397122752 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:46,113 140433378248384 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:47,459 140433415997120 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:47,503 140433571186368 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:47,589 140433214670528 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:47,762 140433590060736 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:47,788 140433176921792 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:52,293 140433378248384 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:52,310 140433397122752 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:53,539 140433195796160 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:58,895 140433378248384 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:59,028 140433571186368 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:59,193 140433176921792 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:07:59,425 140433214670528 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:00,572 140433590060736 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHE | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:00,641 140433397122752 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHZ | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:00,710 140433415997120 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHN | 2012-01-01T00:00:00.000000Z - 2012-01-01T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.

 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.
 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.

2025-01-14 01:08:06,682 140435400354688 INFO utils.log_raw(): TIMING CC Main: 22.730 secs for Preprocess: 18 channels

2025-01-14 01:08:06,683 140435400354688 INFO correlate.check_memory(): Require  0.22gb memory for cross correlations

2025-01-14 01:08:08,282 140435400354688 INFO utils.log_raw(): TIMING CC Main:  1.598 secs for Compute FFTs: 18 channels

2025-01-14 01:08:08,288 140435400354688 INFO correlate.cc_timespan(): Starting CC with 17 station pairs

2025-01-14 01:08:13,451 140435400354688 INFO utils.log_raw(): TIMING CC Main:  5.163 secs for Correlate and write to store

2025-01-14 01:08:13,624 140435400354688 INFO utils.log_raw(): TIMING CC Main: 32.857 secs for Process the chunk of 2012-01-01T00:00:00+0000 - 2012-01-02T00:00:00+0000

2025-01-14 01:08:13,637 140435400354688 INFO s3store.get_channels(): Getting 18 channels for 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:08:13,645 140435400354688 INFO s3store.get_channels(): Getting 9 channels for 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:08:13,648 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.012 secs for get 27 channels

2025-01-14 01:08:13,658 140435400354688 INFO correlate.cc_timespan(): Checking for stations already done: 17 pairs

2025-01-14 01:08:13,660 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.011 secs for check for 6 stations already done (warm up cache)

2025-01-14 01:08:13,662 140435400354688 INFO utils.log_raw(): TIMING CC Main:  0.002 secs for check for stations already done

2025-01-14 01:08:13,662 140435400354688 INFO correlate.cc_timespan(): Still need to process: 6/6 stations, 27/27 channels, 17/17 pairs for 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:08:16,267 140435400354688 INFO correlate._filter_channel_data(): Filtered to 18/27 channels with sampling rate >= 20.0

2025-01-14 01:08:16,270 140435400354688 INFO utils.log_raw(): TIMING CC Main:  2.608 secs for Read channel data: 18 channels

2025-01-14 01:08:18,058 140433590060736 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:18,121 140433214670528 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:18,299 140433176921792 INFO noise_module.preprocess_raw(): removing response for CI.BBR..BHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:19,667 140433195796160 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:19,684 140433378248384 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:19,904 140433397122752 INFO noise_module.preprocess_raw(): removing response for NC.BBGB..HHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:20,053 140433415997120 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:20,099 140433571186368 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:24,487 140433214670528 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:24,702 140433176921792 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:25,667 140433590060736 INFO noise_module.preprocess_raw(): removing response for NC.AFD..HHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:31,053 140433214670528 INFO noise_module.preprocess_raw(): removing response for CI.VES..BHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:31,374 140433397122752 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:31,469 140433571186368 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:31,690 140433415997120 INFO noise_module.preprocess_raw(): removing response for CI.SVD..BHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

2025-01-14 01:08:32,249 140433378248384 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHN | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.

2025-01-14 01:08:32,558 140433195796160 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHE | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.
2025-01-14 01:08:32,776 140433176921792 INFO noise_module.preprocess_raw(): removing response for NC.GDXB..HHZ | 2012-01-02T00:00:00.000000Z - 2012-01-02T23:59:59.950000Z | 20.0 Hz, 1728000 samples using inv

 WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent. 
	 Execution continuing.

2025-01-14 01:08:38,850 140435400354688 INFO utils.log_raw(): TIMING CC Main: 22.580 secs for Preprocess: 18 channels

2025-01-14 01:08:38,851 140435400354688 INFO correlate.check_memory(): Require  0.22gb memory for cross correlations

2025-01-14 01:08:40,707 140435400354688 INFO utils.log_raw(): TIMING CC Main:  1.855 secs for Compute FFTs: 18 channels

2025-01-14 01:08:40,713 140435400354688 INFO correlate.cc_timespan(): Starting CC with 17 station pairs

2025-01-14 01:08:45,639 140435400354688 INFO utils.log_raw(): TIMING CC Main:  4.926 secs for Correlate and write to store

2025-01-14 01:08:45,807 140435400354688 INFO utils.log_raw(): TIMING CC Main: 32.171 secs for Process the chunk of 2012-01-02T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:08:45,817 140435400354688 INFO utils.log_raw(): TIMING CC Main: 65.051 secs for Step 1 in total with 4 cores

The cross correlations are saved as a single file for each channel pair and each increment of inc_hours. We now will stack all the cross correlations over all time chunk and look at all station pairs results.

Step 2: Stack the cross correlation#

We now create the stack stores. Because this tutorial runs locally, we will use an ASDF stack store to output the data. ASDF is a data container in HDF5 that is used in full waveform modeling and inversion. H5 behaves well locally.

Each station pair will have 1 H5 file with all components of the cross correlations. While this produces many H5 files, it has come down to the noisepy team’s favorite option:

the thread-safe installation of HDF5 is not trivial
the choice of grouping station pairs within a single file is not flexible to a broad audience of users.

# open a new cc store in read-only mode since we will be doing parallel access for stacking
cc_store = ASDFCCStore(cc_data_path, mode="r")
stack_store = ASDFStackStore(stack_data_path)

Configure the stacking#

There are various methods for optimal stacking. We refern to Yang et al (2022) for a discussion and comparison of the methods:

Yang X, Bryan J, Okubo K, Jiang C, Clements T, Denolle MA. Optimal stacking of noise cross-correlation functions. Geophysical Journal International. 2023 Mar;232(3):1600-18. https://doi.org/10.1093/gji/ggac410

Users have the choice to implement linear, phase weighted stacks pws (Schimmel et al, 1997), robust stacking (Yang et al, 2022), acf autocovariance filter (Nakata et al, 2019), nroot stacking. Users may choose the stacking method of their choice by entering the strings in config.stack_method.

If chosen all, the current code only ouputs linear, pws, robust since nroot is less common and acf is computationally expensive.

config.stack_method = StackMethod.LINEAR

method_list = [method for method in dir(StackMethod) if not method.startswith("__")]
print(method_list)

['ALL', 'AUTO_COVARIANCE', 'LINEAR', 'NROOT', 'PWS', 'ROBUST', 'SELECTIVE']

cc_store.get_station_pairs()

[(CI.BBR, CI.SVD),
 (CI.VES, NC.GDXB),
 (NC.AFD, NC.GDXB),
 (CI.BBR, NC.BBGB),
 (CI.BBR, CI.VES),
 (NC.AFD, NC.AFD),
 (NC.AFD, CI.VES),
 (NC.BBGB, NC.AFD),
 (CI.VES, CI.VES),
 (NC.BBGB, CI.SVD),
 (CI.VES, CI.SVD),
 (CI.BBR, CI.BBR),
 (NC.BBGB, NC.BBGB),
 (NC.BBGB, CI.VES),
 (CI.SVD, CI.SVD),
 (NC.BBGB, NC.GDXB),
 (NC.GDXB, NC.GDXB)]

stack_cross_correlations(cc_store, stack_store, config)

2025-01-14 01:08:45,953 140435400354688 INFO stack.initializer(): Station pairs: 17

2025-01-14 01:08:49,541 140051491511168 INFO stack.stack_store_pair(): Stacking CI.BBR_CI.VES/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,541 139880003394432 INFO stack.stack_store_pair(): Stacking CI.BBR_CI.SVD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,573 139626855713664 INFO stack.stack_store_pair(): Stacking CI.BBR_CI.BBR/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,620 140051491511168 INFO utils.log_raw(): TIMING:  0.075 secs for loading CCF data
2025-01-14 01:08:49,628 139880003394432 INFO utils.log_raw(): TIMING:  0.082 secs for loading CCF data
2025-01-14 01:08:49,632 139626855713664 INFO utils.log_raw(): TIMING:  0.053 secs for loading CCF data
2025-01-14 01:08:49,638 140234178608000 INFO stack.stack_store_pair(): Stacking CI.BBR_NC.BBGB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,640 140051491511168 INFO utils.log_raw(): TIMING:  0.020 secs for stack/rotate all station pairs (CI.BBR, CI.VES)
2025-01-14 01:08:49,646 139626855713664 INFO utils.log_raw(): TIMING:  0.014 secs for stack/rotate all station pairs (CI.BBR, CI.BBR)
2025-01-14 01:08:49,647 139880003394432 INFO utils.log_raw(): TIMING:  0.019 secs for stack/rotate all station pairs (CI.BBR, CI.SVD)
2025-01-14 01:08:49,664 139626855713664 INFO utils.log_raw(): TIMING:  0.018 secs for writing stack pair (CI.BBR, CI.BBR)
2025-01-14 01:08:49,665 139626855713664 INFO stack.stack_store_pair(): Stacking CI.SVD_CI.SVD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,705 139880003394432 INFO utils.log_raw(): TIMING:  0.058 secs for writing stack pair (CI.BBR, CI.SVD)
2025-01-14 01:08:49,706 139880003394432 INFO stack.stack_store_pair(): Stacking CI.VES_CI.SVD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,711 140051491511168 INFO utils.log_raw(): TIMING:  0.068 secs for writing stack pair (CI.BBR, CI.VES)
2025-01-14 01:08:49,712 140051491511168 INFO stack.stack_store_pair(): Stacking CI.VES_CI.VES/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,719 139626855713664 INFO utils.log_raw(): TIMING:  0.052 secs for loading CCF data
2025-01-14 01:08:49,731 140234178608000 INFO utils.log_raw(): TIMING:  0.087 secs for loading CCF data
2025-01-14 01:08:49,732 139626855713664 INFO utils.log_raw(): TIMING:  0.013 secs for stack/rotate all station pairs (CI.SVD, CI.SVD)

2025-01-14 01:08:49,749 140234178608000 INFO utils.log_raw(): TIMING:  0.018 secs for stack/rotate all station pairs (CI.BBR, NC.BBGB)
2025-01-14 01:08:49,752 139626855713664 INFO utils.log_raw(): TIMING:  0.019 secs for writing stack pair (CI.SVD, CI.SVD)
2025-01-14 01:08:49,752 139626855713664 INFO stack.stack_store_pair(): Stacking CI.VES_NC.GDXB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,760 140051491511168 INFO utils.log_raw(): TIMING:  0.045 secs for loading CCF data
2025-01-14 01:08:49,771 140051491511168 INFO utils.log_raw(): TIMING:  0.011 secs for stack/rotate all station pairs (CI.VES, CI.VES)
2025-01-14 01:08:49,784 139880003394432 INFO utils.log_raw(): TIMING:  0.074 secs for loading CCF data
2025-01-14 01:08:49,790 140051491511168 INFO utils.log_raw(): TIMING:  0.018 secs for writing stack pair (CI.VES, CI.VES)
2025-01-14 01:08:49,791 140051491511168 INFO stack.stack_store_pair(): Stacking NC.AFD_CI.VES/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,803 139880003394432 INFO utils.log_raw(): TIMING:  0.018 secs for stack/rotate all station pairs (CI.VES, CI.SVD)
2025-01-14 01:08:49,811 140234178608000 INFO utils.log_raw(): TIMING:  0.061 secs for writing stack pair (CI.BBR, NC.BBGB)
2025-01-14 01:08:49,812 140234178608000 INFO stack.stack_store_pair(): Stacking NC.AFD_NC.AFD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,832 139626855713664 INFO utils.log_raw(): TIMING:  0.075 secs for loading CCF data
2025-01-14 01:08:49,851 139626855713664 INFO utils.log_raw(): TIMING:  0.019 secs for stack/rotate all station pairs (CI.VES, NC.GDXB)
2025-01-14 01:08:49,856 139880003394432 INFO utils.log_raw(): TIMING:  0.053 secs for writing stack pair (CI.VES, CI.SVD)
2025-01-14 01:08:49,857 139880003394432 INFO stack.stack_store_pair(): Stacking NC.AFD_NC.GDXB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,860 140051491511168 INFO utils.log_raw(): TIMING:  0.066 secs for loading CCF data
2025-01-14 01:08:49,865 140234178608000 INFO utils.log_raw(): TIMING:  0.047 secs for loading CCF data
2025-01-14 01:08:49,877 140234178608000 INFO utils.log_raw(): TIMING:  0.012 secs for stack/rotate all station pairs (NC.AFD, NC.AFD)
2025-01-14 01:08:49,879 140051491511168 INFO utils.log_raw(): TIMING:  0.019 secs for stack/rotate all station pairs (NC.AFD, CI.VES)
2025-01-14 01:08:49,896 140234178608000 INFO utils.log_raw(): TIMING:  0.019 secs for writing stack pair (NC.AFD, NC.AFD)
2025-01-14 01:08:49,898 140234178608000 INFO stack.stack_store_pair(): Stacking NC.BBGB_CI.SVD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,907 139626855713664 INFO utils.log_raw(): TIMING:  0.055 secs for writing stack pair (CI.VES, NC.GDXB)
2025-01-14 01:08:49,907 139626855713664 INFO stack.stack_store_pair(): Stacking NC.BBGB_CI.VES/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:49,920 139880003394432 INFO utils.log_raw(): TIMING:  0.060 secs for loading CCF data
2025-01-14 01:08:49,936 139880003394432 INFO utils.log_raw(): TIMING:  0.016 secs for stack/rotate all station pairs (NC.AFD, NC.GDXB)
2025-01-14 01:08:49,939 140051491511168 INFO utils.log_raw(): TIMING:  0.060 secs for writing stack pair (NC.AFD, CI.VES)
2025-01-14 01:08:49,940 140051491511168 INFO stack.stack_store_pair(): Stacking NC.BBGB_NC.AFD/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000

2025-01-14 01:08:49,973 139626855713664 INFO utils.log_raw(): TIMING:  0.062 secs for loading CCF data
2025-01-14 01:08:49,976 140234178608000 INFO utils.log_raw(): TIMING:  0.074 secs for loading CCF data
2025-01-14 01:08:49,989 139626855713664 INFO utils.log_raw(): TIMING:  0.016 secs for stack/rotate all station pairs (NC.BBGB, CI.VES)
2025-01-14 01:08:49,991 140234178608000 INFO utils.log_raw(): TIMING:  0.015 secs for stack/rotate all station pairs (NC.BBGB, CI.SVD)
2025-01-14 01:08:49,992 139880003394432 INFO utils.log_raw(): TIMING:  0.056 secs for writing stack pair (NC.AFD, NC.GDXB)
2025-01-14 01:08:49,995 139880003394432 INFO stack.stack_store_pair(): Stacking NC.BBGB_NC.BBGB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:50,008 140051491511168 INFO utils.log_raw(): TIMING:  0.064 secs for loading CCF data
2025-01-14 01:08:50,023 140051491511168 INFO utils.log_raw(): TIMING:  0.016 secs for stack/rotate all station pairs (NC.BBGB, NC.AFD)
2025-01-14 01:08:50,040 139880003394432 INFO utils.log_raw(): TIMING:  0.042 secs for loading CCF data
2025-01-14 01:08:50,048 139626855713664 INFO utils.log_raw(): TIMING:  0.058 secs for writing stack pair (NC.BBGB, CI.VES)
2025-01-14 01:08:50,049 139626855713664 INFO stack.stack_store_pair(): Stacking NC.BBGB_NC.GDXB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:50,051 139880003394432 INFO utils.log_raw(): TIMING:  0.012 secs for stack/rotate all station pairs (NC.BBGB, NC.BBGB)
2025-01-14 01:08:50,053 140234178608000 INFO utils.log_raw(): TIMING:  0.062 secs for writing stack pair (NC.BBGB, CI.SVD)
2025-01-14 01:08:50,056 140234178608000 INFO stack.stack_store_pair(): Stacking NC.GDXB_NC.GDXB/2012-01-01T00:00:00+0000 - 2012-01-03T00:00:00+0000
2025-01-14 01:08:50,071 139880003394432 INFO utils.log_raw(): TIMING:  0.019 secs for writing stack pair (NC.BBGB, NC.BBGB)
2025-01-14 01:08:50,081 140051491511168 INFO utils.log_raw(): TIMING:  0.057 secs for writing stack pair (NC.BBGB, NC.AFD)
2025-01-14 01:08:50,102 140234178608000 INFO utils.log_raw(): TIMING:  0.042 secs for loading CCF data
2025-01-14 01:08:50,108 139626855713664 INFO utils.log_raw(): TIMING:  0.055 secs for loading CCF data
2025-01-14 01:08:50,110 140234178608000 INFO utils.log_raw(): TIMING:  0.008 secs for stack/rotate all station pairs (NC.GDXB, NC.GDXB)
2025-01-14 01:08:50,118 139626855713664 INFO utils.log_raw(): TIMING:  0.010 secs for stack/rotate all station pairs (NC.BBGB, NC.GDXB)
2025-01-14 01:08:50,122 140234178608000 INFO utils.log_raw(): TIMING:  0.011 secs for writing stack pair (NC.GDXB, NC.GDXB)
2025-01-14 01:08:50,150 139626855713664 INFO utils.log_raw(): TIMING:  0.031 secs for writing stack pair (NC.BBGB, NC.GDXB)

2025-01-14 01:08:50,724 140435400354688 INFO utils.log_raw(): TIMING:  4.773 secs for step 2 in total

QC_1 of the cross correlations for Imaging#

We now explore the quality of the cross correlations. We plot the moveout of the cross correlations, filtered in some frequency band.

cc_store.get_station_pairs()

[(CI.BBR, CI.SVD),
 (CI.VES, NC.GDXB),
 (NC.AFD, NC.GDXB),
 (CI.BBR, NC.BBGB),
 (CI.BBR, CI.VES),
 (NC.AFD, NC.AFD),
 (NC.AFD, CI.VES),
 (NC.BBGB, NC.AFD),
 (CI.VES, CI.VES),
 (NC.BBGB, CI.SVD),
 (CI.VES, CI.SVD),
 (CI.BBR, CI.BBR),
 (NC.BBGB, NC.BBGB),
 (NC.BBGB, CI.VES),
 (CI.SVD, CI.SVD),
 (NC.BBGB, NC.GDXB),
 (NC.GDXB, NC.GDXB)]

pairs = stack_store.get_station_pairs()
print(f"Found {len(pairs)} station pairs")
sta_stacks = stack_store.read_bulk(timerange, pairs) # no timestamp used in ASDFStackStore

Found 17 station pairs

2025-01-14 01:08:51,148 140435400354688 INFO utils.log_raw(): TIMING:  0.370 secs for loading 17 stacks

plot_all_moveout(sta_stacks, 'Allstack_linear', 0.1, 0.2, 'ZZ', 1)

2025-01-14 01:08:51,174 140435400354688 INFO plotting_modules.plot_all_moveout(): Plottting: Allstack_linear, 17 station pairs

_images/3b558c574cd298c94797e2210e5561351f35b0d294f3c7507502b7da4302025e.png