Add Seagauge waves (#265)

* Resolved rebasing conflicts

* Added __init__.py to sg folder

* Fixed setup.py conflict

* Add read_wb and wb_to_cdf

* Fixed typo in script name

* Add python script to sg __init__.py

* Fixed more conflicts

* Fixed more conflicts

* Made wave bursts

* Fixed nanoprecision error

* Added tide burst processing

* edited gatts for waves files

* Fixed wave burst timestamp

* edited metadata dropped

* Added tests

* Fixed tide wb interval issue

* Added documentation for sg

* Added error if no atmpres.cdf

* Added atmpres.cdf to data folder

* Added water level

* Added water level utility

* Force added atmpres.cdf

* Changed to new runots format

* Fixed issues

* Fixed drop sample issue

* streamlined run cmd

* Responded to Dan's comments

---------

Co-authored-by: De Meo <[email protected]>

doc/config.rst

@@ -226,4 +226,17 @@ SBE 37 MicroCAT
 
 SBE 26plus Seagauge
 -------------------
+For Seagauge tides:
+- ``file_type``: set to ``.tid`` or ``.wb`` to specify raw data file type
+- ``calculated_tide_interval``: enter the desired tide interval when using the .wb file
+- ``calculated_tide_interval_units``: tide interval units
+- ``calculated_tide_duration``: enter the desired tide duration when using the .wb file
+- ``calculated_tide_duration_units``: tide duration units
 - All the _min, _max, _bad_ens, etc. options available to the EXO.
+
+For Seagauge waves:
+- ``calculated_wave_interval``: enter the desired wave interval
+- ``calculated_wave_interval_units``: wave interval units
+- ``wp_min``, ``wp_max``: min/max allowable wave period, in seconds
+- ``wh_min``, ``wh_max``: min/max allowable wave height, in meters
+- ``wp_ratio``: maximum allowable ratio between peak period (``wp_peak``) and mean period (``wp_4060``)

doc/sgtid.rst

@@ -0,0 +1,4 @@
+SBE 26plus Seagauge (tides)
+***************************
+
+Starting from exported .tid or .wb file, use :doc:`runots.py </runots>` to process using the two :doc:`configuration files </config>`.

doc/sgwvs.rst

@@ -0,0 +1,4 @@
+SBE 26plus Seagauge (waves)
+***************************
+
+Starting from exported .wb file, use :doc:`runots.py </runots>` to process using the two :doc:`configuration files </config>`.

setup.py

@@ -80,9 +80,6 @@
             "runtcmcdf2nc.py=stglib.core.runcmd:runtcmcdf2nc",
             "runmcasc2cdf.py=stglib.core.runcmd:runmcasc2cdf",
             "runmccdf2nc.py=stglib.core.runcmd:runmccdf2nc",
-            "runsgtid2cdf.py=stglib.core.runcmd:runsgtid2cdf",
-            "runsgcdf2nc.py=stglib.core.runcmd:runsgcdf2nc",
-            "runsgnc2waves.py=stglib.core.runcmd:runsgnc2waves",
             "runots.py=stglib.core.runcmd:runots",
         ],
     },

stglib/__init__.py

@@ -23,5 +23,6 @@
 from .aqd import aqdutils
 from .core import cmd, filter, qaqc, utils, waves
 from .core.utils import read_globalatts
+from .sg import sgutils
 
 __version__ = _version.get_versions()["version"]

stglib/core/cmd.py

@@ -120,9 +120,13 @@ def runots_parser():
     addinst2cdf(instsp, "asc2cdf")
     addcdf2nc(instsp)
 
-    instsp = add_instrument(subparsers, "sg", "Seabird Seagauge")
+    instsp = add_instrument(subparsers, "sgtid", "Seabird Seagauge Tides")
     addinst2cdf(instsp, "tid2cdf")
     addcdf2nc(instsp)
+
+    instsp = add_instrument(subparsers, "sgwvs", "Seabird Seagauge Waves")
+    addinst2cdf(instsp, "wb2cdf")
+    addcdf2nc(instsp)
     addnc2waves(instsp)
 
     instsp = add_instrument(subparsers, "tcm", "Lowell Tilt Current Meter")
@@ -543,29 +547,3 @@ def mccdf2nc_parser():
     cdfarg(parser)
 
     return parser
-
-
-def sgtid2cdf_parser():
-    description = "Convert SBE 26plus Seagauge .tid files to raw .cdf format. Run this script from the directory containing Seagauge files."
-    parser = argparse.ArgumentParser(description=description)
-    gattsarg(parser)
-    yamlarg(parser)
-
-    return parser
-
-
-def sgcdf2nc_parser():
-    description = "Convert raw SBE 26plus Seagauge .cdf format to processed .nc files, optionally compensating for atmospheric pressure"
-    parser = argparse.ArgumentParser(description=description)
-    cdfarg(parser)
-    atmarg(parser)
-
-    return parser
-
-
-def sgnc2waves_parser():
-    description = "Generate SBE 26plus Seagauge waves statistics file"
-    parser = argparse.ArgumentParser(description=description)
-    ncarg(parser)
-
-    return parser

stglib/core/runcmd.py

@@ -349,32 +349,12 @@ def runmcasc2cdf(args=None):
     stglib.mc.asc_to_cdf(metadata)
 
 
-def runsgcdf2nc(args=None):
-    if not args:
-        args = stglib.cmd.sgcdf2nc_parser().parse_args()
-
-    run_cdf_to_nc(stglib.sg.cdf_to_nc, args)
-
-
-def runsgtid2cdf(args=None):
-    if not args:
-        args = stglib.cmd.sgtid2cdf_parser().parse_args()
-
-    metadata = get_metadata(args)
-
-    stglib.sg.tid_to_cdf(metadata)
-
-
-def runsgnc2waves(args=None):
-    if not args:
-        args = stglib.cmd.sgnc2waves_parser().parse_args()
-
-    stglib.sg.nc_to_waves(args.ncname)
-
-
 def runots():
     args = stglib.cmd.runots_parser().parse_args()
 
+    if "2cdf" in args.step:
+        metadata = get_metadata(args)
+
     if args.instrument == "aqd":
         if args.step == "hdr2cdf":
             runaqdhdr2cdf(args)
@@ -446,13 +426,18 @@ def runots():
             runmcasc2cdf(args)
         elif args.step == "cdf2nc":
             runmccdf2nc(args)
-    elif args.instrument == "sg":
+    elif args.instrument == "sgtid":
         if args.step == "tid2cdf":
-            runsgtid2cdf(args)
+            stglib.sg.tid2cdf.tid_to_cdf(metadata)
+        elif args.step == "cdf2nc":
+            run_cdf_to_nc(stglib.sg.cdf2nc.cdf_to_nc, args)
+    elif args.instrument == "sgwvs":
+        if args.step == "wb2cdf":
+            stglib.sg.wvswb2cdf.wb_to_cdf(metadata)
         elif args.step == "cdf2nc":
-            runsgcdf2nc(args)
+            run_cdf_to_nc(stglib.sg.wvscdf2nc.cdf_to_nc, args)
         elif args.step == "nc2waves":
-            runsgnc2waves(args)
+            stglib.sg.wvsnc2waves.nc_to_waves(args.ncname)
     elif args.instrument == "tcm":
         if args.step == "csv2cdf":
             runtcmcsv2cdf(args)

stglib/core/utils.py

@@ -1398,3 +1398,30 @@ def _todict(matobj):
         else:
             dic[strg] = elem
     return dic
+
+
+def create_water_level_var(ds):
+    """
+    Create water level variable from NAVD88 sensor height
+    """
+
+    if ds.z.attrs["geopotential_datum_name"] == "NAVD88":
+
+        if "sample" in ds.dims:
+            ds["water_level"] = xr.DataArray(
+                ds["P_1ac"].squeeze().mean(dim="sample") + ds["z"].values
+            )
+        else:
+            ds["water_level"] = ds["P_1ac"] + ds["z"].values
+
+        ds["water_level"].attrs["long_name"] = "Water level NAVD88"
+        ds["water_level"].attrs["units"] = "m"
+        ds["water_level"].attrs[
+            "standard_name"
+        ] = "sea_surface_height_above_geopotential_datum"
+        ds["water_level"].attrs["geopotential_datum_name"] = "NAVD88"
+    else:
+        print(
+            "Cannot create water_level variable without height_above_geopotential_datum relative to NAVD88 in global attributes file."
+        )
+    return ds

stglib/sg/__init__.py

@@ -0,0 +1 @@
+from . import cdf2nc, sgutils, tid2cdf, wvscdf2nc, wvsnc2waves, wvswb2cdf

stglib/sg/cdf2nc.py

@@ -0,0 +1,208 @@
+import math
+
+import numpy as np
+import pandas as pd
+import xarray as xr
+
+from ..core import utils
+from . import sgutils
+
+
+def cdf_to_nc(cdf_filename, atmpres=None):
+    """
+    Load a raw .cdf file and generate a processed .nc file
+    """
+
+    # Check for atmpres correction file
+    # Atmpres file is required for seagauge because pressure is measured as absolute pressure
+    if atmpres is None:
+        raise FileNotFoundError(
+            "The atmpres file does not exist. Atmpres file is required for Seagauge because pressure is measured as absolute pressure."
+        )
+
+    # Load raw .cdf data
+    ds = xr.open_dataset(cdf_filename)
+
+    # Atmospheric pressure correction
+    if ds.attrs["file_type"] == ".tid":
+        ds = utils.atmos_correct(ds, atmpres)
+    elif ds.attrs["file_type"] == ".wb":
+        ds = sgutils.atmos_correct_burst(ds, atmpres)
+
+    # If .wb file, split into tide bursts
+    if ds.attrs["file_type"] == ".wb":
+        ds = avg_tide_bursts(ds)
+
+    # Rename variables to CF compliant names
+    ds = ds_rename_vars(ds)
+
+    # Add attributes
+    ds = sgutils.ds_add_attrs(ds)
+
+    # Edit metadata depending on file type
+    if ds.attrs["file_type"] == ".tid":
+        ds = ds_drop_tid(ds)
+        ds = ds.drop("sample")  # Drop sample variable
+    elif ds.attrs["file_type"] == ".wb":
+        ds = ds_drop_wb(ds)
+
+    # Call QAQC
+    ds = sgutils.sg_qaqc(ds)
+
+    # Run utilities
+    ds = utils.clip_ds(ds)
+    ds = utils.create_nominal_instrument_depth(ds)
+    ds = utils.create_z(ds)
+    ds = utils.create_water_level_var(ds)
+    ds = utils.create_water_depth_var(ds)
+    ds = utils.ds_add_lat_lon(ds)
+    ds = utils.add_start_stop_time(ds)
+    ds = utils.add_delta_t(ds)
+    ds = utils.add_min_max(ds)
+
+    # Write to .nc file
+    print("Writing cleaned/trimmed data to .nc file")
+    nc_filename = ds.attrs["filename"] + "s-tide-a.nc"
+
+    ds.to_netcdf(
+        nc_filename, unlimited_dims=["time"], encoding={"time": {"dtype": "i4"}}
+    )
+    utils.check_compliance(nc_filename, conventions=ds.attrs["Conventions"])
+
+    print(f"Done writing netCDF file {nc_filename}")
+
+
+def ds_rename_vars(ds):
+    """
+    Rename variables to be CF compliant
+    """
+    varnames = {"Temp": "T_28"}
+
+    # Check to make sure they exist before trying to rename
+    newvars = {}
+    for k in varnames:
+        if k in ds:
+            newvars[k] = varnames[k]
+    return ds.rename(newvars)
+
+
+def ds_drop_tid(ds):
+    """
+    Drop global attribute metadata not needed for .tid file
+    """
+    gatts = [
+        "WaveInterval",
+        "WaveIntervalUnits",
+        "WaveSamples",
+        "sample_rate",
+        "sample_rate_units",
+        "BurstDuration",
+        "BurstDurationUnits",
+        "WaveBurstsPerDay",
+        "NumberOfWaveBursts",
+    ]
+
+    # Check to make sure they exist
+    for k in gatts:
+        if k in ds.attrs:
+            del ds.attrs[k]
+    return ds
+
+
+def ds_drop_wb(ds):
+    """
+    Drop global attribute metadata not needed for .wb file
+    """
+    gatts = [
+        "WaveInterval",
+        "WaveIntervalUnits",
+        "WaveSamples",
+        "sample_rate",
+        "sample_rate_units",
+        "BurstDuration",
+        "BurstDurationUnits",
+        "WaveBurstsPerDay",
+        "NumberOfWaveBursts",
+        "TideInterval",
+        "TideIntervalUnits",
+        "TideDuration",
+        "TideDurationUnits",
+        "TideSamplesPerDay",
+        "NumberOfTideMeasurements",
+    ]
+
+    # Check to make sure they exist
+    for k in gatts:
+        if k in ds.attrs:
+            del ds.attrs[k]
+    return ds
+
+
+def avg_tide_bursts(ds):
+
+    # Calculate how many rows to subdivide each wave burst (round up)
+    rows = math.ceil(
+        float(ds.attrs["BurstDuration"]) / float(ds.attrs["calculated_tide_interval"])
+    )
+
+    # Calculate how many columns to subdivide each wave burst
+    cols = int(
+        float(ds.attrs["calculated_tide_interval"]) * float(ds.attrs["sample_rate"])
+    )
+
+    # Calculate number of values to average based on tide duration
+    values_avg = int(
+        float(ds.attrs["calculated_tide_duration"]) * float(ds.attrs["sample_rate"])
+    )
+
+    # Define time interval
+    delta_t = int(ds.attrs["calculated_tide_interval"])
+    delta_t = f"{delta_t}s"
+
+    # Reshape P_1 and P_1ac pressure bursts and average
+    P1_avg = []
+    P1ac_avg = []
+    timestamp = []
+    for t in range(0, ds["time"].size):
+
+        # Make timestamp
+        new_time = np.array(
+            pd.date_range(start=ds.time[t].values, periods=rows, freq=delta_t)
+        )
+        timestamp.append(new_time)
+
+        # Reshape adding nans to end of array
+        no_pads = rows * cols - int(ds.attrs["WaveSamples"])
+        new_P1_burst = np.pad(
+            ds.P_1[t].values, (0, no_pads), mode="constant", constant_values=np.nan
+        ).reshape(rows, cols)
+        new_P1ac_burst = np.pad(
+            ds.P_1ac[t].values, (0, no_pads), mode="constant", constant_values=np.nan
+        ).reshape(rows, cols)
+
+        # Average burst for tide duration
+        for j in range(0, new_P1_burst.shape[0]):
+            avg_pres_P1 = np.mean(new_P1_burst[j, slice(0, values_avg)])
+            avg_pres_P1ac = np.mean(new_P1ac_burst[j, slice(0, values_avg)])
+            P1_avg.append(avg_pres_P1)
+            P1ac_avg.append(avg_pres_P1ac)
+
+    # Append all timestamps to array
+    date = np.concatenate(timestamp)
+
+    # Copy global attributes
+    attr_copy = ds.attrs.copy()
+
+    # Delete old xarray because no longer need burst data
+    del ds
+
+    # Make new xarray with averaged data
+    ds = xr.Dataset(
+        data_vars=dict(P_1=(["time"], P1_avg), P_1ac=(["time"], P1ac_avg)),
+        coords=dict(time=date),
+    )
+
+    # Add back global attributes
+    ds.attrs = attr_copy
+
+    return ds