I am trying to download some seismic data in SAC format. I have used the following python code for the download.
# Loop over the events
for event in events:
# Get event information
event_time = event.origins[0].time
event_latitude = event.origins[0].latitude
event_longitude = event.origins[0].longitude
event_magnitude = event.magnitudes[0].mag
event_time_str = event_time.strftime('%Y%m%d%H%M%S')
# Loop over available networks
for network in event_client.get_stations(network=network_code, starttime=start_time, endtime=end_time, minlatitude=min_lat, maxlatitude=max_lat, minlongitude=min_lon, maxlongitude=max_lon):
for station in network:
station_code = station.code
station_latitude = station.latitude
station_longitude = station.longitude
# Calculate the back azimuth, azimuth, distance, delta, and gcarc
distance, azimuth, back_azimuth = gps2dist_azimuth(
event_latitude, event_longitude, station_latitude, station_longitude
)
distance = distance/1000
gcarc = locations2degrees(
event_latitude, event_longitude, station_latitude, station_longitude
)
if not (3330 <= distance <= 11132):
print(f"Distance between {station_code} and event {event_time} not in required range, skipping...")
continue
try:
# Download waveform data for the station
st = waveform_client.get_waveforms(
network=network_code,
station=station_code,
location="*",
channel=channel_code,
starttime=event_time - 300, # Adjust for the desired time window
endtime=event_time + 1500, # Adjust for the desired time window
attach_response=True,
)
except obspy.clients.fdsn.header.FDSNNoDataException:
print(f"No data available on {event_time} for station {station_code}, skipping...")
continue # Skip this station and move on to the next one
# Check if the Stream is empty or has no traces
if not st or len(st) == 0:
print(f"No valid data for station {station_code}, skipping...")
continue # Skip this station and move on to the next one
try:
# Attempt to remove the response
st.remove_response(output="VEL")
except Exception as e:
# Handle the exception (e.g., skip the trace)
print(f"Skipping trace for station {station_code}: {e}")
continue
# Create a directory for each event
event_directory = os.path.join(output_directory, event_time_str)
os.makedirs(event_directory, exist_ok=True)
print(f"Data available on {event_time} for station {station_code}, proceeding...")
print(st)
st.detrend("demean")
st.taper(max_percentage=0.05, type='hann')
st.filter('bandpass', freqmin=0.005, freqmax=1)
st.remove_response(output="VEL")
# Save the processed data as SAC files in the event-specific folder
for trace in st:
network = trace.stats.network
station = trace.stats.station
channel = trace.stats.channel
location = trace.stats.location
sample_rate = trace.stats.sampling_rate
sac_trace = SACTrace.from_obspy_trace(trace)
sac_trace.az = azimuth
sac_trace.baz = back_azimuth
sac_trace.dist = distance
sac_trace.delta = 1/sample_rate
sac_trace.gcarc = gcarc
sac_trace.mag = event_magnitude # Set event magnitude
sac_trace.o = 300
sac_trace.stla = station_latitude # Set station latitude
sac_trace.stlo = station_longitude # Set station longitude
sac_trace.evla = event_latitude # Set event latitude
sac_trace.evlo = event_longitude # Set event longitude
sac_file_path = os.path.join(
event_directory,
f"{event_time_str}_{network}.{station}.{location}.{channel}.SAC"
)
sac_trace.write(sac_file_path, byteorder="little")
print(f"Data on {event_time} for station {station_code} saved, proceeding...")
# Print a message indicating the folder where all the data is saved
print(f"All data saved in the folder: {output_directory}")
This is just the main part of the code. I understand that removing response in VEL results in m/s as the units. But the values are in the 1e-16 range. Is this correct? This data is then used in some program do_mft to extract amplitude vs time period data. I am self-taught in Python and ObsPy. Any and all help is appreciated. I have attached an example plot below.
