Add wave-statistics support for Vector instruments (#186)

* Add vector wave-statistics support

* More wave stats for vec

* add Vector test data

* add more Vec data

* Add Vector tests

* First try at implementing PUV

* Allow for pressure_sensor_height if present

* Doc updates

* Fix docs

* Change cutoff. This should be an option

* Add standard name and make it not be in an if-block.

* Enable setting of frequency cutoffs for PUV

* cleanup

* Tiny

* Doc update

doc/config.rst

@@ -192,9 +192,14 @@ Onset Hobo
 - ``names``: option for user specified column names (only recommended when code will not read names using automated/default method)
 
 Lowell TCM Hobo
-----------
+---------------
 
 - All the _min, _max, _bad_ens, etc. options available to the EXO.
 - ``skipfooter``: number of lines to skip in the CSV file at the end of the file
 - ``ncols``: number of columns of data to read, starting at first
 - ``names``: option for user specified column names (only recommended when code will not read names using automated/default method)
+
+Vector
+------
+- ``pressure_sensor_height`` and ``velocity_sample_volume_height`` to specify the elevations of these two sensors.
+- ``puv``: set to ``true`` to compute PUV wave statistics. **(EXPERIMENTAL)**

doc/vec.rst

@@ -29,3 +29,17 @@ runveccdf2nc.py
 .. argparse::
    :ref: stglib.core.cmd.veccdf2nc_parser
    :prog: runveccdf2nc.py
+
+Compute wave statistics
+=======================
+
+Use stglib's built-in wave-statistics toolbox to compute a wave-statistics file.
+
+**Experimental** PUV support is also present. Testing welcome!
+
+runvecnc2waves.py
+-----------------
+
+.. argparse::
+   :ref: stglib.core.cmd.vecnc2waves_parser
+   :prog: runvecnc2waves.py

doc/wvs.rst

@@ -31,12 +31,12 @@ runwvscdf2nc.py
   :ref: stglib.core.cmd.wvscdf2nc_parser
   :prog: runwvscdf2nc.py
 
-Compute waves statistics
-========================
+Compute wave statistics
+=======================
 
 Using DIWASP: This is run separately in Matlab using the DIWASP toolbox (docs TODO).
 
-Use stglib's built-in waves statistics toolbox to compute a wave statistics file.
+Use stglib's built-in wave-statistics toolbox to compute a wave-statistics file.
 
 runwvsnc2waves.py
 -----------------

setup.py

@@ -70,6 +70,7 @@
             "runeofecdf2nc.py=stglib.core.runcmd:runeofecdf2nc",
             "runvecdat2cdf.py=stglib.core.runcmd:runvecdat2cdf",
             "runveccdf2nc.py=stglib.core.runcmd:runveccdf2nc",
+            "runvecnc2waves.py=stglib.core.runcmd:runvecnc2waves",
             "runsigmat2cdf.py=stglib.core.runcmd:runsigmat2cdf",
             "runsigraw2cdf.py=stglib.core.runcmd:runsigraw2cdf",
             "runsigcdf2nc.py=stglib.core.runcmd:runsigcdf2nc",

stglib/aqd/aqdutils.py

@@ -1091,30 +1091,33 @@ def add_attributes(var, dsattrs):
         if "w_1204" in ds:
             ds["w_1204"].attrs["units"] = "m s-1"
 
-        if "vel1_1277" in ds:
-            ds["vel1_1277"].attrs.update(
-                {
-                    "units": "m s-1",
-                    "long_name": "Beam 1 Velocity",
-                    "epic_code": 1277,
-                }
-            )
-        if "vel2_1278" in ds:
-            ds["vel2_1278"].attrs.update(
-                {
-                    "units": "m s-1",
-                    "long_name": "Beam 2 Velocity",
-                    "epic_code": 1278,
-                }
-            )
-        if "vel3_1279" in ds:
-            ds["vel3_1279"].attrs.update(
-                {
-                    "units": "m s-1",
-                    "long_name": "Beam 3 Velocity",
-                    "epic_code": 1279,
-                }
-            )
+    if "vel1_1277" in ds:
+        ds["vel1_1277"].attrs.update(
+            {
+                "units": "m s-1",
+                "long_name": "Beam 1 Velocity",
+                "epic_code": 1277,
+                "standard_name": "radial_sea_water_velocity_away_from_instrument",
+            }
+        )
+    if "vel2_1278" in ds:
+        ds["vel2_1278"].attrs.update(
+            {
+                "units": "m s-1",
+                "long_name": "Beam 2 Velocity",
+                "epic_code": 1278,
+                "standard_name": "radial_sea_water_velocity_away_from_instrument",
+            }
+        )
+    if "vel3_1279" in ds:
+        ds["vel3_1279"].attrs.update(
+            {
+                "units": "m s-1",
+                "long_name": "Beam 3 Velocity",
+                "epic_code": 1279,
+                "standard_name": "radial_sea_water_velocity_away_from_instrument",
+            }
+        )
 
     if "AGC1_1221" in ds:
         ds["AGC1_1221"].attrs.update(

stglib/core/cmd.py

@@ -325,6 +325,14 @@ def veccdf2nc_parser():
     return parser
 
 
+def vecnc2waves_parser():
+    description = "Generate Vector waves statistics file"
+    parser = argparse.ArgumentParser(description=description)
+    ncarg(parser)
+
+    return parser
+
+
 def sigmat2cdf_parser():
     description = "Convert Signature files exported in Matlab format to raw .cdf format. Run this script from the directory containing Signature files."
     parser = argparse.ArgumentParser(description=description)

stglib/core/runcmd.py

@@ -256,6 +256,12 @@ def runvecdat2cdf():
     stglib.vec.dat2cdf.dat_to_cdf(metadata)
 
 
+def runvecnc2waves():
+    args = stglib.cmd.vecnc2waves_parser().parse_args()
+
+    stglib.vec.nc2waves.nc_to_waves(args.ncname)
+
+
 def runwvscdf2nc():
     args = stglib.cmd.wvscdf2nc_parser().parse_args()
 

stglib/core/waves.py

@@ -5,16 +5,18 @@
 
 
 def make_waves_ds(ds, noise=0.75):
-
-    print("Computing waves statistics")
+    print("Computing wave statistics")
     if "P_1ac" in ds:
         presvar = "P_1ac"
     else:
         presvar = "P_1"
 
     f, Pxx = pressure_spectra(ds[presvar].squeeze(), fs=1 / ds.attrs["sample_interval"])
 
-    z = ds.attrs["initial_instrument_height"]
+    if "pressure_sensor_height" in ds.attrs:
+        z = ds.attrs["pressure_sensor_height"]
+    else:
+        z = ds.attrs["initial_instrument_height"]
     h = ds[presvar].squeeze().mean(dim="sample") + z
 
     k = np.asarray([qkfs(2 * np.pi * f, x) for x in h.values])

stglib/tests/data/.gitattributes

@@ -86,3 +86,10 @@ L0221705_v2_rs.csv filter=lfs diff=lfs merge=lfs -text
 Bang_--_MBW.csv filter=lfs diff=lfs merge=lfs -text
 Barometer_w_RainGauge.csv filter=lfs diff=lfs merge=lfs -text
 1711426_Bang--MBW_0_CR.txt filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.hdr filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.pck filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.sen filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.ssl filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.vhd filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.VEC filter=lfs diff=lfs merge=lfs -text
+NBMCSB02.dat filter=lfs diff=lfs merge=lfs -text

stglib/tests/data/config_NBM22CSB.yaml

@@ -0,0 +1,7 @@
+Conventions: CF-1.8
+basefile: NBMCSB02
+filename: NBMCSBaq01
+initial_instrument_height: 0.5 # DJN guessed value
+orientation: UP # for Vector instruments, UP is when the head is pointing down
+AnalogInput1: "true"
+AnalogInput1_warmup_samples: 1

stglib/tests/data/gatts_NBM22CSB.txt

@@ -0,0 +1,33 @@
+MOORING;	NBM22CSB
+PROGRAM;	Sediment Transport in Coastal Environments
+PROJECT;	Temporal variability in sediment delivery to a North and a Central San Francisco Bay salt marsh 
+EXPERIMENT;	NBM22
+Region;	Northern and Central San Francisco Bay 
+CruiseID;	2022-619-FA
+SciPi;	Jessica Lacy, Samantha McGill
+Year;	2022
+Site;	CSB
+platform_type;	Tripod
+latitude;	37.92865
+longitude;  -122.48343
+WATER_DEPTH;	2.2
+WATER_DEPTH_source; Ships echosounder	
+WATER_DEPTH_note; Water depth at time of deployment
+NAVD88_elevation_ref;	999
+NAVD88_elevation_source;  RTK survey
+elevation_datum;	NAVD88
+Deployment_date;	14-JUN-2022 19:13
+Recovery_date;	10-AUG-2022 18:19
+OFAFunding;	 San Francisco Bay Regional Monitoring Program, San Francisco Estuary Institute
+CollaboratingAgency;	USGS WERC
+DESCRIPTION;	Coffee Table
+magnetic_variation;	13.24
+DATA_TYPE;	TIME
+DATA_SUBTYPE;	MOORED
+DATA_ORIGIN;	USGS/PCMSC
+COORD_SYSTEM;	GEOGRAPHIC
+Conventions;	PMEL/EPIC
+POS_CONST;	0
+DEPTH_CONST;	0
+history;	
+DATA_CMNT; none

stglib/tests/test_scripts.py

@@ -36,6 +36,18 @@ def exo_nc(nc_file, atmpres=None):
     assert "Done writing netCDF file" in result.stdout.decode("utf8")
 
 
+def test_exo():
+    exo_raw("glob_attbel5C.txt", "config_bel5C.yaml")
+    exo_nc("bel53Cexo-raw.cdf")
+    exo_raw("glob_att1119a.txt", "1119Aexo_config.yaml")
+    exo_nc("1119Aexo-raw.cdf")
+    exo_raw("glob_att1151b.txt", "1151Bexo_config.yaml")
+    exo_nc("1151Bexo-raw.cdf")
+    # test for atmospheric correction
+    exo_raw("glob_attbel3C.txt", "config_bel3C.yaml")
+    exo_nc("BEL19B3C03exo-raw.cdf", "atmpres_BEL19B3C03exo.cdf")
+
+
 def aqd_raw(glob_att, config_yaml):
     result = subprocess.run(
         [scripts / "runaqdhdr2cdf.py", glob_att, config_yaml],
@@ -76,18 +88,6 @@ def aqdhr_nc(nc_file):
     assert "Done writing netCDF file" in result.stdout.decode("utf8")
 
 
-def test_exo():
-    exo_raw("glob_attbel5C.txt", "config_bel5C.yaml")
-    exo_nc("bel53Cexo-raw.cdf")
-    exo_raw("glob_att1119a.txt", "1119Aexo_config.yaml")
-    exo_nc("1119Aexo-raw.cdf")
-    exo_raw("glob_att1151b.txt", "1151Bexo_config.yaml")
-    exo_nc("1151Bexo-raw.cdf")
-    # test for atmospheric correction
-    exo_raw("glob_attbel3C.txt", "config_bel3C.yaml")
-    exo_nc("BEL19B3C03exo-raw.cdf", "atmpres_BEL19B3C03exo.cdf")
-
-
 def test_aqd():
     aqd_raw("glob_att1118a_b.txt", "aqd1118A_config.yaml")
     aqd_nc("1118ABaqd-raw.cdf")
@@ -123,6 +123,43 @@ def test_aqdturnaround():
     aqdturnaround("BEL503")
 
 
+def vec_raw(glob_att, config_yaml):
+    result = subprocess.run(
+        [scripts / "runvecdat2cdf.py", glob_att, config_yaml],
+        capture_output=True,
+        cwd=cwd,
+    )
+    assert "Finished writing data" in result.stdout.decode("utf8")
+
+
+def vec_nc(nc_file, atmpres=None):
+    if atmpres is not None:
+        runlist = [scripts / "runveccdf2nc.py", nc_file, "--atmpres", atmpres]
+    else:
+        runlist = [scripts / "runveccdf2nc.py", nc_file]
+    result = subprocess.run(
+        runlist,
+        capture_output=True,
+        cwd=cwd,
+    )
+    assert "Done writing netCDF file" in result.stdout.decode("utf8")
+
+
+def vec_wvs(nc_file):
+    result = subprocess.run(
+        [scripts / "runvecnc2waves.py", nc_file],
+        capture_output=True,
+        cwd=cwd,
+    )
+    assert "Done writing netCDF file" in result.stdout.decode("utf8")
+
+
+def test_vec():
+    vec_raw("gatts_NBM22CSB.txt", "config_NBM22CSB.yaml")
+    vec_nc("NBMCSBaq01-raw.cdf")
+    vec_wvs("NBMCSBaq01-a.nc")
+
+
 def wxt_raw(glob_att, config_yaml):
     result = subprocess.run(
         [scripts / "runwxtcsv2cdf.py", glob_att, config_yaml],

stglib/vec/__init__.py

@@ -1 +1 @@
-from . import dat2cdf, cdf2nc
+from . import cdf2nc, dat2cdf, nc2waves

stglib/vec/dat2cdf.py

@@ -40,6 +40,7 @@ def dat_to_cdf(metadata):
 
     r = np.shape(dssen.Heading)[0]
     senburstlen = int(ds.attrs["VECSamplesPerBurst"] / ds.attrs["VECSamplingRate"] + 1)
+    ds.attrs["sample_interval"] = 1 / ds.attrs["VECSamplingRate"]
 
     mod = r % senburstlen
     if mod: