Thermal lag effects

src/glide/assets/config.yml

@@ -48,8 +48,9 @@ netcdf_attributes:
 # variable in it from the L2 output.
 
 include:
-  flight: true   # fin, battpos, heading, ballast_pumped
-  thermo: true   # salinity, density, rho0, SA, CT, N2, depth, z, sound_speed
+  flight: true            # fin, battpos, heading, ballast_pumped
+  thermo: true            # salinity, density, rho0, SA, CT, N2, depth, z, sound_speed
+  ctd_corrections: true   # rbrctd_* and ctd41cp_* raw variables
 
 
 # =============================================================================
@@ -139,6 +140,19 @@ qc: {}
 l1_variables: {}
 
 
+# =============================================================================
+# CTD CORRECTIONS
+# =============================================================================
+# Sensor-specific lag corrections applied to the CTD temperature before the
+# science dataset is merged with flight. Requires
+# `include: ctd_corrections: true` so the sensor-native source variables are
+# pulled out of the raw science file. 
+
+ctd:
+  rbrctd:
+    temperature_lag: 0.9   # seconds; positive means sensor lags water
+
+
 # =============================================================================
 # MERGED VARIABLES
 # =============================================================================

src/glide/assets/core.yml

@@ -540,6 +540,44 @@ suites:
         units: s-2
         observation_type: calculated
 
+  # ---------------------------------------------------------------------------
+  ctd_corrections:  # CTD-native variables used by ctd.correct_ctd
+  # ---------------------------------------------------------------------------
+  # Used to correct for thermal mass and lag effects in temperature and salinity.
+
+    rbrctd_time:
+      source: sci_rbrctd_timestamp
+      drop_from_l2: True
+      CF:
+        long_name: RBR CTD native timestamp
+        standard_name: time
+        units: seconds since 1970-01-01T00:00:00Z
+        observation_type: measured
+
+    rbrctd_temperature:
+      source: sci_rbrctd_temperature_00
+      drop_from_l2: True
+      CF:
+        long_name: RBR CTD temperature
+        standard_name: sea_water_temperature
+        units: celsius
+        instrument: instrument_ctd
+        observation_type: measured
+        valid_min: -5.0
+        valid_max: 50.0
+
+    rbrctd_conductivity:
+      source: sci_rbrctd_conductivity_00
+      drop_from_l2: True
+      CF:
+        long_name: RBR CTD conductivity
+        standard_name: sea_water_electrical_conductivity
+        units: mS cm-1
+        instrument: instrument_ctd
+        observation_type: measured
+        valid_min: 0.1
+        valid_max: 10.0
+
   # ---------------------------------------------------------------------------
   flight_model:  # output variables from the steady-state flight model
   # ---------------------------------------------------------------------------

src/glide/cli.py

@@ -14,6 +14,7 @@
 from . import (
     ancillery,
     config,
+    ctd,
     gliderdac,
     hotel,
     process_l1,
@@ -227,6 +228,8 @@ def l2(
     flt = process_l1.apply_qc(flt, conf)
     sci = process_l1.apply_qc(sci, conf)
 
+    sci = ctd.correct_ctd(sci, conf)
+
     merged = process_l1.merge(flt, sci, conf, "science")
 
     merged = process_l1.calculate_thermodynamics(merged, conf)

src/glide/config.py

@@ -19,6 +19,7 @@
     "qc",
     "l1_variables",
     "merged_variables",
+    "ctd",
 )
 
 # Helper functions
@@ -202,6 +203,7 @@ def _section(name: str) -> dict:
     l1_vars = _section("l1_variables")
     merged_vars = _section("merged_variables")
     flight_model_cfg = _section("flight")
+    ctd_cfg = _section("ctd")
 
     # Build variable set: core + enabled optional suites
     variables = _merge_suites(core, suites, include)
@@ -244,6 +246,7 @@ def _section(name: str) -> dict:
         ngdac=ngdac,
         include=include,
         flight=flight_model_cfg,
+        ctd=ctd_cfg,
     )
 
     return config

src/glide/ctd.py

@@ -0,0 +1,126 @@
+# Applies sensor-specific lag corrections to the CTD temperature, then
+# re-interpolates temperature and conductivity from the CTD-native timestamp
+# grid back onto the science time grid.
+
+from __future__ import annotations
+
+import logging
+
+import numpy as np
+import xarray as xr
+
+_log = logging.getLogger(__name__)
+
+# Canonical names of the rbrctd suite variables (after format_l1 renaming).
+_RBRCTD_VARS = (
+    "rbrctd_time",
+    "rbrctd_temperature",
+    "rbrctd_conductivity",
+)
+
+_TS_SENTINEL = 946684800  # 2000-01-01T00:00:00Z
+
+
+def correct_ctd(sci: xr.Dataset, config: dict) -> xr.Dataset:
+    """Apply CTD lag correction to the science dataset.
+
+    For the RBR Concerto: shift temperature earlier in time on the native
+    ``rbrctd_time`` grid by ``ctd.rbrctd.temperature_lag`` seconds, then
+    interpolate the lag-shifted temperature and the unadjusted conductivity
+    back onto ``sci.time``, overwriting ``temperature`` and ``conductivity``.
+
+    Returns ``sci`` unchanged when the rbrctd variables are not present.
+    Other CTDs are not yet supported.
+    """
+    if not all(v in sci.variables for v in _RBRCTD_VARS):
+        _log.debug(
+            "rbrctd variables not present in science dataset; skipping CTD correction"
+        )
+        return sci
+
+    lag_s = float(
+        ((config.get("ctd") or {}).get("rbrctd") or {}).get("temperature_lag", 0.9)
+    )
+
+    t_native, T, C = _build_native_rbrctd(sci)
+    if t_native.size < 2:
+        _log.warning("rbrctd has fewer than 2 valid samples; skipping correction")
+        return sci
+
+    T_lag = _apply_lag(t_native, T, lag_s)
+
+    sci_t = _time_as_seconds(sci["time"])
+    T_on_sci = np.interp(sci_t, t_native, T_lag, left=np.nan, right=np.nan)
+    C_on_sci = np.interp(sci_t, t_native, C, left=np.nan, right=np.nan)
+
+    sci = _overwrite(
+        sci,
+        "temperature",
+        T_on_sci,
+        f"lag-corrected (lag={lag_s}s) via ctd.correct_ctd",
+    )
+    sci = _overwrite(
+        sci,
+        "conductivity",
+        C_on_sci,
+        "interpolated from rbrctd native grid via ctd.correct_ctd",
+    )
+    return sci
+
+
+def _build_native_rbrctd(
+    sci: xr.Dataset,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Return (time_s, T, C) on the native rbrctd timestamp grid.
+
+    Drops rows where the native timestamp is missing or below the sentinel,
+    sorts ascending, and removes duplicate timestamps.
+    """
+    t = np.asarray(sci["rbrctd_time"].values, dtype="f8")
+    T = np.asarray(sci["rbrctd_temperature"].values, dtype="f8")
+    C = np.asarray(sci["rbrctd_conductivity"].values, dtype="f8")
+
+    valid = np.isfinite(t) & (t > _TS_SENTINEL)
+    t, T, C = t[valid], T[valid], C[valid]
+
+    order = np.argsort(t, kind="stable")
+    t, T, C = t[order], T[order], C[order]
+
+    _, keep = np.unique(t, return_index=True)
+    keep.sort()
+    return t[keep], T[keep], C[keep]
+
+
+def _apply_lag(time_s: np.ndarray, T: np.ndarray, lag_s: float) -> np.ndarray:
+    """Shift ``T`` earlier in time by ``lag_s`` seconds.
+
+    The sensor reports later than the water it samples, so the "true"
+    temperature at time ``t`` equals the value sampled at ``t + lag``.
+    Equivalently, interpolate ``T`` from ``time - lag`` back onto ``time``.
+    """
+    if lag_s == 0.0:
+        return T.copy()
+    shifted = time_s - lag_s
+    return np.interp(time_s, shifted, T)
+
+
+def _overwrite(
+    sci: xr.Dataset, name: str, values: np.ndarray, comment: str
+) -> xr.Dataset:
+    """Overwrite ``sci[name]`` while preserving its CF attrs and appending a comment."""
+    if name not in sci.variables:
+        _log.warning("Cannot overwrite %s: not present in science dataset", name)
+        return sci
+    attrs = dict(sci[name].attrs)
+    existing = attrs.get("comment", "")
+    attrs["comment"] = f"{existing}; {comment}".lstrip("; ")
+    sci[name] = (sci[name].dims, values, attrs)
+    return sci
+
+
+def _time_as_seconds(t: xr.DataArray) -> np.ndarray:
+    """Return ``t`` as float64 posix seconds regardless of dtype."""
+    vals = t.values
+    if np.issubdtype(vals.dtype, np.datetime64):
+        return vals.astype("datetime64[ns]").astype("f8") / 1e9
+    return np.asarray(vals, dtype="f8")

tests/test_config.py

@@ -95,6 +95,7 @@ def test_load_config_exclude_thermo():
         qc: {}
         l1_variables: {}
         merged_variables: {}
+        ctd: {}
     """)
     try:
         conf = config.load_config(path)
@@ -123,6 +124,7 @@ def test_load_config_qc_override_applied():
             max_gap: 900
         l1_variables: {}
         merged_variables: {}
+        ctd: {}
     """)
     try:
         conf = config.load_config(path)
@@ -154,6 +156,7 @@ def test_load_config_companion_variable():
               long_name: GPS horizontal dilution of precision
               units: "1"
         merged_variables: {}
+        ctd: {}
     """)
     try:
         conf = config.load_config(path)