|
| 1 | +# =============================================================================== |
| 2 | +# Trend statistics for DIE products. |
| 3 | +# |
| 4 | +# Pure analysis (no serialization / I/O): daily aggregation, the qualification |
| 5 | +# gate, and the Mann-Kendall + Theil-Sen trend test. Kept separate from |
| 6 | +# backend/persisters/ogc_features.py (which only builds GeoJSON) because this is |
| 7 | +# statistics, is independently testable, and pulls heavier deps (scipy, |
| 8 | +# pymannkendall) lazily. |
| 9 | +# =============================================================================== |
| 10 | +from datetime import datetime, timezone |
| 11 | +from typing import Optional |
| 12 | + |
| 13 | +# Seconds per Julian year (365.25 days) — used to express the Theil-Sen slope |
| 14 | +# per year. |
| 15 | +_SECONDS_PER_YEAR = 31557600.0 |
| 16 | + |
| 17 | +# A well is classified only when it has enough data for a meaningful |
| 18 | +# Mann-Kendall test: at least 10 daily points, or at least 4 spanning >= 2 years. |
| 19 | +_TREND_MIN_RECORDS = 10 |
| 20 | +_TREND_MIN_RECORDS_WITH_SPAN = 4 |
| 21 | +_TREND_MIN_SPAN_YEARS = 2.0 |
| 22 | +_TREND_ALPHA = 0.05 # significance level for the Mann-Kendall test |
| 23 | + |
| 24 | +_MK_COMMON = ( |
| 25 | + "Monotonic trend is tested with the non-parametric Mann-Kendall test " |
| 26 | + "(pymannkendall.original_test, alpha=0.05) on the daily series ordered by " |
| 27 | + "time; the rate is the Theil-Sen slope vs time. A well is classified only " |
| 28 | + "when it has at least 10 daily points, or at least 4 daily points spanning " |
| 29 | + "at least 2 years; otherwise 'not enough data'." |
| 30 | +) |
| 31 | + |
| 32 | +# Human-readable description of each trend product's method, embedded in the |
| 33 | +# collection so consumers know how the classification was derived. |
| 34 | +TREND_METHOD_DESCRIPTION = ( |
| 35 | + "Depth-to-water trend per well. Observations are downsampled to one point " |
| 36 | + "per calendar day, keeping the daily MINIMUM depth-to-water (shallowest " |
| 37 | + f"reading). {_MK_COMMON} A significant increasing slope is 'increasing' " |
| 38 | + "(water level getting DEEPER, i.e. a declining water table), a significant " |
| 39 | + "decreasing slope is 'decreasing' (water level getting SHALLOWER), else " |
| 40 | + "'stable'. Mirrors the Ocotillo API ogc_depth_to_water_trend_wells " |
| 41 | + "materialized view, using Mann-Kendall + Theil-Sen instead of OLS." |
| 42 | +) |
| 43 | + |
| 44 | +ANALYTE_TREND_METHOD_DESCRIPTION = ( |
| 45 | + "Analyte concentration trend per well. Observations are downsampled to one " |
| 46 | + f"point per calendar day, keeping the daily MEAN concentration. {_MK_COMMON} " |
| 47 | + "A significant increasing slope is 'increasing' (concentration rising), a " |
| 48 | + "significant decreasing slope is 'decreasing', else 'stable'." |
| 49 | +) |
| 50 | + |
| 51 | + |
| 52 | +def parse_epoch_seconds(date, time) -> Optional[float]: |
| 53 | + """Best-effort parse of a DIE date (+optional time) to POSIX seconds (UTC).""" |
| 54 | + if not date: |
| 55 | + return None |
| 56 | + text = f"{date}T{time}" if time else str(date) |
| 57 | + text = text.replace("Z", "") |
| 58 | + try: |
| 59 | + dt = datetime.fromisoformat(text) |
| 60 | + if dt.tzinfo is None: |
| 61 | + dt = dt.replace(tzinfo=timezone.utc) |
| 62 | + return dt.timestamp() |
| 63 | + except (ValueError, TypeError): |
| 64 | + pass |
| 65 | + try: |
| 66 | + dt = datetime.fromisoformat(str(date)).replace(tzinfo=timezone.utc) |
| 67 | + return dt.timestamp() |
| 68 | + except (ValueError, TypeError): |
| 69 | + return None |
| 70 | + |
| 71 | + |
| 72 | +def daily_series(obs_list: list, reducer: str = "min") -> tuple[int, list]: |
| 73 | + """Reduce a well's observations to one point per calendar day, keyed at the |
| 74 | + day's UTC midnight epoch. *reducer* selects the daily aggregate: "min" |
| 75 | + (shallowest depth-to-water), "max", or "mean" (e.g. analyte concentration). |
| 76 | +
|
| 77 | + *obs_list* is a list of observation payload dicts (parameter_value, |
| 78 | + date_measured, time_measured). Operating on dicts avoids rebuilding |
| 79 | + ParameterRecord objects for what can be millions of observations. |
| 80 | +
|
| 81 | + Downsampling bounds the O(n^2) Mann-Kendall cost for high-frequency wells |
| 82 | + (e.g. continuous loggers) and removes within-day sampling noise. Returns |
| 83 | + (raw_observation_count, [(day_epoch_seconds, value), ...] sorted by day). |
| 84 | + """ |
| 85 | + raw_count = 0 |
| 86 | + buckets: dict = {} # date -> (day_epoch, [values]) |
| 87 | + for obs in obs_list: |
| 88 | + value = obs.get("parameter_value") |
| 89 | + epoch = parse_epoch_seconds(obs.get("date_measured"), obs.get("time_measured")) |
| 90 | + if value is None or epoch is None: |
| 91 | + continue |
| 92 | + try: |
| 93 | + v = float(value) |
| 94 | + except (TypeError, ValueError): |
| 95 | + continue |
| 96 | + raw_count += 1 |
| 97 | + day = datetime.fromtimestamp(epoch, tz=timezone.utc).date() |
| 98 | + day_epoch = datetime( |
| 99 | + day.year, day.month, day.day, tzinfo=timezone.utc |
| 100 | + ).timestamp() |
| 101 | + if day not in buckets: |
| 102 | + buckets[day] = (day_epoch, []) |
| 103 | + buckets[day][1].append(v) |
| 104 | + |
| 105 | + reduce_fn = { |
| 106 | + "min": min, |
| 107 | + "max": max, |
| 108 | + "mean": lambda vs: sum(vs) / len(vs), |
| 109 | + }[reducer] |
| 110 | + |
| 111 | + pairs = sorted( |
| 112 | + ((day_epoch, reduce_fn(vals)) for day_epoch, vals in buckets.values()), |
| 113 | + key=lambda p: p[0], |
| 114 | + ) |
| 115 | + return raw_count, pairs |
| 116 | + |
| 117 | + |
| 118 | +def qualifies_for_trend(record_count, span_years) -> bool: |
| 119 | + return record_count >= _TREND_MIN_RECORDS or ( |
| 120 | + record_count >= _TREND_MIN_RECORDS_WITH_SPAN |
| 121 | + and span_years >= _TREND_MIN_SPAN_YEARS |
| 122 | + ) |
| 123 | + |
| 124 | + |
| 125 | +def mann_kendall_trend(years: list, values: list): |
| 126 | + """Run the Mann-Kendall trend test + Theil-Sen slope. |
| 127 | +
|
| 128 | + Returns (trend_category, slope_per_year, p_value, tau). *years* are decimal |
| 129 | + years, *values* the measured quantity, both ordered by time. trend_category |
| 130 | + is one of 'increasing' / 'decreasing' / 'stable'. |
| 131 | + """ |
| 132 | + import pymannkendall as mk |
| 133 | + from scipy.stats import theilslopes |
| 134 | + |
| 135 | + result = mk.original_test(values, alpha=_TREND_ALPHA) |
| 136 | + # Time-aware Theil-Sen slope (per year) — robust and correct for the |
| 137 | + # irregular sampling typical of these records, unlike MK's index-based slope |
| 138 | + # which assumes unit spacing. |
| 139 | + slope_per_year = float(theilslopes(values, years)[0]) |
| 140 | + |
| 141 | + # mk trend is 'increasing' / 'decreasing' / 'no trend'. |
| 142 | + category = "stable" if result.trend == "no trend" else result.trend |
| 143 | + return category, slope_per_year, float(result.p), float(result.Tau) |
0 commit comments