Response stage numbering with gaps

Development
1

Hi, I am a developer working with Blacknest (https://www.blacknest.gov.uk) on the Blacknest Data System (BDS), https://portal.beam.ltd.uk/bdspublic.
The BDS stores seismic data and Metadata from various sources including a lot of legacy data and Metadata. The seismologists are using ObsPy for data analysis work and are having a slight issue with some of the Metadata from certain sources. Some of the Metadata has response stage numbers that, although in order, have gaps in the numbering sequence. I believe this is because the stage numbers were used to represent a particular type of stage that may or may not have been present. The gaps cause issues in ObsPy as it requires the stage numbers to be a complete integer sequence with no gaps.
In the BDS we try an keep all of the data and Metadata in its original state with as little manipulation as possible as one of its purposes is an archive. Now we can re-order the response stage numbering on export, to say SEED or StationXml formats, but we wondered if people thought that it would be worthwhile for ObsPy to handle response stage numbering with gaps in the integer sequence so it could directly read some of these Metadata file directly even if not supplied from the BDS ?
A brief look at the ObsPy code seems to suggest that this feature could be added fairly easily either by re-ordering on import or allowing gaps in the response stage numbers.

2

Both the SEED standard and the StationXML schema are very explicit about this. Responses should have a sequence of stages with sequence numbers starting at 1 and iterating sequentially from there.
I’d argue that renumbering stages properly on export is the way to go here.

3

Certainly we can renumber on export from the BDS, but I suspect some users will have Metadata from various platforms where this might happen and where there “might” be reasons for it.

The response_stagenumber_reorder function in the following code will renumber the stages, if someone is in this boat so it is easier to use ObsPy with the Metadata sources they have:

#! /usr/bin/python3
# Read in and deconvolve seismogram.
#
import sys
import getopt
from obspy import read, read_inventory

def response_stagenumber_reorder(inv):
    """ Re-order the response stage numbers starting at 1, assumes list is in response stage order """
    for network in inv:
        for station in network:
            for channel in station:
                sn = 1;
                sn_last = 0;
                for stage in channel.response.response_stages:
                    if(stage.stage_sequence_number <= sn_last):
                        raise ValueError("Response stage numbering sequence issue");

                    sn_last = stage.stage_sequence_number;
                    stage.stage_sequence_number = sn;
                    sn += 1;

def main():
    st = read('202410310505-0510I49H1.mseed')
    for tr in st:
        print (tr.stats.calib)

    inv = read_inventory('20241031I49H1.xml')
    #inv = read_inventory('20241031I49H1-ordered.xml')

    # Re-order response stage numbers
    response_stagenumber_reorder(inv);
    
    #print(inv.get_contents());
    #print(inv[0][0][0]);

    print("Overal response");
    print(inv[0][0][0].response)
    
    prefilt=(0.1, 0.2, 7.5, 9.9)
    st.remove_response(inventory=inv, pre_filt=prefilt, output="DISP", water_level=60, plot=True)
    st.plot()

if __name__ == "__main__":
    main()