feat(waterdata.xarray): record vertical_datum to distinguish stage parameters

Gage height (00065) and stream water level (63160) share the CF standard_name
water_surface_height_above_reference_datum, which left them indistinguishable.
Attach a vertical_datum attribute (matching USGS's column name) recording the
reference datum: "local site datum" for gage height (00065), "NAVD88" for
stream water level (63160).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: thodson-usgs

dataretrieval/waterdata/types.py

@@ -123,3 +123,13 @@
     "63160": "water_surface_height_above_reference_datum",
     "00045": "lwe_thickness_of_precipitation_amount",
 }
+
+# USGS parameter code -> vertical reference datum, attached as a
+# ``vertical_datum`` attribute. The two water-surface-height parameters share
+# the CF standard_name water_surface_height_above_reference_datum, so the datum
+# distinguishes them: gage height (00065) is measured from a local site (gage)
+# datum, while stream water level (63160) is referenced to NAVD88.
+CF_VERTICAL_DATUM = {
+    "00065": "local site datum",
+    "63160": "NAVD88",
+}

dataretrieval/waterdata/xarray.py

@@ -59,7 +59,12 @@
 
 from . import api as _api
 from .nearest import get_nearest_continuous as _get_nearest_continuous
-from .types import CF_CELL_METHODS, CF_STANDARD_NAMES, CF_UNIT_MAP
+from .types import (
+    CF_CELL_METHODS,
+    CF_STANDARD_NAMES,
+    CF_UNIT_MAP,
+    CF_VERTICAL_DATUM,
+)
 
 __all__ = [
     "get_continuous",
@@ -250,6 +255,9 @@ def _var_attrs(desc, *, unit, pcode, stat, default_cell_method, ancillary, name)
         sn = CF_STANDARD_NAMES.get(str(pcode))
         if sn:
             attrs["standard_name"] = sn
+        datum = CF_VERTICAL_DATUM.get(str(pcode))
+        if datum:
+            attrs["vertical_datum"] = datum
         attrs["usgs_parameter_code"] = str(pcode)
 
     if stat is not None and _pd.notna(stat):

tests/waterdata_xarray_test.py

@@ -100,6 +100,24 @@ def test_missing_standard_name_is_omitted():
     assert ds["mystery"].attrs["usgs_parameter_code"] == "99999"
 
 
+def test_vertical_datum_distinguishes_stage_parameters():
+    # 00065 (gage height) and 63160 (water level above NAVD88) share the CF
+    # standard_name water_surface_height_above_reference_datum; the vertical_datum
+    # attribute records the differing reference datum so they're distinguishable.
+    for pcode, datum in (("00065", "local site datum"), ("63160", "NAVD88")):
+        df = _daily_frame()
+        df["parameter_code"] = pcode
+        ds = wdx._build_ragged(df, _meta(), service="continuous", series_meta={})
+        v = ds["value"]
+        assert v.attrs["standard_name"] == "water_surface_height_above_reference_datum"
+        assert v.attrs["vertical_datum"] == datum
+    # a parameter with no datum mapping gets no vertical_datum attr
+    ds = wdx._build_ragged(
+        _daily_frame(), _meta(), service="daily", series_meta=_DISCHARGE_META
+    )
+    assert "vertical_datum" not in ds["value"].attrs
+
+
 def test_multiple_parameters_outer_join_on_time():
     # discharge at t1,t2 ; temperature at t2,t3 -> union time, NaN fill
     q = _daily_frame(values=(100, 110), times=("2024-06-01", "2024-06-02"))