Dear Stormchaser,
Thank you for your assistance previously. Following your guidance, I have made modifications to my data processing code, but I’ve encountered some issues, particularly with the PGV and PGD calculations. Specifically, I’ve noticed that the PGV values seem to be lower than expected, while the PGD values are unusually high. This has led me to suspect that there might be some problems in my data processing workflow, especially in the step where the instrument response is removed.
I have tried to check and debug my code to the best of my ability, yet the issue persists. I am reaching out for further advice and help from you. Below is an overview of my code location, output examples, and the response removal file I’ve used:
import socket
from obspy import Stream, Trace, UTCDateTime
from obspy.clients.fdsn import Client
from obspy.core.inventory import read_inventory
import numpy as np
# UDP socket configuration
UDP_IP = "0.0.0.0"
UDP_PORT = 8888
BUFFER_SIZE = 1024
# Function to parse UDP data packet
def parse_udp_packet(data):
# Parse data packet
channel, timestamp, *counts = data.strip('{}').split(',')
# Convert timestamp to float
timestamp = float(timestamp)
# Convert counts to integers
counts = list(map(int, counts))
return channel, timestamp, counts
# Create a socket and bind to address
sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
sock.bind((UDP_IP, UDP_PORT))
# Load instrument response from local XML file (SeisComP format)
inv = read_inventory(r".\R0EB1.xml")
# Initialize ObsPy Stream
stream = Stream()
# Conversion factor: counts to micrometers
sensitivity = 1 # Sensitivity for RS1D(is wrong)
conversion_factor_acceleration = sensitivity * 1000000 # Convert from counts/g to micrometers/second^2
conversion_factor_velocity = sensitivity * 1000 # Convert from counts/g to micrometers/second
# Bandpass filter parameters
freqmin = 0.1 # Hz
freqmax = 20 # Hz
#Wait time
wait_time=0.25#(s)
max_times=wait_time*4
times=0
# Read data packets, filter and calculate PGA, PGV, and PGD
while True:
stream = Stream()
# Receive data packet
data, addr = sock.recvfrom(BUFFER_SIZE)
# Parse data packet
channel, timestamp, counts = parse_udp_packet(data.decode('utf-8'))
# print(counts)
# Create ObsPy Trace object
trace = Trace(data=np.array(counts), header={"starttime": timestamp})
trace.stats.network = "AM"
trace.stats.station = "R0EB1"
trace.stats.location = "00"
trace.stats.channel = "EHZ"
# Apply bandpass filter
trace.filter('bandpass', freqmin=freqmin, freqmax=freqmax, corners=4)
# Attach instrument response
# print(trace)
trace.attach_response(inv)
# Append Trace object to Stream
stream += trace
# Remove instrument response
stream_resp_removed = stream.remove_response(output="VEL")
# Calculate PGA, PGV, and PGD
pga_val = np.max(np.abs(stream_resp_removed[0].data)) * conversion_factor_acceleration
pgv_val = np.max(np.abs(np.diff(stream_resp_removed[0].data))) * conversion_factor_velocity
pgd_val = np.max(np.abs(np.cumsum(stream_resp_removed[0].data))) * conversion_factor_velocity
# pga_val = np.max(np.abs(stream[0].data)) * conversion_factor_acceleration
# pgv_val = np.max(np.abs(np.diff(stream[0].data))) * conversion_factor_velocity
# pgd_val = np.max(np.abs(np.cumsum(stream[0].data))) * conversion_factor_velocity
# Print PGA, PGV, and PGD
times=times+1
if times>=max_times:
print("Channel: {}, Timestamp: {}".format(channel, UTCDateTime(timestamp)))
print("PGA: {:.2f} µm/s²".format(pga_val))
print("PGV: {:.2f} µm/s".format(pgv_val))
print("PGD: {:.2f} µm".format(pgd_val))
print("-----------------------------")
times=0
Channel: 'EHZ', Timestamp: 2024-02-20T14:24:15.377000Z
PGA: 1780.53 µm/s²
PGV: 0.55 µm/s
PGD: 23.73 µm
-----------------------------
Channel: 'EHZ', Timestamp: 2024-02-20T14:24:15.627000Z
PGA: 1800.66 µm/s²
PGV: 0.55 µm/s
PGD: 23.42 µm
-----------------------------
Channel: 'EHZ', Timestamp: 2024-02-20T14:24:15.877000Z
PGA: 1784.98 µm/s²
PGV: 0.54 µm/s
PGD: 23.50 µm
-----------------------------
Channel: 'EHZ', Timestamp: 2024-02-20T14:24:16.127000Z
PGA: 1789.50 µm/s²
PGV: 0.55 µm/s
PGD: 23.51 µm
-----------------------------
Channel: 'EHZ', Timestamp: 2024-02-20T14:24:16.377000Z
PGA: 1800.50 µm/s²
PGV: 0.55 µm/s
PGD: 23.24 µm
-
PGA: The PGA values seem reasonable.
-
PGV: The results appear to be too low, far below what was expected.
-
PGD: The calculated values are abnormally high, which does not align with my understanding of seismic data.
-
Response Removal File: The file I used for removing the instrument response is
<?xml version="1.0" encoding="UTF-8"?>
<seiscomp xmlns="http://geofon.gfz-potsdam.de/ns/seiscomp3-schema/0.11" version="0.11">
<Inventory>
<sensor publicID="Sensor-1D-RS-V7-VEL" name="S-1D-RS-V7-VEL" response="ResponsePAZ-1D-RS-V7-VEL">
<description>Velocity</description>
<manufacturer>Raspberry Shake</manufacturer>
<unit>M/S</unit>
<remark>{"unit":"Velocity in Meters Per Second"}</remark>
</sensor>
<datalogger publicID="Datalogger-SHAKE-100hz" name="DL-RS-100">
<gain>1</gain>
<maxClockDrift>0</maxClockDrift>
<decimation sampleRateNumerator="100" sampleRateDenominator="1">
<digitalFilterChain>ResponseFIR-DL100</digitalFilterChain>
</decimation>
</datalogger>
<responsePAZ publicID="ResponsePAZ-1D-RS-V7-VEL" name="RP-1D-RS-V7-VEL">
<type>A</type>
<gain>399650000</gain>
<gainFrequency>5</gainFrequency>
<normalizationFactor>673.744</normalizationFactor>
<normalizationFrequency>5</normalizationFrequency>
<numberOfZeros>3</numberOfZeros>
<numberOfPoles>4</numberOfPoles>
<zeros>(0,0) (0,0) (0,0)</zeros>
<poles>(-1,0) (-3.03,0) (-3.03,0) (-666.67,0)</poles>
</responsePAZ>
<responseFIR publicID="ResponseFIR-DL100" name="RF-DL100.stage_2">
<gain>1</gain>
<decimationFactor>2</decimationFactor>
<delay>0</delay>
<correction>0</correction>
<numberOfCoefficients>38</numberOfCoefficients>
<symmetry>B</symmetry>
<coefficients>-0.000151873 0.000602538 -9.55046e-05 -0.000672574 -8.57924e-05 0.00101667 0.000443862 -0.00137953 -0.0010729 0.00164772 0.0020288 -0.00165744 -0.00330916 0.00120522 0.00482924 -7.07612e-05 -0.00640028 -0.00194821 0.0077177 0.0049948 -0.00836173 -0.00911258 0.00780703 0.0142144 -0.00542901 -0.0200678 0.000475993 0.0263035 0.00807014 -0.0324465 -0.0219791 0.0379693 0.0455535 -0.0423585 -0.0946718 0.0451835 0.314401 0.453841</coefficients>
</responseFIR>
<network publicID="AM.Network" code="AM">
<start>2024-02-08T00:00:00.00Z</start>
<description>RASPBERRY SHAKE NETWORK</description>
<netClass>p</netClass>
<archive>CAPS</archive>
<restricted>false</restricted>
<shared>true</shared>
<station publicID="R0EB1.Station" code="R0EB1">
<start>2024-02-08T00:00:00.00Z</start>
<description>Raspberry Shake Citizen Science Station</description>
<latitude>29.8413</latitude>
<longitude>106.0663</longitude>
<elevation>300</elevation>
<affiliation>RaspShake Network</affiliation>
<restricted>false</restricted>
<shared>true</shared>
<sensorLocation publicID="R0EB1.00.Location" code="00">
<start>2024-02-08T00:00:00.00Z</start>
<latitude>29.8413</latitude>
<longitude>106.0663</longitude>
<elevation>300</elevation>
<stream publicID="Stream/1D-EHZ" code="EHZ" datalogger="Datalogger-SHAKE-100hz" sensor="Sensor-1D-RS-V7-VEL">
<start>2024-02-08T00:00:00.00Z</start>
<dataloggerChannel>0</dataloggerChannel>
<sensorChannel>0</sensorChannel>
<sampleRateNumerator>100</sampleRateNumerator>
<sampleRateDenominator>1</sampleRateDenominator>
<depth>0</depth>
<azimuth>0</azimuth>
<dip>-90</dip>
<gain>399650000</gain>
<gainFrequency>5</gainFrequency>
<gainUnit>M/S</gainUnit>
<format>Steim2</format>
<flags>CG</flags>
<restricted>false</restricted>
<shared>true</shared>
</stream>
</sensorLocation>
</station>
</network>
</Inventory>
</seiscomp>
Given these circumstances, I suspect there might be a mistake in how I’m handling the instrument response or in the data preprocessing stage. I would greatly appreciate your insight on what I might be overlooking or misunderstanding. Any advice on how to improve my code or data processing workflow would be invaluable.
Thank you once again for your help and support.
Best regards,
