From bf7533431bd26fbf6a668a620b27c75705bd1472 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 13:16:56 -0700 Subject: [PATCH 01/26] New format for sensor JSON file MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit I don't have the model or frequency of sensitivity measurement for all sensors — where these are missing, I've left the fields as "?"" for sensor model and -9999 for frequency. These must be filled in later! --- sensor_sensitivities.json | 96 ++++++++++++++++++++++++++++++++------- 1 file changed, 80 insertions(+), 16 deletions(-) diff --git a/sensor_sensitivities.json b/sensor_sensitivities.json index d29aff7..a7ba515 100644 --- a/sensor_sensitivities.json +++ b/sensor_sensitivities.json @@ -1,18 +1,82 @@ { -"SN46" : 0.00902, -"SN47" : 0.00898, -"SN48" : 0.00909, -"SN49" : 0.00904, -"SN50" : 0.00905, -"SN51" : 0.00902, -"SN73" : 0.00898, -"SN74" : 0.00903, -"SN80" : 0.00905, -"SN81" : 0.00898, -"SN82" : 0.00901, -"SN83" : 0.00902, -"SN84" : 0.00902, -"SN85" : 0.00902, -"SN86" : 0.00898, -"SN87" : 0.00903 + "SN46": { + "model": "?", + "sensitivity": 0.00902, + "frequency": -9999 + }, + "SN47": { + "model": "?", + "sensitivity": 0.00898, + "frequency": -9999 + }, + "SN48": { + "model": "60-UHP", + "sensitivity": 0.00909, + "frequency": -9999 + }, + "SN49": { + "model": "60-UHP", + "sensitivity": 0.00904, + "frequency": -9999 + }, + "SN50": { + "model": "60-UHP", + "sensitivity": 0.00905, + "frequency": -9999 + }, + "SN51": { + "model": "60-UHP", + "sensitivity": 0.00902, + "frequency": -9999 + }, + "SN73": { + "model": "?", + "sensitivity": 0.00898, + "frequency": -9999 + }, + "SN74": { + "model": "60-UHP", + "sensitivity": 0.00903, + "frequency": -9999 + }, + "SN80": { + "model": "?", + "sensitivity": 0.00905, + "frequency": -9999 + }, + "SN81": { + "model": "?", + "sensitivity": 0.00898, + "frequency": -9999 + }, + "SN82": { + "model": "?", + "sensitivity": 0.00901, + "frequency": -9999 + }, + "SN83": { + "model": "?", + "sensitivity": 0.00902, + "frequency": -9999 + }, + "SN84": { + "model": "?", + "sensitivity": 0.00902, + "frequency": -9999 + }, + "SN85": { + "model": "?", + "sensitivity": 0.00902, + "frequency": -9999 + }, + "SN86": { + "model": "?", + "sensitivity": 0.00898, + "frequency": -9999 + }, + "SN87": { + "model": "?", + "sensitivity": 0.00903, + "frequency": -9999 + } } From 3b177ea535d948637ca7b630cf8f144599621c71 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 13:21:34 -0700 Subject: [PATCH 02/26] Slight update to reflect changes --- README.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 2d374f0..4afe788 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,8 @@ cube_conversion This command-line tool converts [DiGOS](https://digos.eu/) DATA-CUBE3 files into miniSEED files of a desired length of time with specified metadata. Output miniSEED files have units of Pa, unless the user selects to export the files in -a form suitable for submission to EarthScope (formerly IRIS). The tool +a form suitable for submission to EarthScope (formerly IRIS). A script to produce +StationXML files is also provided, which is useful for EarthScope uploads. The tool can differentiate between channels for 3 channel DATA-CUBE3 files and optionally extract coordinates from the digitizer's GPS. The code only looks for files from digitizers defined in the `digitizer_sensor_pairs.json` file. Therefore, @@ -72,12 +73,13 @@ Supplemental files * `digitizer_offsets.json` — Digitizer offsets in V (We have found that each digitizer has a slight voltage offset from zero) -* `sensor_sensitivities.json` — Infrasound sensor sensitivities in V/Pa +* `sensor_sensitivities.json` — Infrasound sensor model, sensitivity in V/Pa, + and frequency of sensitivity measurement in Hz Usage ----- -To print the script's help menu, execute the following terminal commands: +To print the conversion script's help menu, execute the following terminal commands: ``` conda activate uafinfra # Or your pre-existing env cube_convert --help # Or: python /path/to/cube_conversion/cube_convert.py --help From ad2a8fa7ed7fb0f8fccdf9dcfba77555d6165598 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 13:30:16 -0700 Subject: [PATCH 03/26] Add skeleton of 2nd script --- cube_conversion/cube_stationxml.py | 37 ++++++++++++++++++++++++++++++ pyproject.toml | 1 + 2 files changed, 38 insertions(+) create mode 100755 cube_conversion/cube_stationxml.py diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py new file mode 100755 index 0000000..80c4906 --- /dev/null +++ b/cube_conversion/cube_stationxml.py @@ -0,0 +1,37 @@ +#!/usr/bin/env python + +import argparse +import os + +# ----------------------------------------------------------------------------- +# Advanced configuration options +# ----------------------------------------------------------------------------- + +# ----------------------------------------------------------------------------- + + +# Define callable main function to work with [project.scripts] +def main(): + + # Set up command-line interface + parser = argparse.ArgumentParser( + description='Generate StationXML files from DATA-CUBE³ miniSEED files and metadata.', + allow_abbrev=False, + ) + parser.add_argument( + 'input_dir', + help='directory containing miniSEED files and coordinate files produced by cube_convert', + ) + input_args = parser.parse_args() + + # Check if input directory/ies is/are valid + for input_dir in input_args.input_dir: + if not os.path.exists(input_dir): + raise NotADirectoryError( + f'Input directory \'{input_dir}\' doesn\'t ' 'exist.' + ) + + +# Run the main function if this is called as a script +if __name__ == '__main__': + main() diff --git a/pyproject.toml b/pyproject.toml index fa190c0..5ab6ca4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,3 +8,4 @@ build-backend = "setuptools.build_meta" [project.scripts] cube_convert = "cube_conversion.cube_convert:main" +cube_stationxml = "cube_conversion.cube_stationxml:main" From b7cc7eec7d6715b6b832b60af93213f3c6b884a8 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 13:38:11 -0700 Subject: [PATCH 04/26] Update for single input dir --- cube_conversion/cube_stationxml.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 80c4906..45254d8 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -24,12 +24,11 @@ def main(): ) input_args = parser.parse_args() - # Check if input directory/ies is/are valid - for input_dir in input_args.input_dir: - if not os.path.exists(input_dir): - raise NotADirectoryError( - f'Input directory \'{input_dir}\' doesn\'t ' 'exist.' - ) + # Check if input directory is valid + if not os.path.exists(input_args.input_dir): + raise NotADirectoryError( + f'Input directory \'{input_args.input_dir}\' doesn\'t ' 'exist.' + ) # Run the main function if this is called as a script From 1f07ee235071e97968021eac438f666ac2ebb641 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 15:18:10 -0700 Subject: [PATCH 05/26] Add note on usage and assumptions --- cube_conversion/cube_stationxml.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 45254d8..30fe872 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -1,5 +1,23 @@ #!/usr/bin/env python +""" +Make StationXML files from scratch, given an existing directory of miniSEED files and +metadata generated by cube_convert. Adopted from the ObsPy tutorial here: + https://docs.obspy.org/tutorial/code_snippets/stationxml_file_from_scratch.html + +The code assumes that your sensors are all Chaparral Physics products, and that your +digitizers are all DiGOS DATA-CUBE³s. + +Importantly, this code *will not work properly* unless you're using ObsPy after commit +107666b: + https://github.com/obspy/obspy/tree/107666b9659cbfdc98470a47c003772377551264 +(This commit introduced a change that allows for re-computation of the overall +sensitivity for non-standard units — like Pa.) Until the 1.5.0 release of ObsPy comes +out, you can simply update ObsPy via: + pip install git+https://github.com/obspy/obspy.git +This will fetch and install the newest copy from GitHub. +""" + import argparse import os From 5c04a02fdd175ba44de9abdd69bf9accc8a4021a Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 17:34:20 -0700 Subject: [PATCH 06/26] KEY: Flesh out StationXML creation script Based on https://github.com/liamtoney/yasur_2023/blob/main/code/make_yasur_infra_stationxml.py --- cube_conversion/cube_stationxml.py | 239 +++++++++++++++++++++++++++-- 1 file changed, 230 insertions(+), 9 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 30fe872..6128c5d 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -5,8 +5,9 @@ metadata generated by cube_convert. Adopted from the ObsPy tutorial here: https://docs.obspy.org/tutorial/code_snippets/stationxml_file_from_scratch.html -The code assumes that your sensors are all Chaparral Physics products, and that your -digitizers are all DiGOS DATA-CUBE³s. +The code assumes that your sensors are all Chaparral Physics products, that your +digitizers are all DiGOS DATA-CUBE³s, and — IMPORTANTLY — that you only have one sensor +connected to each digitizer (i.e., one channel per station). Importantly, this code *will not work properly* unless you're using ObsPy after commit 107666b: @@ -19,13 +20,37 @@ """ import argparse -import os +import json +from pathlib import Path -# ----------------------------------------------------------------------------- +from obspy import Stream, read +from obspy.clients.nrl import NRL +from obspy.core.inventory import ( + Channel, + Equipment, + Inventory, + Network, + PolesZerosResponseStage, + Site, + Station, +) + +# -------------------------------------------------------------------------------------- # Advanced configuration options -# ----------------------------------------------------------------------------- +# -------------------------------------------------------------------------------------- +GAIN = 1 # DiGOS DATA-CUBE³ amplifier gain (this should usually be 1!) +BOB_FACTOR = 10 # Breakout box factor for DATA-CUBE³s (this should usually be 10!) +# -------------------------------------------------------------------------------------- -# ----------------------------------------------------------------------------- +# These names must match EXACTLY what is in the NRL. Unforunately, there are differences +# between the NRL accessed through ObsPy (via online) versus a local copy of the NRL. +# These are likely related to NRL v1 versus v2, but it's not clear. I think the v2 ones +# should be preferred, and they can be viewed here: https://ds.iris.edu/ds/nrl/ +# These show up first in the tuples below, the second entries are for the older NRL v1 +# and we include those just to cover the bases. +_SENSOR_MANUFACTURER = 'Chaparral', 'Chaparral Physics' +_DATALOGGER_MANUFACTURER = 'DiGOSOmnirecs', 'DiGOS/Omnirecs' +_DATALOGGER_MODEL = 'DataCube', 'DATACUBE' # Define callable main function to work with [project.scripts] @@ -40,13 +65,209 @@ def main(): 'input_dir', help='directory containing miniSEED files and coordinate files produced by cube_convert', ) + parser.add_argument( + 'station_mappings', + nargs='+', + help='one or more mappings of the form STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2', + ) + parser.add_argument( + 'output_filename', help='filename for the output StationXML file (full path)' + ) + parser.add_argument( + '--nrl-path', + default=None, + help='path to local copy of the NRL (if not provided, use online NRL)', + ) input_args = parser.parse_args() # Check if input directory is valid - if not os.path.exists(input_args.input_dir): - raise NotADirectoryError( - f'Input directory \'{input_args.input_dir}\' doesn\'t ' 'exist.' + input_dir = Path(input_args.input_dir) + if not input_dir.is_dir(): + raise NotADirectoryError(f'Input directory \'{input_dir}\' doesn\'t ' 'exist.') + + # Parse mappings into dictionary with keys being the station codes (TODO: validate!) + station_mappings = {} + for mapping in input_args.station_mappings: + station_code, cube_name, sensor_serial = mapping.split(':') + station_mappings[station_code] = { + 'cube_name': cube_name, + 'sensor_serial': sensor_serial, + } + + # Check if NRL path is a directory if provided + nrl_path = input_args.nrl_path + if nrl_path is not None: + nrl_path = Path(nrl_path).absolute() + if not nrl_path.is_dir(): + raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') + + # Find root directory for cube_conversion repo + root_dir = Path(__file__).parents[1] + + # Load sensor information including sensitivities, model numbers, etc. + with open(root_dir / 'sensor_sensitivities.json') as f: + sensor_info = json.load(f) + + # Read in all data into Stream object for processing + st = Stream() + for mseed_file in sorted(input_dir.glob('??.*.*.???.????.???.??')): + st += read(mseed_file) + + # Process SEED network code + network_codes = set(tr.stats.network for tr in st) + assert len(network_codes) == 1, 'Multiple network codes found in data!' + network_code = network_codes.pop() + + # Process SEED station codes + station_codes = sorted(set(tr.stats.station for tr in st)) + + # Check SEED location codes, as these can indicate if multiple sensors are connected + # to the same digitizer (not supported by this script) + location_codes = sorted(set(tr.stats.location for tr in st)) + assert len(location_codes) == 1, 'Multiple location codes found in data!' + location_code = location_codes.pop() + + # Create Inventory with all required components + inv = Inventory() + net = Network(code=network_code) + for station in station_codes: + # Get station start and end times, SEED channel code, sample rate... + st_station = st.select(station=station) + station_starttime = min(tr.stats.starttime for tr in st_station) + station_endtime = max(tr.stats.endtime for tr in st_station) + channel_codes = sorted(set(tr.stats.channel for tr in st_station)) + assert ( + len(channel_codes) == 1 + ), f'Multiple channel codes found in for station {station}!' + channel_code = channel_codes.pop() + sample_rates = sorted(set(tr.stats.sampling_rate for tr in st_station)) + assert len(sample_rates) == 1, 'Multiple sample rates found in data!' + sample_rate = sample_rates.pop() + # Read in coordinates from JSON file + coord_json_file = sorted( + input_dir.glob( + f'{network_code}.{station}.{location_code}.{channel_code}.json' + ) + ) + if len(coord_json_file) == 0: + raise FileNotFoundError( + f'No coordinate JSON file found for station {station}!' + ) + elif len(coord_json_file) > 1: + raise ValueError # Almost impossible, but just in case + else: + coord_json_file = coord_json_file[0] + with open(coord_json_file) as f: + latitude, longitude, elevation = json.load(f) + # Make Station object + sta = Station( + code=station, + latitude=latitude, + longitude=longitude, + elevation=elevation, + site=Site(name=station), # Bare minumum... could make this more detailed + start_date=station_starttime, + end_date=station_endtime, + ) + # Define sensor Equipment object + serial_number = station_mappings[station]['sensor_serial'] + sensor_model = sensor_info[serial_number]['model'] + sensor = Equipment( + type='Infrasound sensor', + description=f'{_SENSOR_MANUFACTURER[0]} {sensor_model}', # MDA shows this! + manufacturer=_SENSOR_MANUFACTURER[0], + model=sensor_model, + serial_number=serial_number, ) + # Define digitizer Equipment object + cube_name = station_mappings[station]['cube_name'] + data_logger = Equipment( + type='Digitizer', + manufacturer=_DATALOGGER_MANUFACTURER[0], + model=_DATALOGGER_MODEL[0], + serial_number=cube_name, + ) + # Make Channel object + cha = Channel( + code=channel_code, + location_code=location_code, + latitude=latitude, + longitude=longitude, + elevation=elevation, + depth=0, # Required, always 0? + sample_rate=sample_rate, + start_date=station_starttime, + end_date=station_endtime, + sensor=sensor, + data_logger=data_logger, + ) + # Access the NRL to get response information. If the user provided a local path + # to the NRL, use that; otherwise, use the online NRL. + if nrl_path is not None: + nrl = NRL(nrl_path) + option_ind = 0 # Use v2 NRL entries + else: + nrl = NRL() + option_ind = 1 # Use v1 NRL entries + # The contents of the NRL can be explored interactively in a Python prompt, see + # API documentation of NRL submodule: + # http://docs.obspy.org/packages/obspy.clients.nrl.html + # Here we assume that the end point of data logger and sensor are already known. + + # Get the nominal response for this combination of sensor and digitizer + response = nrl.get_response( + sensor_keys=[_SENSOR_MANUFACTURER[option_ind], sensor_model], + datalogger_keys=[ + _DATALOGGER_MANUFACTURER[option_ind], + _DATALOGGER_MODEL[option_ind], + str(GAIN), + str(sample_rate), + ], + ) + # KEY: Add a response stage which applies the breakout box factor, after first + # stage (i.e., right after Pa --> V) — this should be a PZ stage, see: + # https://docs.fdsn.org/projects/stationxml/en/latest/reference.html#response-stage + bob_stage_sequence_number = 2 + bob_stage = PolesZerosResponseStage( + description='DiGOS breakout box voltage step-down', + stage_sequence_number=bob_stage_sequence_number, + stage_gain=1 / BOB_FACTOR, + stage_gain_frequency=1, # [Hz] Anywhere where the response is flat? + input_units='V', + output_units='V', + pz_transfer_function_type='LAPLACE (RADIANS/SECOND)', + normalization_frequency=0, + normalization_factor=1, + poles=[], + zeros=[], + ) + response.response_stages.insert(bob_stage_sequence_number - 1, bob_stage) + for response_stage in response.response_stages[bob_stage_sequence_number:]: + response_stage.stage_sequence_number += 1 # Increment following stage #s + # KEY: Update with our specific measured sensor sensitivity + nominal_sensitivity = response.response_stages[0].stage_gain + measured_sensitivity = sensor_info[serial_number]['sensitivity'] + at_frequency = sensor_info[serial_number]['frequency'] + response.response_stages[0].stage_gain = measured_sensitivity + response.response_stages[0].stage_gain_frequency = at_frequency + print(f'{station}: {nominal_sensitivity} --> {measured_sensitivity} V/Pa') + # KEY: Recalculate overall sensitivity (at the calibration frequency) + response.recalculate_overall_sensitivity(frequency=at_frequency) + + cha.response = response + sta.channels.append(cha) + net.stations.append(sta) + + # Set network start and end dates from station info + net.start_date = min([sta.start_date for sta in net]) + net.end_date = max([sta.end_date for sta in net]) + + inv.networks.append(net) + + # Write to StationXML file + output_filename = Path(output_filename.rstrip('.xml') + '.xml').absolute() + inv.write(output_filename, format='stationxml', validate=True) + print(f'Wrote StationXML file to {output_filename}\n') # Run the main function if this is called as a script From 328eaebae8d02d3b4b021fb15157a055a90dc807 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 17:41:58 -0700 Subject: [PATCH 07/26] Minor docs updates --- cube_conversion/cube_stationxml.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 6128c5d..6e31cb2 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -66,7 +66,7 @@ def main(): help='directory containing miniSEED files and coordinate files produced by cube_convert', ) parser.add_argument( - 'station_mappings', + 'station_mapping', nargs='+', help='one or more mappings of the form STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2', ) @@ -76,7 +76,7 @@ def main(): parser.add_argument( '--nrl-path', default=None, - help='path to local copy of the NRL (if not provided, use online NRL)', + help='path to local copy of the NRL (if not provided, uses online NRL)', ) input_args = parser.parse_args() @@ -87,7 +87,7 @@ def main(): # Parse mappings into dictionary with keys being the station codes (TODO: validate!) station_mappings = {} - for mapping in input_args.station_mappings: + for mapping in input_args.station_mapping: station_code, cube_name, sensor_serial = mapping.split(':') station_mappings[station_code] = { 'cube_name': cube_name, From 5bae49a82486f00c979deffb2ff28d1bcf6d7a75 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 17:46:41 -0700 Subject: [PATCH 08/26] Fix KeyError --- cube_conversion/cube_stationxml.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 6e31cb2..cfe33d3 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -220,8 +220,8 @@ def main(): datalogger_keys=[ _DATALOGGER_MANUFACTURER[option_ind], _DATALOGGER_MODEL[option_ind], - str(GAIN), - str(sample_rate), + f'{GAIN:g}', + f'{sample_rate:g}', ], ) # KEY: Add a response stage which applies the breakout box factor, after first From 87d52c6267bbc02291c4fb8118990268d449c0af Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 17:47:55 -0700 Subject: [PATCH 09/26] Fix filepath def --- cube_conversion/cube_stationxml.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index cfe33d3..7c955b1 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -265,7 +265,9 @@ def main(): inv.networks.append(net) # Write to StationXML file - output_filename = Path(output_filename.rstrip('.xml') + '.xml').absolute() + output_filename = Path( + input_args.output_filename.rstrip('.xml') + '.xml' + ).absolute() inv.write(output_filename, format='stationxml', validate=True) print(f'Wrote StationXML file to {output_filename}\n') From 3c2bea9d8e5b38b7ec18ffd0ceee8c2b26f3dcad Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 17:50:50 -0700 Subject: [PATCH 10/26] More messages --- cube_conversion/cube_stationxml.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 7c955b1..a201a05 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -247,12 +247,14 @@ def main(): # KEY: Update with our specific measured sensor sensitivity nominal_sensitivity = response.response_stages[0].stage_gain measured_sensitivity = sensor_info[serial_number]['sensitivity'] - at_frequency = sensor_info[serial_number]['frequency'] + nominal_frequency = response.response_stages[0].stage_gain_frequency + measured_frequency = sensor_info[serial_number]['frequency'] response.response_stages[0].stage_gain = measured_sensitivity - response.response_stages[0].stage_gain_frequency = at_frequency + response.response_stages[0].stage_gain_frequency = measured_frequency print(f'{station}: {nominal_sensitivity} --> {measured_sensitivity} V/Pa') + print(f'{station}: {nominal_frequency:.2f} --> {measured_frequency:.2f} Hz') # KEY: Recalculate overall sensitivity (at the calibration frequency) - response.recalculate_overall_sensitivity(frequency=at_frequency) + response.recalculate_overall_sensitivity(frequency=measured_frequency) cha.response = response sta.channels.append(cha) From 0f1741d5b19f9e40929ab22d23cb323d1c91f303 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 18:22:29 -0700 Subject: [PATCH 11/26] Sigh... handling multiple NRL versions --- cube_conversion/cube_stationxml.py | 46 ++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 12 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index a201a05..76a510d 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -204,26 +204,48 @@ def main(): # Access the NRL to get response information. If the user provided a local path # to the NRL, use that; otherwise, use the online NRL. if nrl_path is not None: - nrl = NRL(nrl_path) - option_ind = 0 # Use v2 NRL entries + nrl = NRL(str(nrl_path)) + lp_corner = list(nrl.sensors[_SENSOR_MANUFACTURER[0]][sensor_model]) + assert len(lp_corner) == 1, 'Multiple low-pass corner options found!' + lp_corner = lp_corner[0] + hf_corner = list( + nrl.sensors[_SENSOR_MANUFACTURER[0]][sensor_model][lp_corner] + ) + assert len(hf_corner) == 1, 'Multiple high-pass corner options found!' + hf_corner = hf_corner[0] + sensor_keys = [_SENSOR_MANUFACTURER[0], sensor_model, lp_corner, hf_corner] + datalogger_keys = [ + _DATALOGGER_MANUFACTURER[0], + _DATALOGGER_MODEL[0], + f'{GAIN:g}', + f'{sample_rate:g} Hz', + ] else: nrl = NRL() - option_ind = 1 # Use v1 NRL entries + sensor_keys = [_SENSOR_MANUFACTURER[1], sensor_model] + datalogger_keys = [ + _DATALOGGER_MANUFACTURER[1], + _DATALOGGER_MODEL[1], + f'{GAIN:g}', + f'{sample_rate:g}', + ] # The contents of the NRL can be explored interactively in a Python prompt, see # API documentation of NRL submodule: # http://docs.obspy.org/packages/obspy.clients.nrl.html # Here we assume that the end point of data logger and sensor are already known. # Get the nominal response for this combination of sensor and digitizer - response = nrl.get_response( - sensor_keys=[_SENSOR_MANUFACTURER[option_ind], sensor_model], - datalogger_keys=[ - _DATALOGGER_MANUFACTURER[option_ind], - _DATALOGGER_MODEL[option_ind], - f'{GAIN:g}', - f'{sample_rate:g}', - ], - ) + try: + response = nrl.get_response( + sensor_keys=sensor_keys, datalogger_keys=datalogger_keys + ) + except KeyError as e: + msg = ( + f'Could not find response in NRL for sensor keys {sensor_keys} and ' + f'datalogger keys {datalogger_keys}.' + f'\n\n{nrl.sensors[sensor_keys[0]]}' + ) + raise Exception(msg) from e # KEY: Add a response stage which applies the breakout box factor, after first # stage (i.e., right after Pa --> V) — this should be a PZ stage, see: # https://docs.fdsn.org/projects/stationxml/en/latest/reference.html#response-stage From eec6eacadda85cdb4b317151bf45f5a40d0c89ee Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 18:27:13 -0700 Subject: [PATCH 12/26] Add note on multi-channel setups --- cube_conversion/cube_stationxml.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 76a510d..3566982 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -191,13 +191,13 @@ def main(): cha = Channel( code=channel_code, location_code=location_code, - latitude=latitude, - longitude=longitude, + latitude=latitude, # Could be different if multi-channel setup! + longitude=longitude, # Could be different if multi-channel setup! elevation=elevation, depth=0, # Required, always 0? sample_rate=sample_rate, - start_date=station_starttime, - end_date=station_endtime, + start_date=station_starttime, # Could be different if multi-channel setup! + end_date=station_endtime, # Could be different if multi-channel setup! sensor=sensor, data_logger=data_logger, ) From c66187afc6bb10be32877ef1eefc27d9aaa4713a Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 18:36:25 -0700 Subject: [PATCH 13/26] Ignore some warnings --- cube_conversion/cube_stationxml.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 3566982..f01832f 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -21,6 +21,7 @@ import argparse import json +import warnings from pathlib import Path from obspy import Stream, read @@ -34,6 +35,7 @@ Site, Station, ) +from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning # -------------------------------------------------------------------------------------- # Advanced configuration options @@ -52,6 +54,15 @@ _DATALOGGER_MANUFACTURER = 'DiGOSOmnirecs', 'DiGOS/Omnirecs' _DATALOGGER_MODEL = 'DataCube', 'DATACUBE' +# Ignore warning about Pa units +warnings.filterwarnings( + 'ignore', + category=UserWarning, + message="ObsPy can not map unit 'PA' to displacement, velocity, or acceleration", +) +# Ignore deprecation warning about accessing the NRL +warnings.filterwarnings('ignore', category=ObsPyDeprecationWarning) + # Define callable main function to work with [project.scripts] def main(): @@ -299,3 +310,4 @@ def main(): # Run the main function if this is called as a script if __name__ == '__main__': main() + main() From 1627e369a1130ccad2d5cfa7c4286326f752a10d Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 18:52:21 -0700 Subject: [PATCH 14/26] Better printing --- cube_conversion/cube_stationxml.py | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index f01832f..b365aaf 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -142,6 +142,18 @@ def main(): inv = Inventory() net = Network(code=network_code) for station in station_codes: + + # Get info about this station from mappings and sensor info + serial_number = station_mappings[station]['sensor_serial'] + cube_name = station_mappings[station]['cube_name'] + sensor_model = sensor_info[serial_number]['model'] + print('\n------------------------------') + print(f'{network_code}.{station}') + print('------------------------------') + print(f'Cube name: {cube_name}') + print(f'Sensor serial number: {serial_number}') + print(f'Sensor model: {sensor_model}') + # Get station start and end times, SEED channel code, sample rate... st_station = st.select(station=station) station_starttime = min(tr.stats.starttime for tr in st_station) @@ -181,8 +193,6 @@ def main(): end_date=station_endtime, ) # Define sensor Equipment object - serial_number = station_mappings[station]['sensor_serial'] - sensor_model = sensor_info[serial_number]['model'] sensor = Equipment( type='Infrasound sensor', description=f'{_SENSOR_MANUFACTURER[0]} {sensor_model}', # MDA shows this! @@ -191,7 +201,6 @@ def main(): serial_number=serial_number, ) # Define digitizer Equipment object - cube_name = station_mappings[station]['cube_name'] data_logger = Equipment( type='Digitizer', manufacturer=_DATALOGGER_MANUFACTURER[0], @@ -284,8 +293,9 @@ def main(): measured_frequency = sensor_info[serial_number]['frequency'] response.response_stages[0].stage_gain = measured_sensitivity response.response_stages[0].stage_gain_frequency = measured_frequency - print(f'{station}: {nominal_sensitivity} --> {measured_sensitivity} V/Pa') - print(f'{station}: {nominal_frequency:.2f} --> {measured_frequency:.2f} Hz') + print('Updated sensor sensitivity:') + print(f'\t{nominal_sensitivity} --> {measured_sensitivity} V/Pa') + print(f'\t{nominal_frequency:.2f} --> {measured_frequency:.2f} Hz') # KEY: Recalculate overall sensitivity (at the calibration frequency) response.recalculate_overall_sensitivity(frequency=measured_frequency) @@ -304,7 +314,7 @@ def main(): input_args.output_filename.rstrip('.xml') + '.xml' ).absolute() inv.write(output_filename, format='stationxml', validate=True) - print(f'Wrote StationXML file to {output_filename}\n') + print(f'\nWrote StationXML file to {output_filename}\n') # Run the main function if this is called as a script From 3d3339fe3d4a65031a37ec8f5723d152ef50ec13 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 19:03:13 -0700 Subject: [PATCH 15/26] Add StationXML validator from IRIS --- .gitignore | 3 ++ cube_conversion/cube_stationxml.py | 46 +++++++++++++++++++++++++++++- 2 files changed, 48 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 0f95b42..e2d826e 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,6 @@ # Directory that gets created when installing /cube_conversion.egg-info/ + +# IRIS StationXML validator which may be downloaded +/stationxml-validator-*.jar diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index b365aaf..8b63983 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -21,8 +21,10 @@ import argparse import json +import subprocess import warnings from pathlib import Path +from urllib.request import urlretrieve from obspy import Stream, read from obspy.clients.nrl import NRL @@ -89,6 +91,11 @@ def main(): default=None, help='path to local copy of the NRL (if not provided, uses online NRL)', ) + parser.add_argument( + '--validate', + action='store_true', + help='run the IRIS StationXML validator on the output file', + ) input_args = parser.parse_args() # Check if input directory is valid @@ -316,8 +323,45 @@ def main(): inv.write(output_filename, format='stationxml', validate=True) print(f'\nWrote StationXML file to {output_filename}\n') + # Download and run the StationXML validator, if user wants to + if input_args.validate: + # Functions for colored printing... helpful for highlighting errors in validator + # output! + # fmt: off + def print_red(string): + print('\u001b[31m' + string + '\u001b[0m') + def print_yellow(string): + print('\u001b[33m' + string + '\u001b[0m') + # fmt: on + # JAR file will downloaded here if it doesn't already exist + jar_file_path = root_dir / 'stationxml-validator-1.7.5.jar' # Selects version! + if jar_file_path.is_file(): + print(f'Found {jar_file_path.name}. Running command:') + else: + print(f'Downloading {jar_file_path.name}...') + url_base = 'https://github.com/iris-edu/stationxml-validator/releases/' + urlretrieve( + url_base + f'download/{jar_file_path.stem}/{jar_file_path.name}', + jar_file_path, + ) + print('...done. Running command:') + args = ['java', '-jar', jar_file_path, output_filename] + print('\n' + ' '.join([str(arg) for arg in args]) + '\n') + process = subprocess.run(args, capture_output=True, text=True) + if process.stderr: # Something went wrong! + print_red(process.stderr.strip()) + elif process.stdout: # Normal operation + for line in process.stdout.strip().split('\n'): + if ',Error,' in line: + print_red(line) + elif ',Warning,' in line: + print_yellow(line) + else: + print(line) + else: # This will never happen? + raise OSError + # Run the main function if this is called as a script if __name__ == '__main__': main() - main() From 485e016c7a5629c54367b477b26de4b22a96bef4 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 19:09:09 -0700 Subject: [PATCH 16/26] Tweak message printing --- cube_conversion/cube_stationxml.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 8b63983..9069596 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -346,7 +346,10 @@ def print_yellow(string): ) print('...done. Running command:') args = ['java', '-jar', jar_file_path, output_filename] - print('\n' + ' '.join([str(arg) for arg in args]) + '\n') + print('\t' + ' '.join([str(arg) for arg in args])) + print('------------------------------') + print('Output from validator') + print('------------------------------') process = subprocess.run(args, capture_output=True, text=True) if process.stderr: # Something went wrong! print_red(process.stderr.strip()) From fc60f05276aaee2de1c537ad9250c3b690394985 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 19:33:52 -0700 Subject: [PATCH 17/26] Big README update --- README.md | 71 +++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 56 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index 4afe788..e0be7a0 100644 --- a/README.md +++ b/README.md @@ -1,24 +1,25 @@ cube_conversion =============== -This command-line tool converts [DiGOS](https://digos.eu/) DATA-CUBE3 -files into miniSEED files of a desired length of time with specified metadata. -Output miniSEED files have units of Pa, unless the user selects to export the files in -a form suitable for submission to EarthScope (formerly IRIS). A script to produce -StationXML files is also provided, which is useful for EarthScope uploads. The tool -can differentiate between channels for 3 channel DATA-CUBE3 files and -optionally extract coordinates from the digitizer's GPS. The code only looks for -files from digitizers defined in the `digitizer_sensor_pairs.json` file. Therefore, -this file must be updated if pairings change or new pairings are added. The user -can specify a custom "breakout box factor" for setups that modify the signal -voltage via a voltage divider. This tool is currently only set up for conversion -of infrasound data, but future updates will accommodate seismic as well. +These command-line tools convert [DiGOS](https://digos.eu/) DATA-CUBE³ files into +miniSEED files of a desired length of time with specified metadata, and produce +validated StationXML files containing deployment, sensor, and digitizer information +including sensor responses adjusted for individual sensor sensitivities. Output miniSEED +files have units of Pa, unless the user selects to export the files in a form suitable +for submission to EarthScope (formerly IRIS) — namely, integer counts. The miniSEED +conversion tool can differentiate between channels for 3 channel DATA-CUBE³ files and +optionally extract coordinates from the digitizer's GPS. The code only looks for files +from digitizers defined in the `digitizer_sensor_pairs.json` file. Therefore, this file +must be updated if pairings change or new pairings are added. The user can specify a +custom "breakout box factor" for setups that modify the signal voltage via a voltage +divider. This tool is currently only set up for conversion of infrasound data, but +future updates could accommodate seismic as well. Installation ------------ It's recommended that you run this script within a new or pre-existing -[conda](https://docs.conda.io/projects/conda/en/latest/index.html) environment. +[conda](https://docs.conda.io/projects/conda/en/stable/) environment. (If you choose the latter option, ensure that your environment contains all of the packages listed in the [Dependencies](#dependencies) section.) @@ -79,9 +80,15 @@ Supplemental files Usage ----- -To print the conversion script's help menu, execute the following terminal commands: +To run either command-line tool first activate your conda environment, e.g., ``` conda activate uafinfra # Or your pre-existing env +``` + +## `cube_convert` + +To print the conversion script's help menu, run: +``` cube_convert --help # Or: python /path/to/cube_conversion/cube_convert.py --help ``` The help menu is shown below. @@ -124,6 +131,40 @@ means "convert all files in the subdirectories of `~/data/raw/` and place in **01**, and an automatically determined channel code, dividing the sensitivity by 4.5 and extracting coordinates from the digitizer's GPS." +## `cube_stationxml` + +To print the StationXML generation script's help menu, run: +``` +cube_stationxml --help # Or: python /path/to/cube_conversion/cube_stationxml.py --help +``` +The help menu is shown below. +``` +usage: cube_stationxml [-h] [--nrl-path NRL_PATH] [--validate] + input_dir station_mapping [station_mapping ...] output_filename + +Generate StationXML files from DATA-CUBE³ miniSEED files and metadata. + +positional arguments: + input_dir directory containing miniSEED files and coordinate files produced + by cube_convert + station_mapping one or more mappings of the form + STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2 + output_filename filename for the output StationXML file (full path) + +options: + -h, --help show this help message and exit + --nrl-path NRL_PATH path to local copy of the NRL (if not provided, uses online NRL) + --validate run the IRIS StationXML validator on the output file +``` +For example, the command +``` +cube_stationxml ~/data/miniseed/ UAF1:B52:903V2 ~/data/station.xml --validate +``` +means "process all miniSEED files in `~/data/miniseed/` assuming that station **UAF1** +consisted of DATA-CUBE³ **B52** connected to sensor serial number **903V2**, and output +a StationXML file to `~/data/station.xml`, validating the output file with the IRIS +StationXML validator." + A note on SEED band codes ------------------------- @@ -149,7 +190,7 @@ letter of the channel code; e.g., the "B" in "BDF". | R | Extremely long-period | ≥ 0.0001 to < 0.001 | | Note that the band code depends on both the sample rate of the digitizer and -the corner period of the sensor. In `cube_convert.py` we allow for "B", "H", or +the corner period of the sensor. In `cube_convert` we allow for "B", "H", or "C", which covers a range of sample rates from 10 to 1000 Hz, all for corner periods of 10 s or greater. While this covers most infrasound sensors, please confirm that your digitizer sample rate and sensor corner period fit into the From 6bd6b7bf90644731ec36f2ed650a4e15d00cbdd3 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 19:35:35 -0700 Subject: [PATCH 18/26] Adjust width of cube_convert menu to match cube_stationxml --- README.md | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index e0be7a0..8044d37 100644 --- a/README.md +++ b/README.md @@ -93,33 +93,32 @@ cube_convert --help # Or: python /path/to/cube_conversion/cube_convert.py --hel ``` The help menu is shown below. ``` -usage: cube_convert [-h] [-v] [--grab-gps] - [--bob-factor BREAKOUT_BOX_FACTOR] [--earthscope] +usage: cube_convert [-h] [-v] [--grab-gps] [--bob-factor BREAKOUT_BOX_FACTOR] + [--earthscope] input_dir [input_dir ...] output_dir network station {01,02,03,04,AUTO} {AUTO,BDF,HDF,CDF} -Convert DATA-CUBE files to miniSEED files while trimming, adding metadata, and -renaming. Optionally extract coordinates from digitizer GPS. +Convert DATA-CUBE files to miniSEED files while trimming, adding metadata, and renaming. +Optionally extract coordinates from digitizer GPS. positional arguments: - input_dir one or more directories containing raw DATA-CUBE files - (all files must originate from a single digitizer) - [wildcards (*) supported] + input_dir one or more directories containing raw DATA-CUBE files (all files + must originate from a single digitizer) [wildcards (*) supported] output_dir directory for output miniSEED and GPS-related files network desired SEED network code (2 characters, A-Z) station desired SEED station code (3-5 characters, A-Z & 0-9) - {01,02,03,04,AUTO} desired SEED location code (if AUTO, choose - automatically for 3 channel DATA-CUBE files) - {AUTO,BDF,HDF,CDF} desired SEED channel code (if AUTO, determine - automatically using SEED convention [preferred]) + {01,02,03,04,AUTO} desired SEED location code (if AUTO, choose automatically for 3 + channel DATA-CUBE files) + {AUTO,BDF,HDF,CDF} desired SEED channel code (if AUTO, determine automatically using + SEED convention [preferred]) options: -h, --help show this help message and exit -v, --verbose enable verbosity for GIPPtools commands --grab-gps additionally extract coordinates from digitizer GPS --bob-factor BREAKOUT_BOX_FACTOR - factor by which to divide sensitivity values (for - custom breakout boxes [4.5 for UAF DATA-CUBEs]) + factor by which to divide sensitivity values (for custom breakout + boxes [4.5 for UAF DATA-CUBEs]) --earthscope format miniSEED files for EarthScope (formerly IRIS) data upload ``` For example, the command From 7f7d6cfa79ad455130d9e2aa07e8404a1fdcad40 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Wed, 14 Jan 2026 19:55:59 -0700 Subject: [PATCH 19/26] Add note on ObsPy version --- README.md | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 8044d37..136d647 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,14 @@ packages, execute the following terminal command: conda create -n uafinfra -c conda-forge obspy ``` This creates a new environment called `uafinfra` with ObsPy and its dependencies -installed. +installed. **Note:** if you intend to run the StationXML generation tool +(`cube_stationxml`), you must install development ObsPy (after activating your target +environment) via: +``` +pip install git+https://github.com/obspy/obspy.git +``` +Once ObsPy 1.5.0 is released, this step will no longer be necessary. See the docstring +at the top of `cube_stationxml.py` for more information. You must also install GIPPtools and add it to your path as described in the [Dependencies](#dependencies) section. From 30989a07afe9b68894aa46862628e6bacf2d7853 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Fri, 16 Jan 2026 11:46:12 -0700 Subject: [PATCH 20/26] KEY: Remove option to access old NRL since it is not updated! --- README.md | 31 ++++++----- cube_conversion/cube_stationxml.py | 83 +++++++++++------------------- 2 files changed, 49 insertions(+), 65 deletions(-) diff --git a/README.md b/README.md index 136d647..451aa68 100644 --- a/README.md +++ b/README.md @@ -84,6 +84,12 @@ Supplemental files * `sensor_sensitivities.json` — Infrasound sensor model, sensitivity in V/Pa, and frequency of sensitivity measurement in Hz +To use `cube_stationxml`, you also must have a local copy of the Nominal Response +Library (NRL) downloaded. It can be obtained by following the instructions +[here](https://ds.iris.edu/ds/nrl/). Be sure to unzip the downloaded file to create an +`NRL/` directory. **Note:** Obviously, since this is offline, you may need to update +your local copy if a key response change is made, or if a new sensor is added! + Usage ----- @@ -145,31 +151,32 @@ cube_stationxml --help # Or: python /path/to/cube_conversion/cube_stationxml.py ``` The help menu is shown below. ``` -usage: cube_stationxml [-h] [--nrl-path NRL_PATH] [--validate] +usage: cube_stationxml [-h] [--validate] input_dir station_mapping [station_mapping ...] output_filename + nrl_path Generate StationXML files from DATA-CUBE³ miniSEED files and metadata. positional arguments: - input_dir directory containing miniSEED files and coordinate files produced - by cube_convert - station_mapping one or more mappings of the form - STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2 - output_filename filename for the output StationXML file (full path) + input_dir directory containing miniSEED files and coordinate files produced + by cube_convert + station_mapping one or more mappings of the form + STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2 + output_filename filename for the output StationXML file (full path) + nrl_path path to local copy of the NRL (Nominal Response Library) directory options: - -h, --help show this help message and exit - --nrl-path NRL_PATH path to local copy of the NRL (if not provided, uses online NRL) - --validate run the IRIS StationXML validator on the output file + -h, --help show this help message and exit + --validate run the IRIS StationXML validator on the output file ``` For example, the command ``` -cube_stationxml ~/data/miniseed/ UAF1:B52:903V2 ~/data/station.xml --validate +cube_stationxml ~/data/miniseed/ UAF1:B52:903V2 ~/data/station.xml ~/NRL/ --validate ``` means "process all miniSEED files in `~/data/miniseed/` assuming that station **UAF1** consisted of DATA-CUBE³ **B52** connected to sensor serial number **903V2**, and output -a StationXML file to `~/data/station.xml`, validating the output file with the IRIS -StationXML validator." +a StationXML file to `~/data/station.xml`, using the local copy of the NRL at `~/NRL/`, +and validating the output file with the IRIS StationXML validator." A note on SEED band codes ------------------------- diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 9069596..00ce0bd 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -37,7 +37,6 @@ Site, Station, ) -from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning # -------------------------------------------------------------------------------------- # Advanced configuration options @@ -46,15 +45,11 @@ BOB_FACTOR = 10 # Breakout box factor for DATA-CUBE³s (this should usually be 10!) # -------------------------------------------------------------------------------------- -# These names must match EXACTLY what is in the NRL. Unforunately, there are differences -# between the NRL accessed through ObsPy (via online) versus a local copy of the NRL. -# These are likely related to NRL v1 versus v2, but it's not clear. I think the v2 ones -# should be preferred, and they can be viewed here: https://ds.iris.edu/ds/nrl/ -# These show up first in the tuples below, the second entries are for the older NRL v1 -# and we include those just to cover the bases. -_SENSOR_MANUFACTURER = 'Chaparral', 'Chaparral Physics' -_DATALOGGER_MANUFACTURER = 'DiGOSOmnirecs', 'DiGOS/Omnirecs' -_DATALOGGER_MODEL = 'DataCube', 'DATACUBE' +# These names must match EXACTLY what is in the NRL (v2). See here: +# https://ds.iris.edu/ds/nrl/ +_SENSOR_MANUFACTURER = 'Chaparral' +_DATALOGGER_MANUFACTURER = 'DiGOSOmnirecs' +_DATALOGGER_MODEL = 'DataCube' # Ignore warning about Pa units warnings.filterwarnings( @@ -62,8 +57,6 @@ category=UserWarning, message="ObsPy can not map unit 'PA' to displacement, velocity, or acceleration", ) -# Ignore deprecation warning about accessing the NRL -warnings.filterwarnings('ignore', category=ObsPyDeprecationWarning) # Define callable main function to work with [project.scripts] @@ -87,9 +80,8 @@ def main(): 'output_filename', help='filename for the output StationXML file (full path)' ) parser.add_argument( - '--nrl-path', - default=None, - help='path to local copy of the NRL (if not provided, uses online NRL)', + 'nrl_path', + help='path to local copy of the NRL (Nominal Response Library) directory', ) parser.add_argument( '--validate', @@ -112,12 +104,10 @@ def main(): 'sensor_serial': sensor_serial, } - # Check if NRL path is a directory if provided - nrl_path = input_args.nrl_path - if nrl_path is not None: - nrl_path = Path(nrl_path).absolute() - if not nrl_path.is_dir(): - raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') + # Check if NRL path is a directory + nrl_path = Path(input_args.nrl_path).expanduser().absolute() + if not nrl_path.is_dir(): + raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') # Find root directory for cube_conversion repo root_dir = Path(__file__).parents[1] @@ -202,16 +192,16 @@ def main(): # Define sensor Equipment object sensor = Equipment( type='Infrasound sensor', - description=f'{_SENSOR_MANUFACTURER[0]} {sensor_model}', # MDA shows this! - manufacturer=_SENSOR_MANUFACTURER[0], + description=f'{_SENSOR_MANUFACTURER} {sensor_model}', # MDA shows this! + manufacturer=_SENSOR_MANUFACTURER, model=sensor_model, serial_number=serial_number, ) # Define digitizer Equipment object data_logger = Equipment( type='Digitizer', - manufacturer=_DATALOGGER_MANUFACTURER[0], - model=_DATALOGGER_MODEL[0], + manufacturer=_DATALOGGER_MANUFACTURER, + model=_DATALOGGER_MODEL, serial_number=cube_name, ) # Make Channel object @@ -228,34 +218,21 @@ def main(): sensor=sensor, data_logger=data_logger, ) - # Access the NRL to get response information. If the user provided a local path - # to the NRL, use that; otherwise, use the online NRL. - if nrl_path is not None: - nrl = NRL(str(nrl_path)) - lp_corner = list(nrl.sensors[_SENSOR_MANUFACTURER[0]][sensor_model]) - assert len(lp_corner) == 1, 'Multiple low-pass corner options found!' - lp_corner = lp_corner[0] - hf_corner = list( - nrl.sensors[_SENSOR_MANUFACTURER[0]][sensor_model][lp_corner] - ) - assert len(hf_corner) == 1, 'Multiple high-pass corner options found!' - hf_corner = hf_corner[0] - sensor_keys = [_SENSOR_MANUFACTURER[0], sensor_model, lp_corner, hf_corner] - datalogger_keys = [ - _DATALOGGER_MANUFACTURER[0], - _DATALOGGER_MODEL[0], - f'{GAIN:g}', - f'{sample_rate:g} Hz', - ] - else: - nrl = NRL() - sensor_keys = [_SENSOR_MANUFACTURER[1], sensor_model] - datalogger_keys = [ - _DATALOGGER_MANUFACTURER[1], - _DATALOGGER_MODEL[1], - f'{GAIN:g}', - f'{sample_rate:g}', - ] + # Access the local NRL to get response information + nrl = NRL(str(nrl_path)) + lp_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model]) + assert len(lp_corner) == 1, 'Multiple low-pass corner options found!' + lp_corner = lp_corner[0] + hf_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model][lp_corner]) + assert len(hf_corner) == 1, 'Multiple high-pass corner options found!' + hf_corner = hf_corner[0] + sensor_keys = [_SENSOR_MANUFACTURER, sensor_model, lp_corner, hf_corner] + datalogger_keys = [ + _DATALOGGER_MANUFACTURER, + _DATALOGGER_MODEL, + f'{GAIN:g}', + f'{sample_rate:g} Hz', + ] # The contents of the NRL can be explored interactively in a Python prompt, see # API documentation of NRL submodule: # http://docs.obspy.org/packages/obspy.clients.nrl.html From e5443d9e799ddba57b92baf0fd807da376a78346 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Fri, 16 Jan 2026 11:46:56 -0700 Subject: [PATCH 21/26] Use expanduser for StationXML path too --- cube_conversion/cube_stationxml.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 00ce0bd..984d729 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -294,9 +294,9 @@ def main(): inv.networks.append(net) # Write to StationXML file - output_filename = Path( - input_args.output_filename.rstrip('.xml') + '.xml' - ).absolute() + output_filename = ( + Path(input_args.output_filename.rstrip('.xml') + '.xml').expanduser().absolute() + ) inv.write(output_filename, format='stationxml', validate=True) print(f'\nWrote StationXML file to {output_filename}\n') From bed3eeea5d885c4c8e6acea94bec58669e5f647e Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Fri, 16 Jan 2026 11:49:32 -0700 Subject: [PATCH 22/26] Print location being used for local NRL --- cube_conversion/cube_stationxml.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 984d729..6bde976 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -21,6 +21,7 @@ import argparse import json +import os import subprocess import warnings from pathlib import Path @@ -108,6 +109,7 @@ def main(): nrl_path = Path(input_args.nrl_path).expanduser().absolute() if not nrl_path.is_dir(): raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') + print(f'Using local copy of NRL at {nrl_path}{os.sep}') # Find root directory for cube_conversion repo root_dir = Path(__file__).parents[1] From da3d8339a189e9c329a8d505622a0eaae398b111 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Fri, 16 Jan 2026 11:53:21 -0700 Subject: [PATCH 23/26] Format messages using backticks where appropriate --- cube_conversion/cube_stationxml.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 6bde976..8843729 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -109,7 +109,7 @@ def main(): nrl_path = Path(input_args.nrl_path).expanduser().absolute() if not nrl_path.is_dir(): raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') - print(f'Using local copy of NRL at {nrl_path}{os.sep}') + print(f'Using local copy of NRL at `{nrl_path}{os.sep}`') # Find root directory for cube_conversion repo root_dir = Path(__file__).parents[1] @@ -300,7 +300,7 @@ def main(): Path(input_args.output_filename.rstrip('.xml') + '.xml').expanduser().absolute() ) inv.write(output_filename, format='stationxml', validate=True) - print(f'\nWrote StationXML file to {output_filename}\n') + print(f'\nWrote StationXML file to `{output_filename}`\n') # Download and run the StationXML validator, if user wants to if input_args.validate: @@ -315,9 +315,9 @@ def print_yellow(string): # JAR file will downloaded here if it doesn't already exist jar_file_path = root_dir / 'stationxml-validator-1.7.5.jar' # Selects version! if jar_file_path.is_file(): - print(f'Found {jar_file_path.name}. Running command:') + print(f'Found `{jar_file_path.name}` — running command:') else: - print(f'Downloading {jar_file_path.name}...') + print(f'Downloading `{jar_file_path.name}`...') url_base = 'https://github.com/iris-edu/stationxml-validator/releases/' urlretrieve( url_base + f'download/{jar_file_path.stem}/{jar_file_path.name}', @@ -325,7 +325,9 @@ def print_yellow(string): ) print('...done. Running command:') args = ['java', '-jar', jar_file_path, output_filename] - print('\t' + ' '.join([str(arg) for arg in args])) + print('```') + print(' '.join([str(arg) for arg in args])) + print('```') print('------------------------------') print('Output from validator') print('------------------------------') From 2b828d74e9ad3f93ad4935cba28c78a42b755399 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Fri, 16 Jan 2026 15:01:13 -0700 Subject: [PATCH 24/26] Add option to use NRL web services, which is still NRL v2 --- README.md | 33 ++++---- cube_conversion/cube_stationxml.py | 116 +++++++++++++++++++---------- 2 files changed, 95 insertions(+), 54 deletions(-) diff --git a/README.md b/README.md index 451aa68..79daeed 100644 --- a/README.md +++ b/README.md @@ -84,11 +84,12 @@ Supplemental files * `sensor_sensitivities.json` — Infrasound sensor model, sensitivity in V/Pa, and frequency of sensitivity measurement in Hz -To use `cube_stationxml`, you also must have a local copy of the Nominal Response -Library (NRL) downloaded. It can be obtained by following the instructions +For `cube_stationxml`, you may optionally use a local copy of the Nominal Response +Library (NRL). It can be obtained by following the instructions [here](https://ds.iris.edu/ds/nrl/). Be sure to unzip the downloaded file to create an `NRL/` directory. **Note:** Obviously, since this is offline, you may need to update -your local copy if a key response change is made, or if a new sensor is added! +your local copy if a key response change is made, or if a new sensor is added! Hence, +we recommend only using this local copy option if you're working offline. Usage ----- @@ -151,32 +152,32 @@ cube_stationxml --help # Or: python /path/to/cube_conversion/cube_stationxml.py ``` The help menu is shown below. ``` -usage: cube_stationxml [-h] [--validate] +usage: cube_stationxml [-h] [--nrl-path NRL_PATH] [--validate] input_dir station_mapping [station_mapping ...] output_filename - nrl_path Generate StationXML files from DATA-CUBE³ miniSEED files and metadata. positional arguments: - input_dir directory containing miniSEED files and coordinate files produced - by cube_convert - station_mapping one or more mappings of the form - STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2 - output_filename filename for the output StationXML file (full path) - nrl_path path to local copy of the NRL (Nominal Response Library) directory + input_dir directory containing miniSEED files and coordinate files produced + by cube_convert + station_mapping one or more mappings of the form + STATION_CODE:CUBE_NAME:SENSOR_SERIAL, for example UAF1:AVJ:903V2 + output_filename filename for the output StationXML file (full path) options: - -h, --help show this help message and exit - --validate run the IRIS StationXML validator on the output file + -h, --help show this help message and exit + --nrl-path NRL_PATH path to local copy of the NRL (Nominal Response Library) + directory, if not provided makes web services call + --validate run the IRIS StationXML validator on the output file ``` For example, the command ``` -cube_stationxml ~/data/miniseed/ UAF1:B52:903V2 ~/data/station.xml ~/NRL/ --validate +cube_stationxml ~/data/miniseed/ UAF1:B52:903V2 ~/data/station.xml --validate ``` means "process all miniSEED files in `~/data/miniseed/` assuming that station **UAF1** consisted of DATA-CUBE³ **B52** connected to sensor serial number **903V2**, and output -a StationXML file to `~/data/station.xml`, using the local copy of the NRL at `~/NRL/`, -and validating the output file with the IRIS StationXML validator." +a StationXML file to `~/data/station.xml`, validating the output file with the IRIS +StationXML validator." A note on SEED band codes ------------------------- diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 8843729..5d2ddf2 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -24,10 +24,12 @@ import os import subprocess import warnings +import xml.etree.ElementTree as ET from pathlib import Path from urllib.request import urlretrieve -from obspy import Stream, read +import requests +from obspy import Stream, read, read_inventory from obspy.clients.nrl import NRL from obspy.core.inventory import ( Channel, @@ -52,6 +54,9 @@ _DATALOGGER_MANUFACTURER = 'DiGOSOmnirecs' _DATALOGGER_MODEL = 'DataCube' +# Base url for all NRL web service queries +_BASE_URL = 'https://service.iris.edu/irisws/nrl/1/' + # Ignore warning about Pa units warnings.filterwarnings( 'ignore', @@ -81,8 +86,9 @@ def main(): 'output_filename', help='filename for the output StationXML file (full path)' ) parser.add_argument( - 'nrl_path', - help='path to local copy of the NRL (Nominal Response Library) directory', + '--nrl-path', + default=None, + help='path to local copy of the NRL (Nominal Response Library) directory, if not provided makes web services call', ) parser.add_argument( '--validate', @@ -105,11 +111,15 @@ def main(): 'sensor_serial': sensor_serial, } - # Check if NRL path is a directory - nrl_path = Path(input_args.nrl_path).expanduser().absolute() - if not nrl_path.is_dir(): - raise NotADirectoryError(f'NRL path {nrl_path} is not a directory!') - print(f'Using local copy of NRL at `{nrl_path}{os.sep}`') + # Check if NRL path is a directory if provided + nrl_path = input_args.nrl_path + if nrl_path is not None: + nrl_path = Path(input_args.nrl_path).expanduser().absolute() + if not nrl_path.is_dir(): + raise NotADirectoryError( + f'NRL path `{nrl_path}{os.sep}` is not a directory!' + ) + print(f'Using local copy of NRL at `{nrl_path}{os.sep}`') # Find root directory for cube_conversion repo root_dir = Path(__file__).parents[1] @@ -220,38 +230,68 @@ def main(): sensor=sensor, data_logger=data_logger, ) - # Access the local NRL to get response information - nrl = NRL(str(nrl_path)) - lp_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model]) - assert len(lp_corner) == 1, 'Multiple low-pass corner options found!' - lp_corner = lp_corner[0] - hf_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model][lp_corner]) - assert len(hf_corner) == 1, 'Multiple high-pass corner options found!' - hf_corner = hf_corner[0] - sensor_keys = [_SENSOR_MANUFACTURER, sensor_model, lp_corner, hf_corner] - datalogger_keys = [ - _DATALOGGER_MANUFACTURER, - _DATALOGGER_MODEL, - f'{GAIN:g}', - f'{sample_rate:g} Hz', - ] - # The contents of the NRL can be explored interactively in a Python prompt, see - # API documentation of NRL submodule: - # http://docs.obspy.org/packages/obspy.clients.nrl.html - # Here we assume that the end point of data logger and sensor are already known. - - # Get the nominal response for this combination of sensor and digitizer - try: - response = nrl.get_response( - sensor_keys=sensor_keys, datalogger_keys=datalogger_keys + # Access the NRL to get response information. If the user provided a local path + # to the NRL, use that; otherwise, use the NRL web service (requires network + # connection). + if nrl_path is not None: + nrl = NRL(str(nrl_path)) + lp_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model]) + assert len(lp_corner) == 1, 'Multiple low-pass corner options found!' + lp_corner = lp_corner[0] + hf_corner = list(nrl.sensors[_SENSOR_MANUFACTURER][sensor_model][lp_corner]) + assert len(hf_corner) == 1, 'Multiple high-pass corner options found!' + hf_corner = hf_corner[0] + sensor_keys = [_SENSOR_MANUFACTURER, sensor_model, lp_corner, hf_corner] + datalogger_keys = [ + _DATALOGGER_MANUFACTURER, + _DATALOGGER_MODEL, + f'{GAIN:g}', + f'{sample_rate:g} Hz', + ] + # The contents of the NRL can be explored interactively in a Python prompt, see + # API documentation of NRL submodule: + # http://docs.obspy.org/packages/obspy.clients.nrl.html + # Here we assume that the end point of data logger and sensor are already known. + + # Get the nominal response for this combination of sensor and digitizer + try: + response = nrl.get_response( + sensor_keys=sensor_keys, datalogger_keys=datalogger_keys + ) + except KeyError as e: + msg = ( + f'Could not find response in NRL for sensor keys {sensor_keys} and ' + f'datalogger keys {datalogger_keys}.' + f'\n\n{nrl.sensors[sensor_keys[0]]}' + ) + raise Exception(msg) from e + else: # Use NRL web service + # (1) Locate the `instconfig` for the Chaparral Physics sensor + params_catalog = dict( + element='sensor', + manufacturer=_SENSOR_MANUFACTURER, + model=sensor_model, + format='xml', + level='configuration', + ) + response = requests.get(url=_BASE_URL + 'catalog', params=params_catalog) + root = ET.fromstring(response.text) + sensor_instconfigs = [element.text for element in root.iter('instconfig')] + assert len(sensor_instconfigs) == 1, 'Multiple instrument configs found!' + sensor_instconfig = sensor_instconfigs[0] + # (2) Define the `instconfig` for the DiGOS DATA-CUBE³ digitizer (easy!) + digitizer_instconfig = f'datalogger_{_DATALOGGER_MANUFACTURER}_{_DATALOGGER_MODEL}_PG{GAIN:g}_FR{sample_rate:g}' + # (3) Combine the two `instconfig` responses into a single StationXML response + params_combine = dict( + instconfig=':'.join([sensor_instconfig, digitizer_instconfig]), + format='stationxml', ) - except KeyError as e: - msg = ( - f'Could not find response in NRL for sensor keys {sensor_keys} and ' - f'datalogger keys {datalogger_keys}.' - f'\n\n{nrl.sensors[sensor_keys[0]]}' + request = requests.Request( + 'GET', url=_BASE_URL + 'combine', params=params_combine ) - raise Exception(msg) from e + stationxml_url = request.prepare().url + # (4) Read the combined StationXML response into ObsPy + response = read_inventory(stationxml_url)[0][0][0].response # KEY: Add a response stage which applies the breakout box factor, after first # stage (i.e., right after Pa --> V) — this should be a PZ stage, see: # https://docs.fdsn.org/projects/stationxml/en/latest/reference.html#response-stage From b99b7a48baa142dfeb3acb60efba93393f28a403 Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Tue, 20 Jan 2026 08:59:01 -0700 Subject: [PATCH 25/26] Update cube conversion script for new JSON format --- cube_conversion/cube_convert.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cube_conversion/cube_convert.py b/cube_conversion/cube_convert.py index e91def5..8aba73b 100755 --- a/cube_conversion/cube_convert.py +++ b/cube_conversion/cube_convert.py @@ -161,7 +161,7 @@ def main(): else: sensor = digitizer_sensor_pairs[digitizer] try: - sensitivity = sensitivities[sensor] + sensitivity = sensitivities[sensor]['sensitivity'] except KeyError: warnings.warn('No matching sensitivities. Using default of ' f'{DEFAULT_SENSITIVITY} V/Pa.') From 4bdd36a15044b437c9a16afb12854326420381de Mon Sep 17 00:00:00 2001 From: Liam Toney Date: Thu, 7 May 2026 17:07:21 -0600 Subject: [PATCH 26/26] Update docs now that ObsPy 1.5.0 is out --- README.md | 9 +-------- cube_conversion/cube_stationxml.py | 9 --------- 2 files changed, 1 insertion(+), 17 deletions(-) diff --git a/README.md b/README.md index 00e9325..4b5cad6 100644 --- a/README.md +++ b/README.md @@ -29,14 +29,7 @@ packages, execute the following terminal command: conda create -n uafinfra -c conda-forge obspy ``` This creates a new environment called `uafinfra` with ObsPy and its dependencies -installed. **Note:** if you intend to run the StationXML generation tool -(`cube_stationxml`), you must install development ObsPy (after activating your target -environment) via: -``` -pip install git+https://github.com/obspy/obspy.git -``` -Once ObsPy 1.5.0 is released, this step will no longer be necessary. See the docstring -at the top of `cube_stationxml.py` for more information. +installed. You must also install GIPPtools and add it to your path as described in the [Dependencies](#dependencies) section. diff --git a/cube_conversion/cube_stationxml.py b/cube_conversion/cube_stationxml.py index 5d2ddf2..c2e3bde 100755 --- a/cube_conversion/cube_stationxml.py +++ b/cube_conversion/cube_stationxml.py @@ -8,15 +8,6 @@ The code assumes that your sensors are all Chaparral Physics products, that your digitizers are all DiGOS DATA-CUBE³s, and — IMPORTANTLY — that you only have one sensor connected to each digitizer (i.e., one channel per station). - -Importantly, this code *will not work properly* unless you're using ObsPy after commit -107666b: - https://github.com/obspy/obspy/tree/107666b9659cbfdc98470a47c003772377551264 -(This commit introduced a change that allows for re-computation of the overall -sensitivity for non-standard units — like Pa.) Until the 1.5.0 release of ObsPy comes -out, you can simply update ObsPy via: - pip install git+https://github.com/obspy/obspy.git -This will fetch and install the newest copy from GitHub. """ import argparse