Instrument correction (instrument response in dataless file)

Hi,

I am quite new to ObsPy, and I am trying to process a miniseed file by applying instrument correction included in a separate seedless file.
The result I am getting does not look quite right, the time history has a strong baseline deviation which the raw signal did not have; I am not sure if I need to define any more input for the function simulate to work or do any additional processing?

I have done the check just on a single trace of the stream. If I want to extend this to the whole stream, do I just apply st.simulate to the stream and call the relevant dataless file for the other stations sotred in the same miniseed, does the routine understand which instrument correction to apply for each channel?

Many thanks for any help you can give me.

Those are the code lines I am using:

import obspy
from obspy import read
from obspy.io.xseed import Parser

st=read("bdsAutodrm-rU8TsD") # contains the input miniseed file
print(st)
tr7= st[6] # assign trace to new variable

tr7.write('2021-12-27-WOL00-BH1-SH_ASC.asc', format='SH_ASC') # extracting one trace from the stream

parser = Parser("202112270145-0245WOL.dataless")
print(parser)

tr7.simulate(seedresp={'filename': parser, 'units': "ACC"}) #Apply instrument response to trace
tr7.write('2021-12-27-WOL00-BH1-SH_ASC_cor.asc', format='SH_ASC') #saves corrected waveform

You should switch to using read_inventory() and just provide that Inventory object to instrument response removal. Also in general you want to use Stream/Trace.remove_response() if you want to go to physical units. simulate is to simulate your data as recorded by another instrument, different from the one used in the original recording (e.g. simulate a Wood Anderson legacy instrument, etc etc), thats likely not what you want.

If you think something is fishy, please plot your data before and after processing, you can directly attach images here, otherwise people wont be able to comment much without more info.

1 Like

Hi, many thanks for your reply.
I modified the script reading the seedless file as an inventory object and using the trace.remove_response, but I am getting a warning:

C:\Users\davim\Miniconda3\envs\obspy\lib\site-packages\obspy\io\xseed\parser.py:1271: UserWarning: Epoch BN.BKNā€¦BHZ [2008-04-05T00:00:00.000000Z - 2099-01-01T00:00:00.000000Z]: Channel has multiple (but identical) blockettes 58 for stage 0. Only one will be used.

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

Here is the script:
from obspy import read
from obspy.core.inventory import read_inventory

st=read(ā€œbdsAutodrm-rU8TsDā€)
print(st)
tr1= st[0] # assign trace to new variable
tr1.plot()
inv=read_inventory(path_or_file_object=ā€œ202112270145-0245BKN.datalessā€, format=ā€œseedā€, level=ā€˜responseā€™)

tr1_cor=tr1.remove_response(inventory=inv, output=ā€œACCā€)
tr1_cor.plot()

Maybe thereā€™s something wrong with the metadata? I would try to get the metadata again from the people running the experiment, if you know a contactā€¦

Ok, thanks, so you can confirm the script is right

I am also having some issues regarding instrument correction. I have my response files in RESP/dataless seed format. while I am running the code I am getting this error. Can you please explain what mistake I am doing? I am using latest ObsPy version as well.

from obspy import read
from obspy import read_inventory
from obspy import Stream
from obspy.core.inventory.response import InventoryResponse
from obspy import UTCDateTime
from glob import glob
import os

# Paths
data_dir = "/home/lalit/noise_tool_kit/OT23_noise_toolkit/Noise_Toolkit/data/OT23/"
resp_dir = "/home/lalit/obs_resp/"
out_dir = "/home/lalit/Raw_Data_Processing/instrument_corrected"


net = "OT"
stn_codes = ["OT23"]

# Loop over station codes
for stn in stn_codes:
    stn_dir = os.path.join(data_dir, stn)
    # Loop over year and Julian day
    for file_path in sorted(glob(os.path.join(stn_dir, "*"))):
        # Extract year and Julian day from file name
        year = int(file_path.split(".")[-3])
        jday = int(file_path.split(".")[-2])
        # Construct start and end time for request
        start_time = UTCDateTime(year=year, julday=jday)
        end_time = start_time + 86400
        # Loop over channels
        for chn in ["HHZ", "HH1", "HH2"]:
            # Construct file name pattern
            file_pattern = f"OT.{stn}..{chn}.{year}.{jday:03d}.*.*.mseed"
            # Find files matching pattern
            file_paths = sorted(glob(os.path.join(stn_dir, file_pattern)))
            # Skip if no matching files found
            if not file_paths:
                continue
            # Read and merge data
            st = Stream()
            for file_path in file_paths:
                st += read(file_path)
            st.merge()
            # Get instrument response
            resp_pattern = f"RESP.{net}.{stn}..{chn}"
            resp_path = os.path.join(resp_dir, resp_pattern)
            inv = read_inventory(resp_path)
            response = InventoryResponse.get_evalresp_response(inv, start_time, output="VEL")
            # Remove response and correct for sensitivity
            st.remove_response(inventory=inv, output="VEL", pre_filt=None, water_level=60, response=response)
            # Write to file
            out_path = os.path.join(out_dir, os.path.basename(file_paths[0]))
            st.write(out_path, format="MSEED")

I am getting this error

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
/tmp/ipykernel_2555740/169522340.py in <module>
      2 from obspy import read_inventory
      3 from obspy import Stream
----> 4 from obspy.core.inventory.response import InventoryResponse
      5 from obspy.clients.fdsn import Client
      6 from obspy import UTCDateTime

ImportError: cannot import name 'InventoryResponse' from 'obspy.core.inventory.response' (/home/userlibs/miniconda3/envs/ispaq/lib/python3.8/site-packages/obspy/core/inventory/response.py)

There is no such object as InventoryResponse in our code, hence the import error. For instrument correction you do not need to go to lower level API, you can stay at the highest level and just use core objects/functions, read (and the resulting Stream object) and read_inventory (and the resulting Inventory object should be all you have to interact with.

see remove_response

1 Like