DOCUMENT CONTENT
文档内容
"""Case 04 core algorithms: NDVI phenology extraction."""
from __future__ import annotations
import numpy as np
from scipy.signal import savgol_filter
def smooth_ndvi(values, window=9):
"""Savitzky–Golay smooth along the time axis (axis 0)."""
x = np.asarray(values, float)
if x.ndim == 0 or x.shape[0] < 3:
raise ValueError("values must have at least three time observations")
w = int(window)
if w % 2 == 0:
w += 1
max_odd = x.shape[0] if x.shape[0] % 2 else x.shape[0] - 1
w = min(w, max_odd)
w = max(3, w)
if w <= 2:
return x.copy()
return savgol_filter(x, w, min(2, w - 1), axis=0, mode="interp")
def extract_phenology(ndvi, fraction=.2):
"""Extract relative-threshold SOS/EOS for one annual curve.
``ndvi`` is ``(time, cell)`` or ``(time,)``. The returned indices are
positions within that annual curve. This function deliberately does not
combine observations from different years.
"""
x = np.asarray(ndvi, float)
if x.ndim not in (1, 2) or x.shape[0] < 2:
raise ValueError("ndvi must be time x cell or a one-dimensional curve")
if not 0 <= float(fraction) <= 1:
raise ValueError("fraction must be between 0 and 1")
was_1d = x.ndim == 1
if was_1d:
x = x[:, None]
base = np.nanmin(x, axis=0)
peak = np.nanmax(x, axis=0)
amplitude = peak - base
threshold = base + float(fraction) * amplitude
valid = np.isfinite(x)
above = valid & (x >= threshold[None, :])
sos = np.full(x.shape[1], -1, dtype=int)
eos = np.full(x.shape[1], -1, dtype=int)
for j in range(x.shape[1]):
inds = np.flatnonzero(above[:, j])
if inds.size:
sos[j] = int(inds[0]); eos[j] = int(inds[-1])
out = {"sos": sos, "eos": eos, "peak": peak, "amplitude": amplitude, "threshold": threshold}
if was_1d:
return {k: np.asarray(v).squeeze() for k, v in out.items()}
return out
def extract_phenology_by_year(ndvi, fraction=.2):
"""Extract phenology independently for every year.
Parameters
----------
ndvi : array-like, shape (year, time, cell)
One complete annual composite curve per year and cell. A 2-D
``(year, time)`` array is accepted for a single cell.
fraction : float
Relative amplitude threshold used for SOS/EOS.
Returns
-------
dict of arrays with shape ``(year, cell)`` (or ``(year,)`` for one cell),
plus ``yearly_mean_sos`` and ``yearly_mean_eos`` convenience series.
"""
x = np.asarray(ndvi, float)
if x.ndim == 2:
x = x[:, :, None]
one_cell = True
elif x.ndim == 3:
one_cell = False
else:
raise ValueError("ndvi must have shape year x time x cell")
results = [extract_phenology(x[i], fraction=fraction) for i in range(x.shape[0])]
out = {k: np.stack([r[k] for r in results], axis=0) for k in results[0]}
out["yearly_mean_sos"] = np.nanmean(np.where(out["sos"] >= 0, out["sos"], np.nan), axis=1)
out["yearly_mean_eos"] = np.nanmean(np.where(out["eos"] >= 0, out["eos"], np.nan), axis=1)
if one_cell:
return {k: (v[:, 0] if np.asarray(v).ndim == 2 and k not in ("yearly_mean_sos", "yearly_mean_eos") else v) for k, v in out.items()}
return out