I am trying to plot spectra (using PPSD) for hydrophone data. Which variables are obligatory for the metadata in this case? How can I use an array for frequency dependent sensitivity values if ‘poles’ and ‘zeros’ are unknown?
All parameters that PPSD has during initialization of the object apply and change the results of the analysis, but you can just try out the default values for a first try and these parameters mostly change the “looks and feel” of the output data (like how wide the bins are and how much the spectra are smoothed).
If I understand correctly, you have a list of amplitude and phase values (or complex response value) for a given set of discrete frequencies (usually referred to as “frequency-amplitude-phase” pairs)? If that’s how it is, you can set up responses for your stations as described in the “StationXML from scratch” tutorial (see below) and just set up your response manually instead of using NRL by putting together each Channel’s response attribute using the related objects in obspy.core.inventory.response (Response, ResponseListResponseStage, ResponseListElement, compare this post)
The frequency response of the hydrophine is given by two arrays and then interpolated:
define hydrophone calibration data
calfreq = [0,250,10000,20100,30100,40200,50200,60200,70300,80300,90400,100400,110400,120500,130500,140500,150600,160600,170700,180700,190700,200000]
calsens = [-177.90,-177.90,-176.80,-176.35,-177.05,-177.35,-177.30,-178.05,-178.00,-178.40,-178.85,-180.25,-180.50,-179.90,-180.15,-180.20,-180.75,-180.90,-181.45,-181.30,-180.75,-180.30]
interpolate to the frequency resolution of the spectrogram
tck = interpolate.splrep(calfreq, calsens, s=0)
isens = interpolate.splev(f, tck, der=0)
plt.figure(dpi=300)
im = plt.plot(calfreq,calsens,‘bo’,f,isens,‘g’)
plt.xlabel(‘Frequency (Hz)’)
plt.ylabel(‘Hydrophone sensitivity (dB re V/uPA)’)
plt.legend([‘Factory, measured’, ‘Interpolated’])
Are the default values equal to the ones defined for paz in the PPSD doc?
paz = {‘gain’: 60077000.0,
‘poles’: [-0.037004+0.037016j, -0.037004-0.037016j,
-251.33+0j, -131.04-467.29j, -131.04+467.29j],
‘sensitivity’: 2516778400.0,
‘zeros’: [0j, 0j]}
The second part of my answer above seems to apply, but you would need a third array with the phase values at those given frequencies. Interpolation would be done internally by evalresp. No need for any poles/zeros.
Instrument Sensitivity will need value and frequency to be 1.0
import numpy as np
import obspy
from obspy.core.inventory import Inventory, Network, Station, Channel, Site
from obspy.core.inventory.response import InstrumentSensitivity,Response,ResponseListElement, ResponseListResponseStage
#%%
inv = Inventory(networks=[],source="MARS-Metadata")
net = Network(code="MARS", stations=[], description="Hydrophone 1",start_date=obspy.UTCDateTime(2017, 6, 13))
sta = Station(code="", #this needs to be the same as in the traces
latitude=36.712468,
longitude=-122.186898,
elevation=0.0,
creation_date=obspy.UTCDateTime(2016, 1, 2),
site=Site(name="First station"))
cha = Channel(code="", #this needs to be the same as in the trace
location_code="", #this needs to be the same as in the traces
latitude=36.712468,
longitude=-122.186898,
elevation=0.0,
depth=891.0,
azimuth=0.0,
dip=-90.0,
sample_rate=256e3)
#%%
# define hydrophone calibration data: deployment 2 since 13 June 2017
calfreq = [250, 10000, 20100, 30100, 40200, 50200, 60200, 70300, 80300, 90400, 100400, 110400, 120500, 130500, 140500, 150600, 160600, 170700, 180700, 190700, 200000]
calsens = [-177.90,-176.80,-176.35,-177.05,-177.35,-177.30,-178.05,-178.00,-178.40,-178.85,-180.25,-180.50,-179.90,-180.15,-180.20,-180.75,-180.90,-181.45,-181.30,-180.75,-180.30]
calfreq = np.array(calfreq)
calsens = np.array(calsens)
phase = np.zeros(calfreq.shape)
elements=list()
i=0
while i < len(calfreq):
elements.append(ResponseListElement(calfreq[i], calsens[i], phase[i]))
i+=1
responseLS = list()
responseLS.append(ResponseListResponseStage(stage_sequence_number=1, stage_gain=1, stage_gain_frequency=1, input_units="Pa", output_units="V",response_list_elements=elements))
insS=InstrumentSensitivity(value=1.0,frequency=1.0,input_units="PA", output_units="V")
response = Response(resource_id=None, instrument_sensitivity=insS, instrument_polynomial=None, response_stages=responseLS)
response.plot(min_freq=1e3, sampling_rate=200000*2)
cha.response = response
sta.channels.append(cha)
net.stations.append(sta)
inv.networks.append(net)
inv.write("stationMARS.xml", format="stationxml", validate=True)
Does ResponseListElement expect the amplitude value to be in dB or as a scalar?
calsens = [-177.90,-177.90,-176.80,-176.35,-177.05,-177.35,-177.30,-178.05,-178.00,-178.40,-178.85,-180.25,-180.50,-179.90,-180.15,-180.20,-180.75,-180.90,-181.45,-181.30,-180.75,-180.30]
I don’t work with hydrophone data and I do not know how your data is stored, but I expect you would want to have one single response stage like you are trying to set up, and then that stage would have input units Pa and output units Counts. And yes, you would want absolute amplitude values, such that when you’re in the flat(ish) part of your response and you look at your raw data in counts, when your amplitude would reach the value of your sensitivity as per the calibration information, your input signal would be at 1 Pa amplitude. Another way to see it is, you could simply divide your raw digital counts saved data by a sensitivity value and you’d be looking at Pa. But this is only valid around that frequency, that’s why the general approach is better.
Like I mentioned before, in general you would pick a calibration frequency in the flat response part and then use the sensitivity at that frequency and put it in the InstrumentSensitivity object. Then with the ResponseListElements, the amplitudes would be relative to that overall sensitivity, so you would divide your array by that overall sensitivity to state relative amplitude values.
EDIT: Alternatively you could set the overall sensitivity to 1 and keep your absolute values in your array, but that’s a bit ugly.
I’m also not sure what is gonna happen with your phase information, if it is appropriate to just use zeros like you plan to since you’re lacking that info it looks like. Again, I don’t know hydrophones, but for a seismometer that would definitely cause problems.