Add flag for EarthScope-ready miniSEED export

README.md

@@ -3,7 +3,8 @@ cube_conversion
 
 This command-line tool converts [DiGOS](https://digos.eu/) DATA-CUBE<sup>3</sup>
 files into miniSEED files of a desired length of time with specified metadata.
-Output miniSEED files are ready for IRIS upload and have units of Pa. The tool
+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
 can differentiate between channels for 3 channel DATA-CUBE<sup>3</sup> 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,
@@ -81,7 +82,7 @@ $ python /path/to/cube_convert.py --help
 The help menu is shown below.
 ```
 usage: cube_convert.py [-h] [-v] [--grab-gps]
-                       [--bob-factor BREAKOUT_BOX_FACTOR]
+                       [--bob-factor BREAKOUT_BOX_FACTOR] [--earthscope]
                        input_dir [input_dir ...] output_dir network station
                        {01,02,03,04,AUTO} {AUTO,BDF,HDF,CDF}
 
@@ -100,13 +101,14 @@ positional arguments:
   {AUTO,BDF,HDF,CDF}    desired SEED channel code (if AUTO, determine
                         automatically using SEED convention [preferred])
 
-optional arguments:
+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])
+  --earthscope          format miniSEED files for EarthScope (formerly IRIS) data upload
 ```
 For example, the command
 ```

cube_convert.py

@@ -60,6 +60,8 @@
                     dest='breakout_box_factor',
                     help='factor by which to divide sensitivity values (for '
                          'custom breakout boxes [4.5 for UAF DATA-CUBEs])')
+parser.add_argument('--earthscope', action='store_true', dest='earthscope',
+                    help='format miniSEED files for EarthScope (formerly IRIS) data upload')
 input_args = parser.parse_args()
 
 # Check if input directory/ies is/are valid
@@ -144,30 +146,39 @@
     offset = DEFAULT_OFFSET
 print(f'    Digitizer: {digitizer} (offset = {offset} V)')
 
-# Get sensor info and sensitivity
-sensor = digitizer_sensor_pairs[digitizer]
-try:
-    sensitivity = sensitivities[sensor]
-except KeyError:
-    warnings.warn('No matching sensitivities. Using default of '
-                  f'{DEFAULT_SENSITIVITY} V/Pa.')
-    sensitivity = DEFAULT_SENSITIVITY
-print(f'       Sensor: {sensor} (sensitivity = {sensitivity} V/Pa)')
-
-# Apply breakout box correction factor if provided
-if input_args.breakout_box_factor:
-    sensitivity = sensitivity / input_args.breakout_box_factor
-    print('       Dividing sensitivity by breakout box factor of '
-          f'{input_args.breakout_box_factor}')
+# Get sensor info and sensitivity (EarthScope submission requires raw counts, so we skip
+# this step if the user specified `--earthscope`)
+if input_args.earthscope:
+    # Remind user that breakout box correction factor is not applied to raw data if
+    # uploading data to EarthScope
+    if input_args.breakout_box_factor is not None:
+        warnings.warn('Breakout box correction factor ignored for EarthScope submission!')
+else:
+    sensor = digitizer_sensor_pairs[digitizer]
+    try:
+        sensitivity = sensitivities[sensor]
+    except KeyError:
+        warnings.warn('No matching sensitivities. Using default of '
+                      f'{DEFAULT_SENSITIVITY} V/Pa.')
+        sensitivity = DEFAULT_SENSITIVITY
+    print(f'       Sensor: {sensor} (sensitivity = {sensitivity} V/Pa)')
+
+    # Apply breakout box correction factor if provided
+    if input_args.breakout_box_factor:
+        sensitivity = sensitivity / input_args.breakout_box_factor
+        print('       Dividing sensitivity by breakout box factor of '
+              f'{input_args.breakout_box_factor}')
 
 print('------------------------------------------------------------------')
 print(f'Running cube2mseed on {len(raw_files)} raw file(s)...')
 print('------------------------------------------------------------------')
 
 for raw_file in raw_files:
     print(os.path.basename(raw_file))
+    # Default encoding is STEIM-1 (compressed integers) but we use uncompressed here
+    # since there is a direct NumPy data type equivalent (np.int32)
     args = ['cube2mseed', '--fringe-samples=NOMINAL', '--resample=SINC', f'--output-dir={tmp_dir}',
-            '--encoding=FLOAT-64', raw_file]
+            '--encoding=INT-32', raw_file]
     if input_args.verbose:
         args.append('--verbose')
     subprocess.call(args)
@@ -263,9 +274,18 @@
 
     tr.stats.channel = channel_id
 
-    tr.data = tr.data * BITWEIGHT    # Convert from counts to V
-    tr.data = tr.data + offset       # Remove voltage offset
-    tr.data = tr.data / sensitivity  # Convert from V to Pa
+    # Confidence check before we (possibly) modify the data type
+    assert tr.data.dtype == np.int32
+
+    # KEY: Modifying the time series of integer counts here!
+    if input_args.earthscope:  # Keep in integer counts!
+        output_encoding = 'INT32'
+    else:  # Convert all the way to Pa (float values)
+        tr.data = tr.data * BITWEIGHT    # Convert from counts to V
+        tr.data = tr.data + offset       # Remove voltage offset
+        tr.data = tr.data / sensitivity  # Convert from V to Pa
+        assert tr.data.dtype == np.float64
+        output_encoding = 'FLOAT64'
 
     if input_args.location == 'AUTO':
         if file.endswith('.pri0'):    # Channel 1
@@ -288,7 +308,7 @@
     t_min = np.min([t_min, tr.stats.starttime])
     t_max = np.max([t_max, tr.stats.endtime])
 
-    st.write(file, format='MSEED')
+    st.write(file, format='MSEED', encoding=output_encoding)
 
     # Define template for miniSEED renaming
     name_template = (f'{input_args.network}.{input_args.station}'