Lab 1: DAS Basics¶

In [ ]:

# Plotting paramaters
import matplotlib
import matplotlib.pyplot as plt
fontsize = 18
params = {
    'image.interpolation': 'nearest',
    'image.cmap': 'gray',
    'savefig.dpi': 300,  # to adjust notebook inline plot size
    'axes.labelsize': fontsize, # fontsize for x and y labels (was 10)
    'axes.titlesize': fontsize,
    'font.size': fontsize,
    'legend.fontsize': 18,
    'xtick.labelsize': fontsize,
    'ytick.labelsize': fontsize,
    "font.weight": 'bold',
    'text.usetex': False,
    "ytick.major.pad": 6.0,
    "xtick.major.pad": 6.0,
    "axes.labelweight": 'bold'
}
matplotlib.rcParams.update(params)
from mpl_toolkits.axes_grid1.inset_locator import inset_axes

# For timing
import datetime
import dateutil.parser
time_format = "%Y-%m-%dT%H:%M:%S.%f+00:00"

# Plotting paramaters import matplotlib import matplotlib.pyplot as plt fontsize = 18 params = { 'image.interpolation': 'nearest', 'image.cmap': 'gray', 'savefig.dpi': 300, # to adjust notebook inline plot size 'axes.labelsize': fontsize, # fontsize for x and y labels (was 10) 'axes.titlesize': fontsize, 'font.size': fontsize, 'legend.fontsize': 18, 'xtick.labelsize': fontsize, 'ytick.labelsize': fontsize, "font.weight": 'bold', 'text.usetex': False, "ytick.major.pad": 6.0, "xtick.major.pad": 6.0, "axes.labelweight": 'bold' } matplotlib.rcParams.update(params) from mpl_toolkits.axes_grid1.inset_locator import inset_axes # For timing import datetime import dateutil.parser time_format = "%Y-%m-%dT%H:%M:%S.%f+00:00"

In [ ]:

# Other necessary modules
import pandas as pd
import numpy as np

# Other necessary modules import pandas as pd import numpy as np

In [ ]:

# Adding DASutils location
pyDAS_path = "../Scripts/DAS-proc/build/"
import os
os.environ['LD_LIBRARY_PATH'] = pyDAS_path
import sys
sys.path.insert(0,pyDAS_path)
sys.path.insert(0,'../Scripts/DAS-proc/python')
import DASutils

sys.path.insert(0,'../Scripts/')
from MapTools import *

# Adding DASutils location pyDAS_path = "../Scripts/DAS-proc/build/" import os os.environ['LD_LIBRARY_PATH'] = pyDAS_path import sys sys.path.insert(0,pyDAS_path) sys.path.insert(0,'../Scripts/DAS-proc/python') import DASutils sys.path.insert(0,'../Scripts/') from MapTools import *

In [ ]:

!gdown --folder "https://drive.google.com/drive/folders/14OFsSMAAiQQ3B4jzutzDPhnKL-8B9eDs?usp=sharing" -O ../Inputs

!gdown --folder "https://drive.google.com/drive/folders/14OFsSMAAiQQ3B4jzutzDPhnKL-8B9eDs?usp=sharing" -O ../Inputs

In [ ]:

# Fiber channel locations
fiberLoc = "../Inputs/Ridgecrest-DAS-100km_South.csv"

# Fiber channel locations fiberLoc = "../Inputs/Ridgecrest-DAS-100km_South.csv"

In [ ]:

# DAS array
array_df = pd.read_csv(fiberLoc)
good_ch = array_df[array_df["status"]=="good"]["channel"].dropna().astype(int).to_numpy()
ch_lat = array_df[array_df["status"]=="good"]["latitude"].dropna().astype(float).to_numpy()
ch_lon = array_df[array_df["status"]=="good"]["longitude"].dropna().astype(float).to_numpy()
ch_z = -array_df[array_df["status"]=="good"]["elevation"].dropna().astype(float).to_numpy()*1e-3

# DAS array array_df = pd.read_csv(fiberLoc) good_ch = array_df[array_df["status"]=="good"]["channel"].dropna().astype(int).to_numpy() ch_lat = array_df[array_df["status"]=="good"]["latitude"].dropna().astype(float).to_numpy() ch_lon = array_df[array_df["status"]=="good"]["longitude"].dropna().astype(float).to_numpy() ch_z = -array_df[array_df["status"]=="good"]["elevation"].dropna().astype(float).to_numpy()*1e-3

In the first part, we will plot a local earthquake and filter the waveforms and see where the instrument noise will take over the signal. The event that we will analyze is the USGS-catalog M3.0 earthquake: link.

In [ ]:

import geopandas as gpd

minlon, maxlon, minlat, maxlat = -118.5, -116.5, 34.25, 36.25
# Fault lines from USGS catalog
faultsUSGS = gpd.GeoDataFrame.from_file("../Inputs/Qfaults_GIS_clipped/Qfaults_clipped.shp")
faultsUSGS = faultsUSGS.cx[minlon:maxlon, minlat:maxlat]

import geopandas as gpd minlon, maxlon, minlat, maxlat = -118.5, -116.5, 34.25, 36.25 # Fault lines from USGS catalog faultsUSGS = gpd.GeoDataFrame.from_file("../Inputs/Qfaults_GIS_clipped/Qfaults_clipped.shp") faultsUSGS = faultsUSGS.cx[minlon:maxlon, minlat:maxlat]

In [ ]:

ev_lat = 35.660
ev_lon = -117.538
ev_dep = 9.2 # [km]
ev_mag = 3.0
ev_time = "2026-01-08 05:28:06Z"

ev_lat = 35.660 ev_lon = -117.538 ev_dep = 9.2 # [km] ev_mag = 3.0 ev_time = "2026-01-08 05:28:06Z"

In [ ]:

aspect_ratio = 110.574 / (111.320 * np.cos(35.8 / 180.0 * np.pi))
extent = [minlon, maxlon, minlat, maxlat]

fig = plt.figure(figsize=(20, 20))
ax = plt.axes(projection=cartopy.crs.PlateCarree())
ax.set_extent(extent,crs=cartopy.crs.PlateCarree())

ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.LAND, edgecolor='black')
ax.add_feature(cartopy.feature.LAKES, edgecolor='black')
ax.add_feature(cartopy.feature.BORDERS)
ax.add_feature(cartopy.feature.STATES)
xticks=np.linspace(extent[0], extent[1], 9)
yticks=np.linspace(extent[2], extent[3], 9)
xticks1=np.linspace(extent[0], extent[1], 9)
yticks1=np.linspace(extent[2], extent[3], 9)
ax.set_xticks(xticks)
ax.set_yticks(yticks)
ax.set_axis_off()

lon_formatter = LongitudeFormatter(number_format='.2f',
                                       degree_symbol='')
lat_formatter = LatitudeFormatter(number_format='.2f',
                                      degree_symbol='')
ax.gridlines(draw_labels=True, dms=False, linestyle='-',
             xlocs=xticks1, ylocs=yticks1, xformatter=lon_formatter, yformatter=lat_formatter)

add_zebra_frame(ax, lw=5, crs=ccrs.PlateCarree(), zorder=3)

scale_bar(ax, length=50,location=(0.25, 0.025))
ax.plot(ch_lon, ch_lat, "k-", lw=4.0, label="DAS channels")
ax.scatter(ev_lon, ev_lat, c="r", s=ev_mag*40.0, edgecolors= "k", label="Earthquake")
faultsUSGS.plot(ax=ax, color="gray", lw=2.0, label="USGS faults")


_ = ax.legend()

aspect_ratio = 110.574 / (111.320 * np.cos(35.8 / 180.0 * np.pi)) extent = [minlon, maxlon, minlat, maxlat] fig = plt.figure(figsize=(20, 20)) ax = plt.axes(projection=cartopy.crs.PlateCarree()) ax.set_extent(extent,crs=cartopy.crs.PlateCarree()) ax.add_feature(cartopy.feature.OCEAN) ax.add_feature(cartopy.feature.LAND, edgecolor='black') ax.add_feature(cartopy.feature.LAKES, edgecolor='black') ax.add_feature(cartopy.feature.BORDERS) ax.add_feature(cartopy.feature.STATES) xticks=np.linspace(extent[0], extent[1], 9) yticks=np.linspace(extent[2], extent[3], 9) xticks1=np.linspace(extent[0], extent[1], 9) yticks1=np.linspace(extent[2], extent[3], 9) ax.set_xticks(xticks) ax.set_yticks(yticks) ax.set_axis_off() lon_formatter = LongitudeFormatter(number_format='.2f', degree_symbol='') lat_formatter = LatitudeFormatter(number_format='.2f', degree_symbol='') ax.gridlines(draw_labels=True, dms=False, linestyle='-', xlocs=xticks1, ylocs=yticks1, xformatter=lon_formatter, yformatter=lat_formatter) add_zebra_frame(ax, lw=5, crs=ccrs.PlateCarree(), zorder=3) scale_bar(ax, length=50,location=(0.25, 0.025)) ax.plot(ch_lon, ch_lat, "k-", lw=4.0, label="DAS channels") ax.scatter(ev_lon, ev_lat, c="r", s=ev_mag*40.0, edgecolors= "k", label="Earthquake") faultsUSGS.plot(ax=ax, color="gray", lw=2.0, label="USGS faults") _ = ax.legend()

In [ ]:

filenameDAS = "../Inputs/RidgeCrest-South-2026-01-08T050503Z.h5"

filenameDAS = "../Inputs/RidgeCrest-South-2026-01-08T050503Z.h5"

In [ ]:

eqTime = dateutil.parser.parse(ev_time)

# Reading data 
startDate = eqTime - datetime.timedelta(seconds=30)
endDate = eqTime + datetime.timedelta(seconds=120)
selectedFiles = [filenameDAS]
if selectedFiles is not None:
    DAS_data, info = DASutils.readFile_HDF(selectedFiles, 0.0, 49.0, verbose=1, diff=True, detrend=False, tapering=False, filter=False, median=False,
                                           minTime=startDate, maxTime=endDate,
                                           desampling=False, nChbuffer=10000, system="OptaSense")
    nt = DAS_data.shape[1]
    fs = info['fs']
    dChEq = info['dx']
    dt = 1.0 / fs
    ot = (dateutil.parser.parse(info['begTime'].strftime(time_format)) - eqTime).total_seconds()
    tAx_sec = np.linspace(ot, ot+(nt-1)*dt, nt)

eqTime = dateutil.parser.parse(ev_time) # Reading data startDate = eqTime - datetime.timedelta(seconds=30) endDate = eqTime + datetime.timedelta(seconds=120) selectedFiles = [filenameDAS] if selectedFiles is not None: DAS_data, info = DASutils.readFile_HDF(selectedFiles, 0.0, 49.0, verbose=1, diff=True, detrend=False, tapering=False, filter=False, median=False, minTime=startDate, maxTime=endDate, desampling=False, nChbuffer=10000, system="OptaSense") nt = DAS_data.shape[1] fs = info['fs'] dChEq = info['dx'] dt = 1.0 / fs ot = (dateutil.parser.parse(info['begTime'].strftime(time_format)) - eqTime).total_seconds() tAx_sec = np.linspace(ot, ot+(nt-1)*dt, nt)

In [ ]:

min_t = -30.0
max_t = min_t + 120
pclip=96.0

DAS_data_proc = DAS_data[good_ch,:]
minCh = 0
maxCh = DAS_data_proc.shape[0]

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
cbar.set_label('strain rate [$\mu$m/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

min_t = -30.0 max_t = min_t + 120 pclip=96.0 DAS_data_proc = DAS_data[good_ch,:] minCh = 0 maxCh = DAS_data_proc.shape[0] # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) cbar.set_label('strain rate [$\mu$m/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

In [ ]:

fmin = 0.0
fmax = 10.0 # Let's try multiple frequency: 10.0, 5.0, 2.0, 1.0, 0.1 Hz
minCh = 0
maxCh = DAS_data_proc.shape[0]
min_t = 0.0
max_t = min_t + 30 # For 0.1 Hz, let's plot longer record 120 s
pclip=96.0

DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
cbar.set_label('strain rate [$\mu$m/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

fmin = 0.0 fmax = 10.0 # Let's try multiple frequency: 10.0, 5.0, 2.0, 1.0, 0.1 Hz minCh = 0 maxCh = DAS_data_proc.shape[0] min_t = 0.0 max_t = min_t + 30 # For 0.1 Hz, let's plot longer record 120 s pclip=96.0 DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True) # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) cbar.set_label('strain rate [$\mu$m/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

Let's now plot a teleseismic event. Specifically, we will plot a M6.5 occurred earlier in 2026 in Mexico: link.

In [ ]:

import obspy
from obspy.taup import TauPyModel

def compute_azim_dist(lat_c,lon_c,lat_ev,lon_ev):
    (dist_km, az, _) = obspy.geodetics.base.calc_vincenty_inverse(lat_ev,lon_ev,lat_c,lon_c) # A(event) -> B(location)
    dist_deg = obspy.geodetics.base.locations2degrees(lat_ev,lon_ev,lat_c,lon_c)
    return az, dist_deg, dist_km

# We will also compute traveltime to find the teleseismic arrival
model = TauPyModel(model="iasp91")
def compute_TT(tele_depth,dist,allPhases=False):
    if allPhases:
        arrivals = model.get_travel_times(source_depth_in_km=tele_depth, distance_in_degree=dist,
                                              receiver_depth_in_km=0.0)
    else:
        arrivals = model.get_travel_times(source_depth_in_km=tele_depth, distance_in_degree=dist,
                                              receiver_depth_in_km=0.0,phase_list=["P", "S"])
    return arrivals

import obspy from obspy.taup import TauPyModel def compute_azim_dist(lat_c,lon_c,lat_ev,lon_ev): (dist_km, az, _) = obspy.geodetics.base.calc_vincenty_inverse(lat_ev,lon_ev,lat_c,lon_c) # A(event) -> B(location) dist_deg = obspy.geodetics.base.locations2degrees(lat_ev,lon_ev,lat_c,lon_c) return az, dist_deg, dist_km # We will also compute traveltime to find the teleseismic arrival model = TauPyModel(model="iasp91") def compute_TT(tele_depth,dist,allPhases=False): if allPhases: arrivals = model.get_travel_times(source_depth_in_km=tele_depth, distance_in_degree=dist, receiver_depth_in_km=0.0) else: arrivals = model.get_travel_times(source_depth_in_km=tele_depth, distance_in_degree=dist, receiver_depth_in_km=0.0,phase_list=["P", "S"]) return arrivals

In [ ]:

ev_lat = 16.798
ev_lon = -99.394
ev_dep = 18.0 # [km]
ev_mag = 6.5
ev_time = "2026-01-02 13:58:15Z"

# Computing distance in degrees
az, dist_deg, dist_km = compute_azim_dist(np.mean(ch_lat), np.mean(ch_lon), ev_lat, ev_lon)

ev_lat = 16.798 ev_lon = -99.394 ev_dep = 18.0 # [km] ev_mag = 6.5 ev_time = "2026-01-02 13:58:15Z" # Computing distance in degrees az, dist_deg, dist_km = compute_azim_dist(np.mean(ch_lat), np.mean(ch_lon), ev_lat, ev_lon)

In [ ]:

from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(20, 8), edgecolor='w')

# Initialize Basemap
m = Basemap(projection='cyl',
             llcrnrlat=-90, urcrnrlat=90,
             llcrnrlon=-180, urcrnrlon=180,
             resolution='c')

# Draw map details
m.drawcoastlines()
m.fillcontinents()
m.drawparallels(np.arange(-90, 90, 30), labels=[1, 1, 0, 1], fontsize=18)
m.drawmeridians(np.arange(-180, 180, 30), labels=[1, 1, 0, 1], rotation=45, fontsize=18)

plt.xlabel('Longitude', labelpad=30)
plt.ylabel('Latitude', labelpad=30)

# Convert lat and lon to map projection coordinates
lons, lats = m(ev_lon, ev_lat)
# Plot points as red dots
m.scatter(lons, lats, marker='o', color='r', zorder=5, label="Teleseismic Event")

# Plotting Mammoth Lakes
# Convert lat and lon to map projection coordinates
lons_array, lats_array = m(np.mean(ch_lon), np.mean(ch_lat))
# Plot points as blue triangles
m.scatter(lons_array, lats_array, marker='^', s=100, color='b', zorder=5, label="DAS Array")

# Add legend
plt.legend(loc=1)

# Uncomment if you want to include the title
plt.title("Azimuth: %.2f [deg]\n Distance: %.2f [deg] \n Magnitude: %s" % 
          (az, dist_deg, ev_mag))

# Show the plot
plt.show()

from mpl_toolkits.basemap import Basemap fig = plt.figure(figsize=(20, 8), edgecolor='w') # Initialize Basemap m = Basemap(projection='cyl', llcrnrlat=-90, urcrnrlat=90, llcrnrlon=-180, urcrnrlon=180, resolution='c') # Draw map details m.drawcoastlines() m.fillcontinents() m.drawparallels(np.arange(-90, 90, 30), labels=[1, 1, 0, 1], fontsize=18) m.drawmeridians(np.arange(-180, 180, 30), labels=[1, 1, 0, 1], rotation=45, fontsize=18) plt.xlabel('Longitude', labelpad=30) plt.ylabel('Latitude', labelpad=30) # Convert lat and lon to map projection coordinates lons, lats = m(ev_lon, ev_lat) # Plot points as red dots m.scatter(lons, lats, marker='o', color='r', zorder=5, label="Teleseismic Event") # Plotting Mammoth Lakes # Convert lat and lon to map projection coordinates lons_array, lats_array = m(np.mean(ch_lon), np.mean(ch_lat)) # Plot points as blue triangles m.scatter(lons_array, lats_array, marker='^', s=100, color='b', zorder=5, label="DAS Array") # Add legend plt.legend(loc=1) # Uncomment if you want to include the title plt.title("Azimuth: %.2f [deg]\n Distance: %.2f [deg] \n Magnitude: %s" % (az, dist_deg, ev_mag)) # Show the plot plt.show()

In [ ]:

# Computing traveltime for p and s body waves using tau-p code for all channels without considering local structure 
ch_list = np.concatenate((np.arange(0,good_ch.shape[0],50), [good_ch.shape[0]-1]))
P_phase_time_taup = np.zeros(ch_list.shape)
S_phase_time_taup = np.zeros(ch_list.shape)
for idx, ich in enumerate(ch_list):
    _, dist_d,_ = compute_azim_dist(ch_lat[ich], ch_lon[ich], ev_lat, ev_lon)
    arrivals = compute_TT(ev_dep, dist_d, allPhases=True)
    # Getting first P and S arrival times
    arr_names =  np.array([arrival.name for arrival in arrivals])
    P_phase_time_taup[idx] = arrivals[np.where(arr_names == "P")[0][0]].time
    S_phase_time_taup[idx] = arrivals[np.where(arr_names == "S")[0][0]].time

# Computing traveltime for p and s body waves using tau-p code for all channels without considering local structure ch_list = np.concatenate((np.arange(0,good_ch.shape[0],50), [good_ch.shape[0]-1])) P_phase_time_taup = np.zeros(ch_list.shape) S_phase_time_taup = np.zeros(ch_list.shape) for idx, ich in enumerate(ch_list): _, dist_d,_ = compute_azim_dist(ch_lat[ich], ch_lon[ich], ev_lat, ev_lon) arrivals = compute_TT(ev_dep, dist_d, allPhases=True) # Getting first P and S arrival times arr_names = np.array([arrival.name for arrival in arrivals]) P_phase_time_taup[idx] = arrivals[np.where(arr_names == "P")[0][0]].time S_phase_time_taup[idx] = arrivals[np.where(arr_names == "S")[0][0]].time

In [ ]:

filenameDAS1 = "../Inputs/RidgeCrest-South-2026-01-02T133503Z.h5"
filenameDAS2 = "../Inputs/RidgeCrest-South-2026-01-02T140503Z.h5"

filenameDAS1 = "../Inputs/RidgeCrest-South-2026-01-02T133503Z.h5" filenameDAS2 = "../Inputs/RidgeCrest-South-2026-01-02T140503Z.h5"

In [ ]:

eqTime = dateutil.parser.parse(ev_time)

# Reading data
P_time = np.mean(P_phase_time_taup) - 120.0
S_time = np.mean(S_phase_time_taup) + 500.0

startDate = eqTime + datetime.timedelta(seconds=P_time)
endDate = eqTime + datetime.timedelta(seconds=S_time)
selectedFiles = [filenameDAS1, filenameDAS2]
if selectedFiles is not None:
    DAS_data, info = DASutils.readFile_HDF(selectedFiles, 0.0, 49.0, verbose=1, diff=True, detrend=False, tapering=False, filter=False, median=False,
                                           minTime=startDate, maxTime=endDate,
                                           desampling=False, nChbuffer=10000, system="OptaSense")
    nt = DAS_data.shape[1]
    fs = info['fs']
    dChEq = info['dx']
    dt = 1.0 / fs
    ot = (dateutil.parser.parse(info['begTime'].strftime(time_format)) - eqTime).total_seconds()
    tAx_sec = np.linspace(ot, ot+(nt-1)*dt, nt)

eqTime = dateutil.parser.parse(ev_time) # Reading data P_time = np.mean(P_phase_time_taup) - 120.0 S_time = np.mean(S_phase_time_taup) + 500.0 startDate = eqTime + datetime.timedelta(seconds=P_time) endDate = eqTime + datetime.timedelta(seconds=S_time) selectedFiles = [filenameDAS1, filenameDAS2] if selectedFiles is not None: DAS_data, info = DASutils.readFile_HDF(selectedFiles, 0.0, 49.0, verbose=1, diff=True, detrend=False, tapering=False, filter=False, median=False, minTime=startDate, maxTime=endDate, desampling=False, nChbuffer=10000, system="OptaSense") nt = DAS_data.shape[1] fs = info['fs'] dChEq = info['dx'] dt = 1.0 / fs ot = (dateutil.parser.parse(info['begTime'].strftime(time_format)) - eqTime).total_seconds() tAx_sec = np.linspace(ot, ot+(nt-1)*dt, nt)

In [ ]:

min_t = tAx_sec[0]
max_t = tAx_sec[-1]
pclip=96.0

DAS_data_proc = DAS_data[good_ch,:]
minCh = 0
maxCh = DAS_data_proc.shape[0]

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave")
ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave")
ax.legend()
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
cbar.set_label('strain rate [$\mu$m/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

min_t = tAx_sec[0] max_t = tAx_sec[-1] pclip=96.0 DAS_data_proc = DAS_data[good_ch,:] minCh = 0 maxCh = DAS_data_proc.shape[0] # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave") ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave") ax.legend() posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) cbar.set_label('strain rate [$\mu$m/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

In [ ]:

fmin = 0.5
fmax = 2.0
min_t = tAx_sec[0]
max_t = tAx_sec[-1]
pclip=90.0

DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3
minCh = 0
maxCh = DAS_data_proc.shape[0]

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave")
ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave")
ax.legend()
cbar.set_label('strain rate [nm/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

fmin = 0.5 fmax = 2.0 min_t = tAx_sec[0] max_t = tAx_sec[-1] pclip=90.0 DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3 minCh = 0 maxCh = DAS_data_proc.shape[0] # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave") ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave") ax.legend() cbar.set_label('strain rate [nm/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

In [ ]:

fmin = 0.5
fmax = 2.0

min_t = np.mean(P_phase_time_taup) - 30.0
max_t = np.mean(P_phase_time_taup) + 60.0
pclip = 90.0

DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3
minCh = 0
maxCh = DAS_data_proc.shape[0]

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave")
ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave")
ax.legend()
cbar.set_label('strain rate [nm/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

fmin = 0.5 fmax = 2.0 min_t = np.mean(P_phase_time_taup) - 30.0 max_t = np.mean(P_phase_time_taup) + 60.0 pclip = 90.0 DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3 minCh = 0 maxCh = DAS_data_proc.shape[0] # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave") ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave") ax.legend() cbar.set_label('strain rate [nm/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

In [ ]:

fmin = 2.0
fmax = 5.0
minCh = 0
maxCh = DAS_data_proc.shape[0]
min_t = np.mean(S_phase_time_taup) - 30.0
max_t = np.mean(S_phase_time_taup) + 60.0
pclip = 85.0

DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)
clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip)
im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, 
           extent=[minCh, maxCh, max_t, min_t],
           aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic'))
ax.set_ylabel("Time from earthquake [s]")
ax.set_xlabel("Channel number")
posn = ax.get_position()
cbar_ax = fig.add_axes([0, 0, 0.1, 0.1])
cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height])
cbar = plt.colorbar(im, cax=cbar_ax)
ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave")
ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave")
ax.legend()
cbar.set_label('strain rate [nm/m/s]')
ax.set_xlim(minCh,maxCh)
ax.set_ylim(max_t,min_t)
ax.grid()

fmin = 2.0 fmax = 5.0 minCh = 0 maxCh = DAS_data_proc.shape[0] min_t = np.mean(S_phase_time_taup) - 30.0 max_t = np.mean(S_phase_time_taup) + 60.0 pclip = 85.0 DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch,:], fmin, fmax, dt, order=14, zerophase=True)*1e3 # Plotting data fig, ax = plt.subplots(figsize=(20,12)) it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) clipVal = np.percentile(np.absolute(DAS_data_proc[minCh:maxCh,it_min:it_max]), pclip) im = ax.imshow(DAS_data_proc[minCh:maxCh,it_min:it_max].T, extent=[minCh, maxCh, max_t, min_t], aspect='auto', vmin=-clipVal, vmax=clipVal, cmap=plt.get_cmap('seismic')) ax.set_ylabel("Time from earthquake [s]") ax.set_xlabel("Channel number") posn = ax.get_position() cbar_ax = fig.add_axes([0, 0, 0.1, 0.1]) cbar_ax.set_position([posn.x0 + posn.width + 0.01, posn.y0, 0.02, posn.height]) cbar = plt.colorbar(im, cax=cbar_ax) ax.plot(ch_list,P_phase_time_taup,"r--",lw=4,label="P-wave") ax.plot(ch_list,S_phase_time_taup,"c--",lw=4,label="S-wave") ax.legend() cbar.set_label('strain rate [nm/m/s]') ax.set_xlim(minCh,maxCh) ax.set_ylim(max_t,min_t) ax.grid()

In [ ]:

# plotting parameters
fmin = 0.5
fmax = 3.0
ich = 3000
ch_buff = 50
minCh = ich - ch_buff
maxCh = ich + ch_buff

min_t = np.mean(P_phase_time_taup) - 60.0
max_t = np.mean(S_phase_time_taup) + 400.0
it_min = int((min_t-ot)/dt+0.5)
it_max = int((max_t-ot)/dt+0.5)

DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch[minCh:maxCh],:], fmin, fmax, dt, order=14, zerophase=True)*1e3

trace = np.mean(DAS_data_proc[:, it_min:it_max], axis=0)

# Plotting data
fig, ax = plt.subplots(figsize=(20,12))
ax.plot(tAx_sec[it_min:it_max], trace, "k-", lw=0.5)

# P arrival
p_time = P_phase_time_taup[ch_list==ich]
ax.axvline(p_time, color='b', linestyle='--', linewidth=2)
ax.text(p_time, ax.get_ylim()[1], 'P', color='b', fontsize=16, ha='center', va='bottom')

# S arrival
s_time = S_phase_time_taup[ch_list==ich]
ax.axvline(s_time, color='r', linestyle='--', linewidth=2)
ax.text(s_time, ax.get_ylim()[1], 'S', color='r', fontsize=16, ha='center', va='bottom')

ax.set_xlabel("Time from earthquake [s]")
ax.set_ylabel('strain rate [nm/m/s]')
ax.set_xlim(min_t,max_t)
ax.grid()

# plotting parameters fmin = 0.5 fmax = 3.0 ich = 3000 ch_buff = 50 minCh = ich - ch_buff maxCh = ich + ch_buff min_t = np.mean(P_phase_time_taup) - 60.0 max_t = np.mean(S_phase_time_taup) + 400.0 it_min = int((min_t-ot)/dt+0.5) it_max = int((max_t-ot)/dt+0.5) DAS_data_proc = DASutils.bandpass2D_c(DAS_data[good_ch[minCh:maxCh],:], fmin, fmax, dt, order=14, zerophase=True)*1e3 trace = np.mean(DAS_data_proc[:, it_min:it_max], axis=0) # Plotting data fig, ax = plt.subplots(figsize=(20,12)) ax.plot(tAx_sec[it_min:it_max], trace, "k-", lw=0.5) # P arrival p_time = P_phase_time_taup[ch_list==ich] ax.axvline(p_time, color='b', linestyle='--', linewidth=2) ax.text(p_time, ax.get_ylim()[1], 'P', color='b', fontsize=16, ha='center', va='bottom') # S arrival s_time = S_phase_time_taup[ch_list==ich] ax.axvline(s_time, color='r', linestyle='--', linewidth=2) ax.text(s_time, ax.get_ylim()[1], 'S', color='r', fontsize=16, ha='center', va='bottom') ax.set_xlabel("Time from earthquake [s]") ax.set_ylabel('strain rate [nm/m/s]') ax.set_xlim(min_t,max_t) ax.grid()

In [ ]:

# ------------------- Parameters -------------------
ich = 3000
ch_buff = 50
minCh = ich - ch_buff
maxCh = ich + ch_buff
min_t = np.mean(P_phase_time_taup) - 60.0
max_t = np.mean(S_phase_time_taup) + 400.0
it_min = int((min_t - ot) / dt + 0.5)
it_max = int((max_t - ot) / dt + 0.5)

p_time = P_phase_time_taup[ch_list == ich]
s_time = S_phase_time_taup[ch_list == ich]

# ------------------- Frequency bands -------------------
# For teleseismic: body waves 0.5-2 Hz, surface waves 0.01-0.1 Hz
freq_bands = [
    (0.01, 0.1,  "0.01–0.1 Hz\nSurface waves"),
    (0.1,  0.5,  "0.1–0.5 Hz\nLong period body"),
    (0.5,  2.0,  "0.5–2.0 Hz\nBody waves"),
    (1.0,  3.0,  "1.0–3.0 Hz\nShort period body"),
    (2.0,  8.0,  "2.0–8.0 Hz\nHigh freq / coda"),
]

# ------------------- Compute and normalize traces -------------------
traces = []
for fmin, fmax, label in freq_bands:
    DAS_proc = DASutils.bandpass2D_c(DAS_data[good_ch[minCh:maxCh], :], fmin, fmax, dt, order=14, zerophase=True) * 1e3
    tr = np.mean(DAS_proc[:, it_min:it_max], axis=0)
    tr /= (np.max(np.abs(tr)) + 1e-10)  # normalize to [-1, 1]
    traces.append((tr, label))

# ------------------- Plot -------------------
spacing = 2.5  # vertical separation between traces

fig, ax = plt.subplots(figsize=(20, 12))

for i, (tr, label) in enumerate(traces):
    offset = i * spacing
    ax.plot(tAx_sec[it_min:it_max], tr + offset, 'k', linewidth=0.8)
    ax.text(min_t+3.0, offset+0.5, label.split('\n')[-1], fontsize=14, va='bottom', ha="left",
            bbox=dict(fc='white', ec='none', alpha=0.7))

# P arrival
ax.axvline(p_time, color='b', linestyle='--', linewidth=2)
ax.text(p_time, 1.02, 'P', color='b', fontsize=16, ha='center', va='bottom',
        transform=ax.get_xaxis_transform(), clip_on=False)

# S arrival
ax.axvline(s_time, color='r', linestyle='--', linewidth=2)
ax.text(s_time, 1.02, 'S', color='r', fontsize=16, ha='center', va='bottom',
        transform=ax.get_xaxis_transform(), clip_on=False)

ax.set_xlabel("Time from earthquake [s]", fontsize=16)
ax.set_ylabel("Frequency band", fontsize=16)
ax.set_xlim(min_t, max_t)
ax.set_yticks([i * spacing for i in range(len(traces))])
ax.set_yticklabels([label.split('\n')[0] for _, label in traces], fontsize=12)
ax.grid(axis='x')
plt.tight_layout()
plt.show()

# ------------------- Parameters ------------------- ich = 3000 ch_buff = 50 minCh = ich - ch_buff maxCh = ich + ch_buff min_t = np.mean(P_phase_time_taup) - 60.0 max_t = np.mean(S_phase_time_taup) + 400.0 it_min = int((min_t - ot) / dt + 0.5) it_max = int((max_t - ot) / dt + 0.5) p_time = P_phase_time_taup[ch_list == ich] s_time = S_phase_time_taup[ch_list == ich] # ------------------- Frequency bands ------------------- # For teleseismic: body waves 0.5-2 Hz, surface waves 0.01-0.1 Hz freq_bands = [ (0.01, 0.1, "0.01–0.1 Hz\nSurface waves"), (0.1, 0.5, "0.1–0.5 Hz\nLong period body"), (0.5, 2.0, "0.5–2.0 Hz\nBody waves"), (1.0, 3.0, "1.0–3.0 Hz\nShort period body"), (2.0, 8.0, "2.0–8.0 Hz\nHigh freq / coda"), ] # ------------------- Compute and normalize traces ------------------- traces = [] for fmin, fmax, label in freq_bands: DAS_proc = DASutils.bandpass2D_c(DAS_data[good_ch[minCh:maxCh], :], fmin, fmax, dt, order=14, zerophase=True) * 1e3 tr = np.mean(DAS_proc[:, it_min:it_max], axis=0) tr /= (np.max(np.abs(tr)) + 1e-10) # normalize to [-1, 1] traces.append((tr, label)) # ------------------- Plot ------------------- spacing = 2.5 # vertical separation between traces fig, ax = plt.subplots(figsize=(20, 12)) for i, (tr, label) in enumerate(traces): offset = i * spacing ax.plot(tAx_sec[it_min:it_max], tr + offset, 'k', linewidth=0.8) ax.text(min_t+3.0, offset+0.5, label.split('\n')[-1], fontsize=14, va='bottom', ha="left", bbox=dict(fc='white', ec='none', alpha=0.7)) # P arrival ax.axvline(p_time, color='b', linestyle='--', linewidth=2) ax.text(p_time, 1.02, 'P', color='b', fontsize=16, ha='center', va='bottom', transform=ax.get_xaxis_transform(), clip_on=False) # S arrival ax.axvline(s_time, color='r', linestyle='--', linewidth=2) ax.text(s_time, 1.02, 'S', color='r', fontsize=16, ha='center', va='bottom', transform=ax.get_xaxis_transform(), clip_on=False) ax.set_xlabel("Time from earthquake [s]", fontsize=16) ax.set_ylabel("Frequency band", fontsize=16) ax.set_xlim(min_t, max_t) ax.set_yticks([i * spacing for i in range(len(traces))]) ax.set_yticklabels([label.split('\n')[0] for _, label in traces], fontsize=12) ax.grid(axis='x') plt.tight_layout() plt.show()