Add waterlevel-trends data product (one feature per well)

New ogc_waterlevel_trend product: nm_waterlevel_trends. One GeoJSON
feature per well with a depth-to-water trend computed from the well's
water-level timeseries, across all NM sources.

Algorithm ported from the Ocotillo API ogc_depth_to_water_trend_wells
materialized view:
- least-squares slope (REGR_SLOPE) of depth-to-water (ft) vs observation
  time, scaled to ft/year (365.25-day year);
- classify only with >=10 records, or >=4 records spanning >=2 years;
  slope > 0.25 ft/yr "increasing" (deeper/declining), < -0.25
  "decreasing", else "stable"; otherwise "not enough data".

The product includes a trend_method description of how the trend was
calculated (TREND_METHOD_DESCRIPTION).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: jirhiker

backend/persisters/ogc_features.py

@@ -11,11 +11,87 @@
 from datetime import datetime, timezone
 from typing import Optional
 
+# Seconds per Julian year (365.25 days) — matches Ocotillo's trend MV, which
+# divides the per-second regression slope by 31557600 to get ft/year.
+_SECONDS_PER_YEAR = 31557600.0
+
+# Trend classification thresholds, ported verbatim from the Ocotillo
+# ogc_depth_to_water_trend_wells materialized view.
+_TREND_SLOPE_THRESHOLD = 0.25  # ft/year
+_TREND_MIN_RECORDS = 10
+_TREND_MIN_RECORDS_WITH_SPAN = 4
+_TREND_MIN_SPAN_YEARS = 2.0
+
+# Human-readable description of the trend method, embedded in the product so
+# consumers know how the classification was derived.
+TREND_METHOD_DESCRIPTION = (
+    "Depth-to-water trend per well. A least-squares linear regression "
+    "(equivalent to PostgreSQL REGR_SLOPE) is fit to depth-to-water-below-"
+    "ground-surface (feet) against observation time; the per-second slope is "
+    "scaled by 31,557,600 s/yr (a 365.25-day year) to slope_ft_per_year. A well "
+    "is classified only when it has at least 10 measurements, or at least 4 "
+    "measurements spanning at least 2 years; otherwise 'not enough data'. When "
+    "classified: slope > 0.25 ft/yr is 'increasing' (water level getting "
+    "DEEPER, i.e. a declining water table), slope < -0.25 ft/yr is 'decreasing' "
+    "(water level getting SHALLOWER), otherwise 'stable'. Ported from the "
+    "Ocotillo API ogc_depth_to_water_trend_wells materialized view."
+)
+
 
 def _timestamp_now() -> str:
     return datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")
 
 
+def _parse_epoch_seconds(date, time) -> Optional[float]:
+    """Best-effort parse of a DIE date (+optional time) to POSIX seconds (UTC)."""
+    if not date:
+        return None
+    text = f"{date}T{time}" if time else str(date)
+    text = text.replace("Z", "")
+    for parse in (datetime.fromisoformat,):
+        try:
+            dt = parse(text)
+            if dt.tzinfo is None:
+                dt = dt.replace(tzinfo=timezone.utc)
+            return dt.timestamp()
+        except (ValueError, TypeError):
+            pass
+    try:
+        dt = datetime.fromisoformat(str(date)).replace(tzinfo=timezone.utc)
+        return dt.timestamp()
+    except (ValueError, TypeError):
+        return None
+
+
+def _regr_slope(xs: list, ys: list) -> Optional[float]:
+    """Least-squares slope of ys on xs (PostgreSQL REGR_SLOPE). None if x has no
+    variance (fewer than 2 distinct x values)."""
+    n = len(xs)
+    if n < 2:
+        return None
+    mx = sum(xs) / n
+    my = sum(ys) / n
+    sxx = sum((x - mx) ** 2 for x in xs)
+    if sxx == 0:
+        return None
+    sxy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
+    return sxy / sxx
+
+
+def _classify_trend(slope_ft_per_year, record_count, span_years) -> str:
+    qualifies = record_count >= _TREND_MIN_RECORDS or (
+        record_count >= _TREND_MIN_RECORDS_WITH_SPAN
+        and span_years >= _TREND_MIN_SPAN_YEARS
+    )
+    if not qualifies or slope_ft_per_year is None:
+        return "not enough data"
+    if slope_ft_per_year > _TREND_SLOPE_THRESHOLD:
+        return "increasing"
+    if slope_ft_per_year < -_TREND_SLOPE_THRESHOLD:
+        return "decreasing"
+    return "stable"
+
+
 def _make_feature(record, collection_id: str) -> dict:
     """Build one OGC-compliant Feature from a SummaryRecord or SiteRecord."""
     source = getattr(record, "source", "")
@@ -174,6 +250,125 @@ def dump_major_chemistry_collection(path: str, records: list, meta: dict) -> dic
     return collection
 
 
+def dump_waterlevel_trend_collection(
+    path: str,
+    site_records: list,
+    timeseries_records: list,
+    meta: dict,
+) -> dict:
+    """
+    Write an OGC FeatureCollection of per-well depth-to-water trends to *path*.
+    One Feature per well.
+
+    *site_records* and *timeseries_records* are index-aligned: ``site_records[i]``
+    is the well and ``timeseries_records[i]`` is its list of ParameterRecord
+    observations (DIE water-level values are already depth-to-water below ground
+    surface in feet, so no measuring-point adjustment is applied here).
+
+    Each feature carries: record_count, first/last_observation_datetime,
+    span_years, slope_ft_per_year, trend_category, well_depth(+units). The
+    collection carries ``trend_method`` describing the calculation. See
+    TREND_METHOD_DESCRIPTION.
+
+    §V: MUST include top-level id, type, numberReturned, timeStamp.
+    §V: Each Feature MUST have top-level id.
+    """
+    collection_id = meta.get("id", "collection")
+
+    features = []
+    for site, obs_list in zip(site_records, timeseries_records):
+        pairs = []
+        for obs in obs_list:
+            value = getattr(obs, "parameter_value", None)
+            epoch = _parse_epoch_seconds(
+                getattr(obs, "date_measured", None), getattr(obs, "time_measured", None)
+            )
+            if value is None or epoch is None:
+                continue
+            try:
+                pairs.append((epoch, float(value)))
+            except (TypeError, ValueError):
+                continue
+
+        pairs.sort(key=lambda p: p[0])
+        record_count = len(pairs)
+        xs = [p[0] for p in pairs]
+        ys = [p[1] for p in pairs]
+
+        if record_count >= 2:
+            span_years = (xs[-1] - xs[0]) / _SECONDS_PER_YEAR
+            slope = _regr_slope(xs, ys)
+            slope_ft_per_year = None if slope is None else slope * _SECONDS_PER_YEAR
+        else:
+            span_years = 0.0
+            slope_ft_per_year = None
+
+        trend_category = _classify_trend(slope_ft_per_year, record_count, span_years)
+
+        source = getattr(site, "source", "") or ""
+        rid = getattr(site, "id", "") or ""
+        feature_id = f"{source}:{rid}" if source and rid else str(rid)
+
+        def _iso(epoch):
+            if epoch is None:
+                return None
+            return datetime.fromtimestamp(epoch, tz=timezone.utc).strftime(
+                "%Y-%m-%dT%H:%M:%SZ"
+            )
+
+        props = {
+            "source": source,
+            "id": rid,
+            "name": getattr(site, "name", None),
+            "well_depth": getattr(site, "well_depth", None),
+            "well_depth_units": getattr(site, "well_depth_units", None),
+            "record_count": record_count,
+            "first_observation_datetime": _iso(xs[0]) if record_count else None,
+            "last_observation_datetime": _iso(xs[-1]) if record_count else None,
+            "span_years": round(span_years, 3),
+            "slope_ft_per_year": (
+                None if slope_ft_per_year is None else round(slope_ft_per_year, 4)
+            ),
+            "trend_category": trend_category,
+        }
+
+        lat = getattr(site, "latitude", None)
+        lon = getattr(site, "longitude", None)
+        elev = getattr(site, "elevation", None)
+        coords = [lon, lat] if elev is None else [lon, lat, elev]
+
+        features.append({
+            "type": "Feature",
+            "id": feature_id,
+            "geometry": {"type": "Point", "coordinates": coords},
+            "properties": props,
+        })
+
+    collection = {
+        "type": "FeatureCollection",
+        "id": collection_id,
+        "title": meta.get("title", collection_id),
+        "description": meta.get("description", ""),
+        "trend_method": TREND_METHOD_DESCRIPTION,
+        "timeStamp": _timestamp_now(),
+        "numberMatched": len(features),
+        "numberReturned": len(features),
+        "links": [
+            {
+                "href": meta.get("href", ""),
+                "rel": "self",
+                "type": "application/geo+json",
+            }
+        ],
+        "features": features,
+    }
+
+    with open(path, "w", encoding="utf-8") as f:
+        json.dump(collection, f, indent=2, default=str)
+
+    return collection
+
+
 def dump_timeseries_collection(
     path: str,
     site_records: list,

orchestration/assets/products.py

@@ -37,6 +37,7 @@
     dump_major_chemistry_collection,
     dump_summary_collection,
     dump_timeseries_collection,
+    dump_waterlevel_trend_collection,
 )
 from backend.record import ParameterRecord, SiteRecord, SummaryRecord
 from backend.unifier import unify_source
@@ -167,8 +168,8 @@ def _build_combine_asset(product: dict, source_keys: list, source_asset_keys: li
 
     Depends on every source asset (wired via ``ins``), merges their
     records/sites/timeseries, writes the OGC GeoJSON collection — summary,
-    timeseries, or major-chemistry depending on ``output_type`` — and uploads it
-    to GCS."""
+    timeseries, major-chemistry, or waterlevel-trend depending on
+    ``output_type`` — and uploads it to GCS."""
     pid = product["id"]
     output_type = product["output_type"]
     ins = {
@@ -204,6 +205,14 @@ def _combine_asset(
             elif output_type == "ogc_summary":
                 records = [SummaryRecord(p) for p in all_records]
                 dump_summary_collection(str(out), records, meta)
+            elif output_type == "ogc_waterlevel_trend":
+                # site_records and the per-site observation lists are index
+                # aligned (see source asset); compute one trend per well.
+                site_records = [SiteRecord(p) for p in all_sites]
+                series = [
+                    [ParameterRecord(o) for o in site_ts] for site_ts in all_timeseries
+                ]
+                dump_waterlevel_trend_collection(str(out), site_records, series, meta)
             else:
                 site_records = [SiteRecord(p) for p in all_sites]
                 flat = [

orchestration/config/products.yaml

@@ -68,3 +68,18 @@ products:
       state: NM
     sources:
       exclude: []
+
+  # One feature per well with a depth-to-water trend (slope_ft_per_year +
+  # category) computed from each well's water-level timeseries. The collection
+  # carries a trend_method description. Algorithm ported from the Ocotillo API
+  # ogc_depth_to_water_trend_wells materialized view (see ogc_features.py).
+  - id: nm_waterlevel_trends
+    parameter: waterlevels
+    output_type: ogc_waterlevel_trend
+    title: "NM Water Level Trends"
+    description: "Per-well depth-to-water trend (slope and category), all NM sources"
+    schedule: "0 12 * * *"
+    spatial_filter:
+      state: NM
+    sources:
+      exclude: []

orchestration/definitions.py

@@ -11,7 +11,12 @@
 
 _PRODUCTS_PATH = Path(__file__).parent / "config" / "products.yaml"
 
-_SUPPORTED_OUTPUT_TYPES = {"ogc_summary", "ogc_timeseries", "ogc_major_chemistry"}
+_SUPPORTED_OUTPUT_TYPES = {
+    "ogc_summary",
+    "ogc_timeseries",
+    "ogc_major_chemistry",
+    "ogc_waterlevel_trend",
+}
 
 
 def _load_products() -> dict:

tests/test_persisters/test_ogc_features.py

@@ -228,3 +228,53 @@ def test_geometry_and_required_fields(self, tmp_path):
         assert "timeStamp" in result
         feat = result["features"][0]
         assert feat["geometry"]["coordinates"] == [-106.0, 34.0]
+
+
+from backend.persisters.ogc_features import dump_waterlevel_trend_collection
+
+
+def _trend_site(source="NMBGMR", rid="W1", well_depth=100.0):
+    return SiteRecord({
+        "source": source, "id": rid, "name": f"Well {rid}",
+        "latitude": 34.0, "longitude": -106.0, "elevation": None,
+        "well_depth": well_depth, "well_depth_units": "ft",
+    })
+
+
+def _trend_obs(date, value):
+    return ParameterRecord({"parameter_value": value, "date_measured": date, "time_measured": None})
+
+
+class TestWaterLevelTrendCollection:
+    def test_classifies_trends_and_carries_method(self, tmp_path):
+        increasing = [_trend_obs(f"{2010 + i}-01-01", 50.0 + 0.5 * i) for i in range(12)]
+        stable = [_trend_obs(f"{2010 + i}-01-01", 50.0) for i in range(12)]
+        sparse = [_trend_obs("2010-01-01", 50.0), _trend_obs("2011-01-01", 51.0), _trend_obs("2012-01-01", 52.0)]
+        decreasing = [_trend_obs(f"{2010 + i}-01-01", 60.0 - 1.0 * i) for i in range(5)]
+
+        sites = [_trend_site(rid="A"), _trend_site(rid="B"), _trend_site("NWIS", "C"), _trend_site("PVACD", "D")]
+        series = [increasing, stable, sparse, decreasing]
+
+        out = tmp_path / "tr.geojson"
+        result = dump_waterlevel_trend_collection(str(out), sites, series, {"id": "nm_waterlevel_trends"})
+
+        assert result["numberReturned"] == 4
+        assert "trend_method" in result and result["trend_method"]
+        by_id = {f["id"]: f["properties"] for f in result["features"]}
+
+        assert by_id["NMBGMR:A"]["trend_category"] == "increasing"
+        assert round(by_id["NMBGMR:A"]["slope_ft_per_year"], 2) == 0.5
+        assert by_id["NMBGMR:B"]["trend_category"] == "stable"
+        assert by_id["NWIS:C"]["trend_category"] == "not enough data"  # only 3 records
+        assert by_id["PVACD:D"]["trend_category"] == "decreasing"      # 5 records / 4 yr span
+
+    def test_required_fields_and_geometry(self, tmp_path):
+        sites = [_trend_site(rid="W1")]
+        series = [[_trend_obs("2010-01-01", 50.0)]]
+        out = tmp_path / "tr.geojson"
+        result = dump_waterlevel_trend_collection(str(out), sites, series, {"id": "nm_waterlevel_trends"})
+        assert result["type"] == "FeatureCollection"
+        assert "timeStamp" in result
+        feat = result["features"][0]
+        assert feat["geometry"]["coordinates"] == [-106.0, 34.0]
+        assert feat["properties"]["trend_category"] == "not enough data"  # single record