Skip to content

Commit

Permalink
extend remove_response for polynomial responses
Browse files Browse the repository at this point in the history
  • Loading branch information
@barsch
barsch committed Feb 16, 2014
1 parent ef9e60d commit a6c5593
Show file tree
Hide file tree
Showing 4 changed files with 198 additions and 2 deletions.
84 changes: 84 additions & 0 deletions obspy/core/tests/data/stationxml_BK.CMB.__.LKS.xml
View file
@@ -0,0 +1,84 @@
<?xml version="1.0" encoding="ISO-8859-1"?>

<FDSNStationXML xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0" xsi:schemaLocation="http://www.fdsn.org/xml/station/1 http://www.fdsn.org/xml/station/fdsn-station-1.0.xsd">
<Source>IRIS-DMC</Source>
<Sender>IRIS-DMC</Sender>
<Module>IRIS WEB SERVICE: fdsnws-station | version: 1.0.11</Module>
<ModuleURI>http://service.iris.edu/fdsnws/station/1/query?channel=LKS&amp;station=CMB&amp;network=BK&amp;starttime=2004-06-15T00%3A00%3A00.000000&amp;level=response</ModuleURI>
<Created>2014-02-16T15:50:00</Created>
<Network code="BK" startDate="1980-01-01T00:00:00" endDate="2500-12-12T23:59:59" restrictedStatus="open">
<Description>Berkeley Digital Seismograph Network</Description>
<TotalNumberStations>80</TotalNumberStations>
<SelectedNumberStations>1</SelectedNumberStations>
<Station code="CMB" startDate="1996-09-25T19:19:00" endDate="2599-12-31T23:59:59" restrictedStatus="open" alternateNetworkCodes=".EARTHSCOPE,_GSN-BROADBAND,_US-TA,_GSN,_FDSN,_US-REGIONAL,.UNRESTRICTED,_ANSS,_US-REF,_US-ALL,_ANSS-BB,.EARTHSCOPE,_FDSN-ALL,_US-ALL,_REALTIME,_STS-1">
<Latitude>38.03455</Latitude>
<Longitude>-120.38651</Longitude>
<Elevation>697.0</Elevation>
<Site>
<Name>Columbia College, Columbia, CA, USA</Name>
</Site>
<CreationDate>1986-10-25T00:00:00</CreationDate>
<TotalNumberChannels>111</TotalNumberChannels>
<SelectedNumberChannels>2</SelectedNumberChannels>
<Channel locationCode=" " startDate="2004-06-15T00:00:00" restrictedStatus="open" endDate="2010-12-17T00:00:00" code="LKS">
<Latitude>38.03455</Latitude>
<Longitude>-120.38651</Longitude>
<Elevation>697.0</Elevation>
<Depth>2.0</Depth>
<Azimuth>0.0</Azimuth>
<Dip>0.0</Dip>
<Type>CONTINUOUS</Type>
<Type>GEOPHYSICAL</Type>
<SampleRate>1.0</SampleRate>
<ClockDrift>1.0E-4</ClockDrift>
<CalibrationUnits>
<Name>A</Name>
<Description>Amperes</Description>
</CalibrationUnits>
<Sensor>
<Type>YSI 44031 Thermistor</Type>
</Sensor>
<Response>
<InstrumentSensitivity>
<InputUnits>
<Name>C</Name>
<Description>Degrees Centigrade</Description>
</InputUnits>
<OutputUnits>
<Name>COUNTS</Name>
<Description>Digital Counts</Description>
</OutputUnits>
</InstrumentSensitivity>
<Stage number="1">
<Polynomial>
<InputUnits>
<Name>C</Name>
<Description>Degrees Centigrade</Description>
</InputUnits>
<OutputUnits>
<Name>V</Name>
<Description>Volts</Description>
</OutputUnits>
<ApproximationType>MACLAURIN</ApproximationType>
<FrequencyLowerBound>0.0</FrequencyLowerBound>
<FrequencyUpperBound>0.016</FrequencyUpperBound>
<ApproximationLowerBound>3.4</ApproximationLowerBound>
<ApproximationUpperBound>68</ApproximationUpperBound>
<MaximumError>.1</MaximumError>
<Coefficient number="2">12.5091</Coefficient>
<Coefficient number="3">13.8687</Coefficient>
<Coefficient number="4">4.31802</Coefficient>
<Coefficient number="5">1.44163</Coefficient>
<Coefficient number="6">.829113</Coefficient>
<Coefficient number="7">.738976</Coefficient>
<Coefficient number="8">.41094</Coefficient>
<Coefficient number="9">.121068</Coefficient>
<Coefficient number="10">.017856</Coefficient>
<Coefficient number="11">.001043</Coefficient>
</Polynomial>
</Stage>
</Response>
</Channel>
</Station>
</Network>
</FDSNStationXML>
52 changes: 52 additions & 0 deletions obspy/core/tests/data/stationxml_IU.ANTO.30.LDO.xml
View file
@@ -0,0 +1,52 @@
<?xml version="1.0" encoding="ISO-8859-1"?>

<FDSNStationXML xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0" xsi:schemaLocation="http://www.fdsn.org/xml/station/1 http://www.fdsn.org/xml/station/fdsn-station-1.0.xsd">
<Source>IRIS-DMC</Source>
<Sender>IRIS-DMC</Sender>
<Module>IRIS WEB SERVICE: fdsnws-station | version: 1.0.11</Module>
<ModuleURI>http://service.iris.edu/fdsnws/station/1/query?network=IU&amp;level=response&amp;station=ANTO&amp;location=30&amp;starttime=2010-07-23T00%3A00%3A00.000000&amp;channel=LDO</ModuleURI>
<Created>2014-02-16T15:49:04</Created>
<Network code="IU" startDate="1988-01-01T00:00:00" endDate="2500-12-12T23:59:59" restrictedStatus="open">
<Description>Global Seismograph Network (GSN - IRIS/USGS)</Description>
<TotalNumberStations>258</TotalNumberStations>
<SelectedNumberStations>1</SelectedNumberStations>
<Station code="ANTO" startDate="2010-07-23T00:00:00" endDate="2599-12-31T23:59:59" restrictedStatus="open" alternateNetworkCodes="_GSN,.UNRESTRICTED,_FDSN-ALL,_FDSN,_GSN-BROADBAND,_CARIBE-EWS,_REALTIME">
<Latitude>39.868</Latitude>
<Longitude>32.7934</Longitude>
<Elevation>1090.0</Elevation>
<Site>
<Name>Ankara, Turkey</Name>
</Site>
<CreationDate>1992-09-25T00:00:00</CreationDate>
<TotalNumberChannels>60</TotalNumberChannels>
<SelectedNumberChannels>1</SelectedNumberChannels>
<Channel locationCode="30" startDate="2010-07-23T00:00:00" restrictedStatus="open" endDate="2012-11-14T21:27:58" code="LDO">
<Latitude>39.868</Latitude>
<Longitude>32.7934</Longitude>
<Elevation>1090.0</Elevation>
<Depth>0.0</Depth>
<Azimuth>0.0</Azimuth>
<Dip>0.0</Dip>
<Type>CONTINUOUS</Type>
<Type>WEATHER</Type>
<SampleRate>1.0</SampleRate>
<ClockDrift>0.0</ClockDrift>
<Sensor>
<Type>CI/PAS pressure sensor</Type>
</Sensor>
<Response>
<InstrumentPolynomial>
<ApproximationType>MACLAURIN</ApproximationType>
<FrequencyLowerBound>0.0</FrequencyLowerBound>
<FrequencyUpperBound>0.5</FrequencyUpperBound>
<ApproximationLowerBound>80000</ApproximationLowerBound>
<ApproximationUpperBound>110000</ApproximationUpperBound>
<MaximumError>0</MaximumError>
<Coefficient number="1">80000</Coefficient>
<Coefficient number="2">.014305</Coefficient>
</InstrumentPolynomial>
</Response>
</Channel>
</Station>
</Network>
</FDSNStationXML>
35 changes: 35 additions & 0 deletions obspy/core/tests/test_trace.py
View file
Expand Up @@ -12,6 +12,7 @@
import numpy as np
import unittest
import warnings
import os


MATPLOTLIB_VERSION = getMatplotlibVersion()
Expand Down Expand Up @@ -1539,6 +1540,40 @@ def test_remove_response(self):
tr2.remove_response(pre_filt=(0.1, 0.5, 30, 50))
np.testing.assert_array_almost_equal(tr1.data, tr2.data)

def test_remove_polynomial_response(self):
"""
"""
from obspy.station import read_inventory
path = os.path.dirname(__file__)

# blockette 62, stage 0
tr = read()[0]
tr.stats.network = 'IU'
tr.stats.station = 'ANTO'
tr.stats.location = '30'
tr.stats.channel = 'LDO'
tr.stats.starttime = UTCDateTime("2010-07-23T00:00:00")
# remove response
del tr.stats.response
filename = os.path.join(path, 'data', 'stationxml_IU.ANTO.30.LDO.xml')
inv = read_inventory(filename, format='StationXML')
tr.attach_response(inv)
tr.remove_response()

# blockette 62, stage 1 + blockette 58, stage 2
tr = read()[0]
tr.stats.network = 'BK'
tr.stats.station = 'CMB'
tr.stats.location = ''
tr.stats.channel = 'LKS'
tr.stats.starttime = UTCDateTime("2004-06-16T00:00:00")
# remove response
del tr.stats.response
filename = os.path.join(path, 'data', 'stationxml_BK.CMB.__.LKS.xml')
inv = read_inventory(filename, format='StationXML')
tr.attach_response(inv)
tr.remove_response()


def suite():
return unittest.makeSuite(TraceTestCase, 'test')
Expand Down
29 changes: 27 additions & 2 deletions obspy/core/trace.py
View file
Expand Up @@ -12,9 +12,9 @@
from obspy.core.utcdatetime import UTCDateTime
from obspy.core.util import AttribDict, createEmptyDataChunk
from obspy.core.util.base import _getFunctionFromEntryPoint
from obspy.core.util.misc import flatnotmaskedContiguous
from obspy.core.util.decorator import raiseIfMasked, skipIfNoData, \
taper_API_change
from obspy.core.util.misc import flatnotmaskedContiguous
import math
import numpy as np
import warnings
Expand Down Expand Up @@ -2159,7 +2159,7 @@ def remove_response(self, output="VEL", water_level=60, pre_filt=None,
:type taper_fraction: float
:param taper_fraction: Taper fraction of cosine taper to use.
"""
from obspy.station import Response
from obspy.station import Response, PolynomialResponseStage
from obspy.signal.invsim import cosTaper, c_sac_taper, specInv

if "response" not in self.stats:
Expand All @@ -2171,6 +2171,31 @@ def remove_response(self, output="VEL", water_level=60, pre_filt=None,
"(but is of type %s).") % type(self.stats.response)
raise TypeError(msg)


response = self.stats.response
# polynomial response using blockette 62 stage 0
if not response.response_stages and response.instrument_polynomial:
coefficients = response.instrument_polynomial.coefficients
self.data = np.poly1d(coefficients[::-1])(self.data)
return self

# polynomial response using blockette 62 stage 1 and no other stages
if len(response.response_stages) == 1 and \
isinstance(response.response_stages[0], PolynomialResponseStage):
# check for gain
if response.response_stages[0].stage_gain is None:
msg = 'Stage gain not defined for %s - setting it to 1.0'
warnings.warn(msg % self.id)
gain = 1
else:
gain = response.response_stages[0].stage_gain
coefficients = response.response_stages[0].coefficients[:]
for i in range(len(coefficients)):
coefficients[i] /= math.pow(gain, i)
self.data = np.poly1d(coefficients[::-1])(self.data)
return self

# use evalresp
data = self.data.astype("float64")
npts = len(data)
# time domain pre-processing
Expand Down

1 comment on commit a6c5593

@obspy-bot
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This commit has been mentioned on ObsPy Forum. There might be relevant details there:

https://discourse.obspy.org/t/how-to-remove-polynomial-response/1834/8

Please sign in to comment.