Obspy: EVRESP ERROR when using tr.remove_response

[Permian]

Hello,

I get this error message using obspy 'tr.remove_response()' command:

EVRESP ERROR (... [File: ; Start date: ; Stage: 9]):
check_channel; required decimation blockette for IIR or FIR filter missing,
skipping to next response now
* ValueError: check_channel: Illegal RESP format

This happens only when I use the permanent stations from the DK network archived at GEOFON where I have also downloaded the StationXML file from.
The GEOFON team assume an obspy bug....

I am using Python 2.7.6 and obspy 1.1.0

I appreciate your help.

[krischer]

Can you post a small reproducible example that also acquires the data from a web service?

[Permian]

from obspy import UTCDateTime, Stream, read_inventory
from obspy.clients.fdsn import Client

client = Client("GFZ")

pre_filter = [0.005, 0.008, 2, 4]

t0 = UTCDateTime(2015,5,5,10,0,0)
t1 = UTCDateTime(2015,5,5,11,0,0)

bulk = [("DK", "COP", "*", "BHZ", t0, t1),
        ("DK", "BSD", "*", "BHZ", t0, t1),
        ("DK", "MUD", "*", "BHZ", t0, t1)]
inv = client.get_stations_bulk(bulk, level="response")

st = Stream()
st += client.get_waveforms("DK", "COP", "*", "BHZ", t0, t1)

st.attach_response(inv)
st.remove_response(output='VEL', pre_filt=pre_filter, taper_fraction=0.005)

[megies]

This stage is causing the problem:

435           <Stage number="9">
436             <PolesZeros resourceId="ResponsePAZ#20170829123231.311077.14" name="GFZ:DK1980:HDR24_iirpaz_1">
437               <InputUnits>
438                 <Name>COUNTS</Name>
439               </InputUnits>
440               <OutputUnits>
441                 <Name>COUNTS</Name>
442               </OutputUnits>
443               <PzTransferFunctionType>DIGITAL (Z-TRANSFORM)</PzTransferFunctionType>
444               <NormalizationFactor>0.999969</NormalizationFactor>
445               <NormalizationFrequency>1</NormalizationFrequency>
446               <Zero number="1">
447                 <Real>1</Real>
448                 <Imaginary>0</Imaginary>
449               </Zero>
450               <Pole number="0">
451                 <Real>0.99937</Real>
452                 <Imaginary>0</Imaginary>
453               </Pole>
454             </PolesZeros>
455             <StageGain>
456               <Value>1</Value>
457               <Frequency>1</Frequency>
458             </StageGain>
459           </Stage>

[Jollyfant]

Seems like an IIR decimation stage and for PolesZeros the required decimation attributes are not present (and generally unsupported by the StationXML schema). The GFZ webservices also omits these values in the SC3ML output, despite having support for it after verison 0.10.

We can raise a warning instead of throwing but that's about it? I think the decimation attributes are needed.

@andres-h can share some wisdom maybe.

[Permian]

I have contacted Andres before. He has written to me: "stage 9 of DK stations is a digital IIR filter (poles&zeros). There is no decimation and decimation blockette is not required by the FDSN standard"

[Jollyfant]

Yeah, it's not required by the metadata schema but EVRESP needs them or it throws that error. One thing we could do when the decimation blockette it is missing is add a bogus one with a factor of 1 internally when reading inventories.

E.g. pushing this in your code makes it pass:

for net in inv:
  for sta in net:
    for cha in sta:
      cha.response.response_stages[8].decimation_factor = 1
      cha.response.response_stages[8].decimation_delay = 0.0
      cha.response.response_stages[8].decimation_correction = 0.0
      cha.response.response_stages[8].decimation_offset = 0
      cha.response.response_stages[8].decimation_input_sample_rate = cha.response.response_stages[9].decimation_input_sample_rate

Thoughts?

[megies]

Nice catch about the decimation stuff @Jollyfant!

when the decimation blockette it is missing is add a bogus one with a factor of 1 internally when reading inventories

Maybe only do it during the evalresp call? In general we should not automagically change user data I think..

[Permian]

From my point of view that should be an acceptable solution.
However, now I get warning:
WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent.
Execution continuing.

[krischer]

Seems like an IIR decimation stage and for PolesZeros the required decimation attributes are not present (and generally unsupported by the StationXML schema)

All response stages can have the decimation attributes - see the ResponseStageType in the schema.

Maybe only do it during the evalresp call? In general we should not automagically change user data I think..

I'll give it a shot. I'm a bit surprised we did not encounter this any earlier given the massive amount of testing we did for the responses.

[Permian]

Lion, me too. I have requested many data from various networks and clients. This is the only case where the error was encountered.

[krischer]

From my point of view that should be an acceptable solution.
However, now I get warning:
WARNING (norm_resp): computed and reported sensitivities differ by more than 5 percent.
Execution continuing.

That is (hopefully) an unrelated problem. evalresp internally calculates the total sensitivity by just multiplying the sensitivities of all stages. It then compares this to the total sensitivity written in the file and raises this warning if it differs too much (but it always uses the computed one). Usually this is no cause to worry if you deconvolve the whole response. Just to be sure: Do you have a chance to compare the response ObsPy computes (e.g. inv.plot_reponse(0.001)) with a reference solution for these stations?

[Permian]

Well, I have compared the response with evalresp from the rdseed package (so it is not really a reference solution I guess). I have requested the RESP file using:
from obspy.clients.iris import Client
client.resp("DK",`` "BSD", "*", "BHZ", dt, filename='RESP.DK.BSD.BHZ.seed')
However, there is this discrepancy for station BSD of network DK where the sensitivity warning also occurs.

Station ANGG does not show any mismatch.

[krischer]

Hmm- yea something is odd here - can you post the code you used to generate the figures? I don't really fully understand what is plotted there.

[Permian]

from obspy import UTCDateTime, Stream, read_inventory
import matplotlib.pyplot as plt
import os, sys
import numpy as np

netw = 'DK'
stat = 'BSD'
chn = 'BHZ'

## try to request a RESP seed file for DK data
## see http://www.eas.slu.edu/eqc/eqc_cps/TUTORIAL/RESPONSE/index.html for a tutorial
## using gsac and evalsresp for removing the instrument response

from obspy.clients.iris import Client
client = Client()

dt = UTCDateTime("2015-05-05T06:30:00.000")

data = client.resp(netw, stat, "*", chn, dt, filename='RESP.%s.%s.%s.seed' % (netw, stat, chn))

## creates AMP.xxx and PHASE.xxx response files using rdseed package
os.system("evalresp %s %s 2015 125 0.001 10 500 -u 'vel' -f RESP.%s.%s.%s.seed" % (stat, chn, netw, stat, chn))

## load amplitude response from seed RESP file
amp_resp = np.genfromtxt('AMP.%s.%s..%s' % (netw, stat, chn))
phs_resp = np.genfromtxt('PHASE.%s.%s..%s' % (netw, stat, chn))


## create figure with two axis to overlay with the obspy response plot
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212, sharex=ax1)
## plot the seed RESP amplitude response
ax1.loglog(amp_resp[:,0], amp_resp[:,1], 'r', linewidth=2, label='seed RESP')
ax1.grid()
ax1.set_ylabel('Amplitude')
ax2.grid()
ax2.semilogx(phs_resp[:,0], phs_resp[:,1], 'r', linewidth=2, label='seed RESP')
ax2.set_ylabel('Phase (deg)')
ax2.set_xlabel('Frequency (Hz)')
ax1.legend(loc=4)
ax1.set_title('%s.%s.%s' % (netw, stat, chn))


## request the response from GFZ as inventory file
from obspy.clients.fdsn import Client
client = Client("GFZ")

pre_filter = [0.005, 0.008, 2, 4]
t0 = UTCDateTime(2015,5,5,10,0,0)
t1 = UTCDateTime(2015,5,5,11,0,0)

bulk = [(netw, stat, "*", chn, t0, t1)]
inv = client.get_stations_bulk(bulk, level="response")

for net in inv:
  for sta in net:
    for cha in sta:
      cha.response.response_stages[8].decimation_factor = 1
      cha.response.response_stages[8].decimation_delay = 0.0
      cha.response.response_stages[8].decimation_correction = 0.0
      cha.response.response_stages[8].decimation_offset = 0
      cha.response.response_stages[8].decimation_input_sample_rate = cha.response.response_stages[9].decimation_input_sample_rate
      ## plot the obspy inventory response for each channel
      ## see https://docs.obspy.org/_modules/obspy/core/inventory/response.html#Response.plot
      cha.plot(min_freq=0.001, output="VEL", plot_degrees=True, label='obspy', axes=(ax1, ax2))      

st = Stream()
st += client.get_waveforms(netw, stat, "*", "BHZ", t0, t1)
st.attach_response(inv)
st.remove_response(output='VEL', pre_filt=pre_filter, taper_fraction=0.005)

[krischer]

OMG - after a 2 hours debug session I was about to doubt my sanity. Turns out that the StationXML and the RESP files are actually not identical. I'm not sure which one is at fault here, but the RESP file has an extra 9 after the comma, e.g. +9.99937E-01 instead of <Real>0.99937</Real> in XML file...

I'm not sure how anyone is supposed to notice that. And I am honestly a bit surprised how incredibly sensitive this particular response is to changes in the poles. Maybe one should try to find a more stable set of paz for this stage if possible?

@Permian Will you contact the GFZ and tell them about it?

Applying this diff to the RESP files results in equal response curves:

696c696
< B053F15-18     0  +9.99937E-01  +0.00000E+00  +0.00000E+00  +0.00000E+00
---
> B053F15-18     0  +9.9937E-01  +0.00000E+00  +0.00000E+00  +0.00000E+00
930c930

[dsentinel]

OMG - after a 2 hours debug session I was about to doubt my sanity. Turns out that the StationXML and the RESP files are actually not identical. I'm not sure which one is at fault here, but the RESP file has an extra 9 after the comma, e.g. +9.99937E-01 instead of 0.99937 in XML file...

Blasted! I burned some time on this too, and came to the same discrepancy.

However, I think that the error is in the XML for station BSD. Unless it is a tailored response, the xml matches up better with (EDIT: ANGG is a Q330, not a HDR24) the NRL:

nexus lines are nominal from the NRL
altered is the xml downloaded per your method with the decimation fields added

[Permian]

Yes, this is a significant discrepancy in the long-period range.
I have contacted the GEOFON guys.

[krischer]

Thanks everyone! And especially @dsentinel for your detective work. Together with #2128 I consider this issue resolved from our side.

[gaoyajian]

I found a similar mistake in a station in Africa from IRIS. AF.MONG

[megies]

I found a similar mistake in a station in Africa from IRIS. AF.MONG

Report it to IRIS, they will know who to contact about it.

[gaoyajian]

I found a similar mistake in a station in Africa from IRIS. AF.MONG

Report it to IRIS, they will know who to contact about it.

yes