Station response issues

nbalfour commented 7 years ago

This is an issue that I haven't had time to get my head around but I will provide as much information about the problem that I can.

Firstly, where I used python I used python3.

My first test was to attach the response metadata to the waveform data I collected try to remove the response.

from obspy import UTCDateTime
from obspy.clients.fdsn import Client as FDSN_Client

client = FDSN_Client("http://beta-service.geonet.org.nz/")
t = UTCDateTime("2016-09-01T16:37:00.000")
st = client.get_waveforms("NZ", "TDHS","20", "?N?", t, t + 300,attach_response=True)
pre_filt = (0.005, 0.006, 30.0, 35.0)
st.remove_response(output='ACC', pre_filt=pre_filt)
st.plot()

This is the error it returned:

---------------------------------------------------------------------------
ObsPyException                            Traceback (most recent call last)
<ipython-input-1-a005388685ac> in <module>()
     10 st = client.get_waveforms("NZ", "TDHS","20", "?N?", t, t + 300,attach_response=True)
     11 pre_filt = (0.005, 0.006, 30.0, 35.0)
---> 12 st.remove_response(output='ACC', pre_filt=pre_filt)
     13 st.plot()
     14 pga = st.max()

C:\Users\natalieb\AppData\Local\Continuum\Anaconda3\lib\site-packages\obspy\core\stream.py in remove_response(self, *args, **kwargs)
   3029         """
   3030         for tr in self:
-> 3031             tr.remove_response(*args, **kwargs)
   3032         return self
   3033 

<decorator-gen-162> in remove_response(self, inventory, output, water_level, pre_filt, zero_mean, taper, taper_fraction, plot, fig, **kwargs)

C:\Users\natalieb\AppData\Local\Continuum\Anaconda3\lib\site-packages\obspy\core\trace.py in _add_processing_info(func, *args, **kwargs)
    230     info = info % "::".join(arguments)
    231     self = args[0]
--> 232     result = func(*args, **kwargs)
    233     # Attach after executing the function to avoid having it attached
    234     # while the operation failed.

C:\Users\natalieb\AppData\Local\Continuum\Anaconda3\lib\site-packages\obspy\core\trace.py in remove_response(self, inventory, output, water_level, pre_filt, zero_mean, taper, taper_fraction, plot, fig, **kwargs)
   2703         freq_response, freqs = \
   2704             response.get_evalresp_response(self.stats.delta, nfft,
-> 2705                                            output=output, **kwargs)
   2706 
   2707         if plot:

C:\Users\natalieb\AppData\Local\Continuum\Anaconda3\lib\site-packages\obspy\core\inventory\response.py in get_evalresp_response(self, t_samp, nfft, output, start_stage, end_stage)
    768             msg = ("Can not use evalresp on response with no response "
    769                    "stages.")
--> 770             raise ObsPyException(msg)
    771 
    772         import obspy.signal.evrespwrapper as ew

ObsPyException: Can not use evalresp on response with no response stages.

This sort of makes sense in that when I do the same request using the URL http://beta-service.geonet.org.nz/fdsnws/station/1/query?station=TDHS&channel=HNZ&level=response It has no response stage information.

[ ] Can someone please check that the input StationXML for beta is actually missing this information. If it is provided in the input StationXML then FDSN is doings something odd that removes it.

I then did the same thing with a different station... KIKS It returns the response information for the URL query (below) but still gets the same error with obspy. 'http://beta-service.geonet.org.nz/fdsnws/station/1/query?station=KIKS&channel=HNZ&level=response'

My impression is that there are a couple of problems going on here.

ozym commented 7 years ago

As far as I can tell it should be there, via the station.xml file in the api bucket. Here's the last response chunk for HNZ.

    <Station code="TDHS" startDate="2002-04-15T00:00:00" endDate="9999-01-01T00:00:00" restrictedStatus="open">
      <Description>National strong motion network</Description>
      <Comment id="1">
        <Value>Location is given in WGS84</Value>
      </Comment>
      <Latitude datum="WGS84">-37.633319923</Latitude>
      <Longitude datum="WGS84">178.365371742</Longitude>
      <Elevation>5</Elevation>
      <Site>
        <Name>Te Araroa District High School</Name>
        <Description>within 5 km of Te Araroa</Description>
      </Site>
      <CreationDate>2002-04-15T00:00:00</CreationDate>
...
      <Channel code="HNZ" startDate="2014-07-21T00:00:05" endDate="9999-01-01T00:00:00" restrictedStatus="open" locationCode="20">
        <Comment id="1">
          <Value>Location estimated from internal GPS clock</Value>
        </Comment>
        <Comment id="2">
          <Value>Location is given in WGS84</Value>
        </Comment>
        <Comment id="3">
          <Value>Sensor orientation not known</Value>
        </Comment>
        <Latitude datum="WGS84">-37.633319923</Latitude>
        <Longitude datum="WGS84">178.365371742</Longitude>
        <Elevation>5</Elevation>
        <Depth>0</Depth>
        <Azimuth>0</Azimuth>
        <Dip>-90</Dip>
        <Type>TRIGGERED</Type>
        <Type>GEOPHYSICAL</Type>
        <SampleRate>200</SampleRate>
        <SampleRateRatio>
          <NumberSamples>200</NumberSamples>
          <NumberSeconds>1</NumberSeconds>
        </SampleRateRatio>
        <StorageFormat>Steim2</StorageFormat>
        <ClockDrift>0.0001</ClockDrift>
        <Sensor resourceId="Sensor#CUSP3D:42017">
          <Type>Strong Motion Sensor</Type>
          <Description>CUSP 3D</Description>
          <Model>CUSP3D</Model>
          <SerialNumber>42017</SerialNumber>
          <InstallationDate>2014-07-21T00:00:05</InstallationDate>
        </Sensor>
        <DataLogger resourceId="Datalogger#CUSP3D:">
          <Type>Strong Motion Recorder</Type>
          <Description>Cusp3D</Description>
          <Manufacturer>CSI</Manufacturer>
          <Model>CUSP3D</Model>
          <InstallationDate>2014-07-21T00:00:05</InstallationDate>
        </DataLogger>
        <Response>
          <InstrumentSensitivity>
            <Value>101971.62129779284</Value>
            <Frequency>1</Frequency>
            <InputUnits>
              <Name>m/s**2</Name>
            </InputUnits>
            <OutputUnits>
              <Name>count</Name>
            </OutputUnits>
          </InstrumentSensitivity>
          <Stage number="1">
            <PolesZeros resourceId="PolesZeros#CUSP SENSOR" name="TDHS.20.HNZ.2014.202.stage_1">
              <InputUnits>
                <Name>m/s**2</Name>
              </InputUnits>
              <OutputUnits>
                <Name>V</Name>
              </OutputUnits>
              <PzTransferFunctionType>LAPLACE (RADIANS/SECOND)</PzTransferFunctionType>
              <NormalizationFactor>1</NormalizationFactor>
              <NormalizationFrequency>1</NormalizationFrequency>
            </PolesZeros>
            <StageGain>
              <Value>0.10197162129779283</Value>
              <Frequency>1</Frequency>
            </StageGain>
          </Stage>
          <Stage number="2">
            <Coefficients resourceId="Coefficients#CUSP-200" name="TDHS.20.HNZ.2014.202.stage_2">
              <InputUnits>
                <Name>V</Name>
              </InputUnits>
              <OutputUnits>
                <Name>count</Name>
              </OutputUnits>
              <CfTransferFunctionType>DIGITAL</CfTransferFunctionType>
            </Coefficients>
            <Decimation>
              <InputSampleRate>200</InputSampleRate>
              <Factor>1</Factor>
              <Offset>0</Offset>
              <Delay>0</Delay>
              <Correction>0</Correction>
            </Decimation>
            <StageGain>
              <Value>1e+06</Value>
              <Frequency>1</Frequency>
            </StageGain>
          </Stage>
        </Response>
      </Channel>
...

nbalfour commented 7 years ago

@ozym please compare this against the output from service.geonet.org.nz. They are quite different

junghao commented 7 years ago

I've tested the service with a station xml file in my computer (copied on 02/Jun) and it does return stages.

ozym commented 7 years ago

@nbalfour here's the output from the old system (c.f. the input above) and see far below for beta.

      <Channel code="HNZ" startDate="2014-07-21T00:00:05" restrictedStatus="open" locationCode="20">
        <Latitude>-37.63332</Latitude>
        <Longitude>178.365372</Longitude>
        <Elevation>5</Elevation>
        <Depth>0</Depth>
        <Azimuth>90</Azimuth>
        <Dip>-90</Dip>
        <SampleRate>200</SampleRate>
        <SampleRateRatio>
          <NumberSamples>200</NumberSamples>
          <NumberSeconds>1</NumberSeconds>
        </SampleRateRatio>
        <StorageFormat>Steim2</StorageFormat>
        <ClockDrift>0.0001</ClockDrift>
        <Sensor resourceId="Sensor#20150130114231.370299.49546">
          <Type>CSI CUSP3D</Type>
          <Description>CUSP3D</Description>
          <Manufacturer>CSI</Manufacturer>
          <Model>CUSP3D</Model>
        </Sensor>
        <DataLogger resourceId="Datalogger#20150130114231.399153.49580">
          <Description>TDHS.2014.202.HNZ20</Description>
        </DataLogger>
        <Response>
          <InstrumentSensitivity>
            <Value>101972</Value>
            <Frequency>1</Frequency>
            <InputUnits>
              <Name>M/S**2</Name>
            </InputUnits>
            <OutputUnits>
              <Name>COUNTS</Name>
            </OutputUnits>
          </InstrumentSensitivity>
          <Stage number="1">
            <PolesZeros resourceId="ResponsePAZ#20150130114231.370651.49547" name="TDHS.2014.202.NZ20">
              <InputUnits>
                <Name>M/S**2</Name>
              </InputUnits>
              <OutputUnits>
                <Name>V</Name>
              </OutputUnits>
              <PzTransferFunctionType>LAPLACE (RADIANS/SECOND)</PzTransferFunctionType>
              <NormalizationFactor>1</NormalizationFactor>
              <NormalizationFrequency>1</NormalizationFrequency>
            </PolesZeros>
            <StageGain>
              <Value>0.101972</Value>
              <Frequency>1</Frequency>
            </StageGain>
          </Stage>
          <Stage number="2">
            <Coefficients>
              <InputUnits>
                <Name>V</Name>
              </InputUnits>
              <OutputUnits>
                <Name>COUNTS</Name>
              </OutputUnits>
              <CfTransferFunctionType>DIGITAL</CfTransferFunctionType>
            </Coefficients>
            <Decimation>
              <InputSampleRate>200</InputSampleRate>
              <Factor>1</Factor>
              <Offset>0</Offset>
              <Delay>0</Delay>
              <Correction>0</Correction>
            </Decimation>
            <StageGain>
              <Value>1000000</Value>
              <Frequency>0</Frequency>
            </StageGain>
          </Stage>
        </Response>
      </Channel>

However, here's the output of the currently running beta which appears to be missing the response stages.

        <Response>
          <InstrumentSensitivity>
            <Value>101971.62129779284</Value>
            <Frequency>1</Frequency>
            <InputUnits>
              <Name>m/s**2</Name>
            </InputUnits>
            <OutputUnits>
              <Name>count</Name>
            </OutputUnits>
          </InstrumentSensitivity>
        </Response>

mabznz commented 7 years ago

That station XML is fine but its whatever the obspy 'get_waveforms()' call gets from the FDSN service, Presume it a dataselect query returning miniSEED? Not knowing a lot, does that miniSEED data have the channel/station data imbedded in it in beta?

st = client.get_waveforms('NZ','TDHS','20','BN1',t1,t2,attach_response=True)
tr = st[0]
print(tr.stats.response)

With present production you get:

Channel Response
    From M/S**2 () to COUNTS ()
    Overall Sensitivity: 101972 defined at 1.000 Hz
    4 stages:
        Stage 1: PolesZerosResponseStage from M/S**2 to V, gain: 0.101972
        Stage 2: CoefficientsTypeResponseStage from V to COUNTS, gain: 1e+06
        Stage 3: FIRResponseStage from COUNTS to COUNTS, gain: 1
        Stage 4: FIRResponseStage from COUNTS to COUNTS, gain: 1

With beta you get:

Channel Response
    From m/s**2 () to count ()
    Overall Sensitivity: 101972 defined at 1.000 Hz
    0 stages:

mabznz commented 7 years ago

Yep, sorry, ignore my last comment. Was looking at wrong channel in station xml. Some channels have stages some don't.

mabznz commented 7 years ago

@junghao @gclitheroe

This may be just because stationXML in beta is outdated?

I ran stationxml out of delta for small subset of stations and ran FDSN locally. Ran tests locally and against beta.

If Beta running latest, code seems fine, input should be the problem.

mabznz commented 7 years ago

Geoff updated stationXML in beta and confirmed it was a data issue.

junghao commented 7 years ago

At this moment the response from beta service doesn't contain <Stage> http://beta-service.geonet.org.nz/fdsnws/station/1/query?station=KARZ&network=NZ&channel=EHZ&level=response

But the production service does: http://service.geonet.org.nz/fdsnws/station/1/query?station=KARZ&network=NZ&channel=EHZ&level=response

And I've downloaded the station source xml from S3 then ran the service locally it contains <Stage>.

@gclitheroe

GeoNet / fdsn

Station response issues #62