This page is a static rendering of the Jupyter notebook 3_tutorial_noisepy_scedc_s3_explained.ipynb. The raw notebook is available at https://github.com/cloudbank-project/cloudbank_showcase/blob/main/cloud_analysis_seismology/3_tutorial_noisepy_scedc_s3_explained.ipynb.
SCEDC Tutorial: Cloud-Powered Seismic Noise Processing with NoisePy on AWS S3
This tutorial demonstrates how to leverage NoisePy, a powerful Python package for ambient seismic noise cross-correlations, within a cloud environment using AWS S3. We will process a toy dataset from the SCEDC (Southern California Earthquake Data Center) directly on S3, showcasing how cloud storage and computational resources enable scalable and efficient seismic data analysis.
The primary goal is to illustrate the workflow for accessing, processing, and storing large seismic datasets entirely in the cloud, minimizing local data transfers and leveraging AWS's robust infrastructure.
The raw seismic waveform data and station metadata are conveniently hosted on AWS S3 as part of the SCEDC Public Data Set: https://scedc.caltech.edu/data/getstarted-pds.html
Why Cloud for Seismic Data Processing?
Processing large volumes of seismic data traditionally requires significant local storage and computational power. Cloud platforms like AWS offer several advantages:
- Scalability: Easily scale compute resources up or down based on your processing needs, without managing physical hardware.
- Accessibility: Access data and run analyses from anywhere with an internet connection.
- Cost-Effectiveness: Pay only for the resources you consume, avoiding large upfront investments.
- Data Proximity: Process data closer to where it resides (e.g., directly on S3), reducing data transfer times and costs.
- Collaboration: Share data and workflows more easily with collaborators.
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
# Ensure datetimerange is available for DateTimeRange usage
import importlib.util
if importlib.util.find_spec("datetimerange") is None:
import subprocess
import sys
subprocess.check_call([sys.executable, "-m", "pip", "install", "datetimerange"])
Understanding NoisePy: A Brief Overview
NoisePy is a specialized Python library designed for processing continuous seismic noise data to extract valuable information about the Earth's subsurface. It primarily focuses on:
- Cross-correlation: This technique compares seismic noise signals recorded at different locations to derive a "Green's function," which is essentially a synthetic seismogram representing the Earth's response between those two points. Think of it like listening to the echoes between two points to understand the path the sound traveled.
- Stacking: After calculating many cross-correlations over time, these are stacked (averaged) to enhance coherent signals and suppress random noise, revealing clearer Earth properties.
By performing these operations, seismologists can image the Earth's interior, monitor changes over time, and study seismic wave propagation without relying on large, infrequent earthquakes.
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
Here we import the various Python modules required for our seismic data processing workflow. Special attention is given to modules that facilitate interaction with AWS S3, enabling us to work with cloud-based data efficiently:
noisepy.seisfunctions: Core NoisePy routines for cross-correlation and stacking.NumpyCCStore,NumpyStackStore: These are storage interfaces, in this case, for saving processed correlation and stacked data.SCEDCS3DataStore: Crucially, this module allows NoisePy to directly read seismic waveform data from the SCEDC Public Dataset on AWS S3, treating it as a local filesystem without requiring explicit downloads.XMLStationChannelCatalog: This module handles the station metadata (like sensor locations and instrument responses) which is also sourced from an S3 bucket.ConfigParameters: The central object to define and manage all processing parameters.boto3: The Amazon Web Services (AWS) SDK for Python. We use this to securely retrieve AWS credentials configured in your environment, ensuring that we can write data back to your S3 buckets without hardcoding sensitive information.
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
from noisepy.seis.io.numpystore import NumpyCCStore, NumpyStackStore
from noisepy.seis.io.s3store import SCEDCS3DataStore # Object to query SCEDC data from on S3
from noisepy.seis.io.channel_filter_store import channel_filter
from noisepy.seis.io.datatypes import CCMethod, ConfigParameters, FreqNorm, RmResp, StackMethod, TimeNorm # Main configuration object
from noisepy.seis.io.channelcatalog import XMLStationChannelCatalog # Required stationXML handling object
import os
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 = "./scedc_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}}
We will work with a single day worth of data on SCEDC. The continuous data is organized with a single day and channel per miniseed (https://scedc.caltech.edu/data/cloud.html). 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.
# S3_DATA points to the public AWS S3 bucket containing the continuous waveform data from SCEDC.
# NoisePy is configured to read directly from this remote location, leveraging cloud data access.
S3_DATA = "s3://scedc-pds/continuous_waveforms/"
# Define the timeframe for our analysis. This tutorial focuses on a single day of data.
# The `datetime` objects specify the start and end of this period in UTC.
start = datetime(2002, 1, 1, tzinfo=timezone.utc)
end = datetime(2002, 1, 2, tzinfo=timezone.utc)
timerange = DateTimeRange(start, end)
print(timerange)
The station information, including the instrumental response, is crucial for accurate seismic data processing. This metadata is stored as StationXML files in a publicly accessible AWS S3 bucket. NoisePy can efficiently access this information directly from S3, eliminating the need for local downloads and streamlining the workflow in a cloud environment.
S3_STATION_XML = "s3://scedc-pds/FDSNstationXML/CI/" # This specifies the S3 bucket path where StationXML files (metadata for seismic stations) are located. NoisePy reads this directly from AWS S3.
Ambient Noise Project Configuration
The ConfigParameters() object in NoisePy is central to defining how our seismic data processing workflow will operate. It holds all the critical settings for data acquisition, preprocessing, cross-correlation, and stacking.
When working in a cloud environment, careful consideration of these parameters can significantly impact performance and cost. For instance:
inc_hours: This parameter defines the chunk of time for which data will be loaded and processed in parallel. In a cloud context, a largerinc_hoursmight mean more data processed per task, potentially reducing the overhead of starting new tasks. However, it also means more memory usage per worker. Optimizing this value is crucial for balancing speed and resource utilization on cloud instances.sampling_rate: Defines the target sampling rate. Higher rates mean more data points and thus higher computational demand.cc_len,step: These control the length and overlap of the time windows for cross-correlation, directly influencing the number of computations.
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 = True # StationXML metadata is required to remove the instrument response
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)
# 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
A core advantage of NoisePy's cloud integration is its ability to directly read and process data residing in remote S3 buckets. In this step, we configure the data store to pull seismic waveforms directly from the SCEDC Public Dataset on AWS S3.
This approach means we are not downloading gigabytes or terabytes of raw data to our local machine or compute instance before processing. Instead, NoisePy intelligently streams or accesses the necessary data chunks from S3 as needed, significantly reducing I/O bottlenecks, storage requirements, and transfer costs in a cloud-native workflow.
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.
Declaring Data and Cross-Correlation Stores
Here, we initialize key components for our seismic data processing workflow, showcasing both cloud-native data access and local storage for intermediate results:
catalog = XMLStationChannelCatalog(...): This object is responsible for managing station metadata (StationXML files). By pointing it toS3_STATION_XML, NoisePy knows where to find essential information about each seismic station, including its location and instrument response characteristics. Thestorage_options=S3_STORAGE_OPTIONSensures it uses anonymous access for public S3 data, a common pattern for public datasets on AWS.raw_store = SCEDCS3DataStore(...): This is our primary interface for reading raw seismic waveform data directly from theS3_DATAbucket. It's configured with thecatalogfor metadata and achannel_filterto select specific networks, stations, and channels within our definedtimerange. This demonstrates a true cloud-native approach, as data is processed in place without prior local download.cc_store = ASDFCCStore(cc_data_path): Whileraw_storereads from S3, the intermediate cross-correlation results are stored locally in the directory specified bycc_data_path. This is a common practice to optimize performance during the computationally intensive cross-correlation phase, writing results to fast local storage (e.g., an attached SSD on a cloud instance) before potentially transferring final products back to S3 for long-term storage or sharing.
# Define the seismic networks, stations, and channels we are interested in.
# This acts as a filter for the vast amount of data available on S3.
config.networks = ["CI"]
config.stations = ["RPV", "SVD", "BBR"]
config.channels = ["BH?"]
# Initialize the station catalog: NoisePy will fetch station metadata (StationXML) from the specified S3 bucket.
catalog = XMLStationChannelCatalog(S3_STATION_XML, storage_options=S3_STORAGE_OPTIONS) # Reads metadata from S3
# Initialize the raw data store: This is where NoisePy will read continuous seismic waveforms from.
# It's configured to directly access the SCEDC public S3 bucket using the defined filters and time range.
raw_store = SCEDCS3DataStore(S3_DATA, catalog,
channel_filter(config.networks, config.stations, config.channels),
timerange, storage_options=S3_STORAGE_OPTIONS) # Reads raw data directly from S3
get the time range of the data in the data store inventory
span = raw_store.get_timespans()
print(span)
Get the channels available during a given time spane
channel_list=raw_store.get_channels(span[0])
print(channel_list)
Perform the cross correlation
The data will be pulled from SCEDC, cross correlated, and stored locally if this notebook is ran locally.
If you are re-calculating, we recommend to clear the old cc_store.
from noisepy.seis.io.asdfstore import ASDFCCStore, ASDFStackStore # Object to store ASDF data within noisepy
# The `cc_store` is set up to store the intermediate cross-correlation results.
# For performance, these are often written to local storage (e.g., attached disk on a cloud instance)
# before eventual transfer of final products to a more permanent cloud storage like S3.
cc_store = ASDFCCStore(cc_data_path)
cross_correlate(raw_store, config, cc_store)
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
After computing individual cross-correlations, the next crucial step is to stack them. Stacking involves averaging multiple cross-correlation functions, which significantly enhances the coherent signals (representing Earth's subsurface properties) and suppresses random noise. This process improves the signal-to-noise ratio, making the subtle seismic signals clearer and more interpretable.
In a cloud workflow, storing the stacked results back to AWS S3 offers several advantages: persistent storage, easy sharing with collaborators, and accessibility for further analysis from any compute resource. Here, we configure the stack stores to write our final products directly to a designated S3 bucket. You must create and own this bucket (or a specific prefix) before continuing and supply its path via the STACK_STORE_PATH environment variable or by editing the notebook.
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.
# Define the S3 bucket path where the stacked cross-correlation results will be saved.
# IMPORTANT: Provide the path to an existing bucket/prefix you control. You can set the
# STACK_STORE_PATH environment variable (e.g., STACK_STORE_PATH=s3://my-bucket/noisepy-stacks)
# before launching Jupyter, or edit the fallback value below.
import os
STACK_STORE_PATH = os.environ.get("STACK_STORE_PATH", "s3://stack-store-bucket-demo").strip()
if STACK_STORE_PATH == "s3://stack-store-bucket-demo":
raise ValueError(
"Set the STACK_STORE_PATH environment variable or edit this notebook to point to an existing S3 bucket "
"that you own. The bucket must be created ahead of time."
)
s3_path = STACK_STORE_PATH
import boto3
# Use boto3 to automatically fetch AWS credentials from your environment (e.g., ~/.aws/credentials, IAM role).
# This avoids hardcoding sensitive access keys directly in the notebook.
session = boto3.Session()
credentials = session.get_credentials()
current_credentials = credentials.get_frozen_credentials()
# Configure storage options for writing to S3, using the retrieved credentials.
# This allows NoisePy to authenticate and write data to your specified S3 bucket.
S3_STORAGE_WRITE_OPTIONS = {"s3": {"key": current_credentials.access_key,
"secret": current_credentials.secret_key,
"token": current_credentials.token}} # Include session token for temporary credentials if applicable
# Open the cross-correlation store in read-only mode, as we will now only be reading these intermediate results for stacking.
cc_store = ASDFCCStore(cc_data_path, mode="r")
# Initialize the stack store. This is configured to write the final stacked results to the S3 bucket
# specified by `s3_path`, using the AWS credentials configured in `S3_STORAGE_WRITE_OPTIONS`.
# This ensures our final processed data is stored persistently and accessibly in the cloud.
stack_store = NumpyStackStore(s3_path, storage_options=S3_STORAGE_WRITE_OPTIONS)
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)
cc_store.get_station_pairs()
stack_cross_correlations(cc_store, stack_store, config)
QC_1: Quality Control of Cross-Correlations for Imaging
After the intensive cross-correlation and stacking processes, it's crucial to perform Quality Control (QC) to ensure the results are meaningful and suitable for further analysis, such as seismic imaging. This section focuses on visualizing the stacked cross-correlation functions (CCFs) to assess their quality.
Understanding the Moveout Plot
We generate a moveout plot (sometimes called a record section). This type of plot displays the stacked cross-correlation functions (CCFs) for various station pairs, ordered by their inter-station distance (distance between the two sensors). It helps us visualize how seismic waves propagate across different distances.
What you see in the plot:
- X-axis (Time Lag): Represents the time difference between signals arriving at the two stations. Positive lags correspond to waves traveling from the source to the receiver, while negative lags correspond to waves traveling from the receiver to the source.
- Y-axis (Inter-station Distance): The physical separation between the pairs of seismic stations.
- Waveforms (Color/Amplitude): The actual stacked cross-correlation functions. Strong, coherent arrivals indicate robustly correlated noise signals, suggesting clear wave propagation paths.
Interpretation of the '3 Lines' (or main features) in a typical plot:
- Positive Lag Coherent Arrival: This prominent feature typically represents the causal part of the Green's function, corresponding to seismic energy traveling from one station to the other. Its slope (change in arrival time with distance) gives insight into the average velocity of seismic waves in the subsurface.
- Negative Lag Coherent Arrival: This feature is the anti-causal part, representing energy traveling in the opposite direction. Ideally, in a perfectly symmetric noise field, it would be a mirror image of the positive lag. Its presence and symmetry are important indicators of data quality and noise source distribution.
- Zero Lag Energy: Often, there's a burst of energy around zero lag. This can correspond to scattering from near-surface heterogeneities or instrument-related noise. Ideally, for imaging purposes, we primarily focus on the coherent arrivals at non-zero lags.
By filtering the cross-correlations in a specific frequency band (e.g., between 0.1 Hz and 0.2 Hz, as configured), we can highlight different seismic phases and investigate frequency-dependent Earth properties. A clear, symmetric moveout pattern with distinct coherent arrivals indicates high-quality cross-correlation results.
cc_store.get_station_pairs()
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
plot_all_moveout(sta_stacks, 'Allstack_linear', 0.1, 0.2, 'ZZ', 1)