diff --git a/metrics_stats.py b/metrics_stats.py index df63002..e22bd57 100644 --- a/metrics_stats.py +++ b/metrics_stats.py @@ -148,21 +148,46 @@ def calculate_phenocam_stats(phenocam_ts): } -def _s2_kept_date_set(base: Path, strategy: str) -> set: +def _s2_gcc_series_from_preselection(base: Path): + """Build the raw S2 GCC series from s2_preselection.json. + + Uses the 3x3 site-window band means stored per raw S2 acquisition and + computes GCC = b03 / (b02 + b03 + b04). Scale cancels, so DN vs + reflectance is irrelevant. Returns (all_gcc, flags) where all_gcc maps + YYYY-MM-DD -> gcc for every row with a positive band sum, and flags maps + the same date key -> (excluded_aggressive, excluded_nonaggressive). + """ path = base / "raw" / "preselection" / "s2_preselection.json" if not path.exists(): - return set() + return {}, {} with open(path) as f: rows = json.load(f) - key = f"excluded_{strategy}" - out = set() + all_gcc: dict = {} + flags: dict = {} for e in rows: - if e.get(key): - continue nk = _norm_date_key(e.get("date")) - if nk: - out.add(nk) - return out + if not nk: + continue + try: + b02 = float(e.get("b02")) + b03 = float(e.get("b03")) + b04 = float(e.get("b04")) + except (TypeError, ValueError): + continue + total = b02 + b03 + b04 + if not np.isfinite(total) or total <= 0: + continue + gcc = b03 / total + if not np.isfinite(gcc): + continue + if nk in all_gcc: + continue + all_gcc[nk] = float(gcc) + flags[nk] = ( + bool(e.get("excluded_aggressive")), + bool(e.get("excluded_nonaggressive")), + ) + return all_gcc, flags def _whittaker_smooth_dict(obs_dates, obs_values, lam: float, n_min: int = 3): @@ -221,33 +246,27 @@ def calculate_all_metrics(season, site_name, site_position): results["phenocam_stats"] = phenocam_stats baseline = {} - s2_ts = {} - for sub in ("processed_aggressive_sigma20", "processed_nonaggressive_sigma20"): - p = base / sub / "gcc" / "s2" / "timeseries.json" - if p.exists(): - s2_ts = load_timeseries(p) - if s2_ts: - break - if s2_ts: - m0 = calculate_temporal_metrics(s2_ts, phenocam_ts) + all_gcc, flags = _s2_gcc_series_from_preselection(base) + if all_gcc: + m0 = calculate_temporal_metrics(all_gcc, phenocam_ts) if m0: baseline["s2"] = m0 - for strategy in ("aggressive", "nonaggressive"): - kept = _s2_kept_date_set(base, strategy) - filtered = [ - (k, v) - for k, v in sorted( - s2_ts.items(), key=lambda x: _norm_date_key(x[0]) or "" - ) - if _norm_date_key(k) in kept and v is not None - ] - if not filtered: + for strategy, flag_idx in (("aggressive", 0), ("nonaggressive", 1)): + kept_items = sorted( + ( + (d, g) + for d, g in all_gcc.items() + if d in flags and not flags[d][flag_idx] + ), + key=lambda x: x[0], + ) + if not kept_items: continue - sub_ts = dict(filtered) - mcf = calculate_temporal_metrics(sub_ts, phenocam_ts) + kept_ts = dict(kept_items) + mcf = calculate_temporal_metrics(kept_ts, phenocam_ts) if mcf: baseline.setdefault("s2_cloudfree", {})[strategy] = mcf - obs_d, obs_v = zip(*filtered) + obs_d, obs_v = zip(*kept_items) smooth = _whittaker_smooth_dict(obs_d, obs_v, WHITTAKER_LAMBDA_DAYS_SQ) if smooth: mw = calculate_temporal_metrics(smooth, phenocam_ts) diff --git a/requirements.txt b/requirements.txt index c04d625..dfb227a 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,6 +4,7 @@ openeo python-dotenv netCDF4 numpy +timesat requests scipy matplotlib diff --git a/run.py b/run.py index f590008..73643a0 100644 --- a/run.py +++ b/run.py @@ -17,6 +17,7 @@ from preparation import ( ) from metrics_indices import create_prepared_fusion_timeseries from metrics_stats import calculate_all_metrics +from phenology_timesat import write_phenocam_phenology_for_site def run_pipeline(season, site_position, site_name): @@ -27,6 +28,9 @@ def run_pipeline(season, site_position, site_name): # download_s3(season, site_position, site_name) # download_phenocam(season, site_position, site_name) + print(f"PhenoCam phenology (50 % amplitude): {site_name}, {season}") + write_phenocam_phenology_for_site(site_name, season) + print(f"Creating preselection timeseries: {site_name}, {season}") create_timeseries(season, site_position, site_name) diff --git a/webapp/fusion.html b/webapp/fusion.html index d3e5879..205cba6 100644 --- a/webapp/fusion.html +++ b/webapp/fusion.html @@ -44,6 +44,7 @@ Fusion Postprocessed Metrics + Phenology

Innsbruck

2024

diff --git a/webapp/index.html b/webapp/index.html index f5adc0d..aec4edf 100644 --- a/webapp/index.html +++ b/webapp/index.html @@ -25,6 +25,7 @@ .site-info h2 { margin: 0 0 20px 0; font-size: 18px; color: #666; } .phenocam-label { font-size: 12px; margin-bottom: 5px; color: #666; } .phenocam-date { font-size: 11px; margin-top: 5px; color: #999; } + .timeseries-phenology { font-size: 11px; color: #444; line-height: 1.45; margin-top: 6px; white-space: pre-line; } .phenocam-image { width: 100%; height: 200px; object-fit: contain; border: 1px solid #ccc; } .sitemap { height: 200px; border: 1px solid #ccc; margin-top: 32px; } #greennesstimeseries { margin-top: 0; } @@ -54,6 +55,7 @@ Fusion Postprocessed Metrics + Phenology
@@ -70,6 +72,7 @@
Greenness Index Timeseries
+
@@ -172,6 +175,8 @@ let timeseries = { s2: [], fusion: [], s3: [] }; let greennessTimeseries = { s2: [], fusion: [], s3: [] }; let phenocamGreennessTimeseries = []; + /** @type {{ green_up_50pct_date: string | null, green_down_50pct_date: string | null } | null} */ + let phenocamPhenology = null; // Add site marker to all maps for (const source of ["s2", "fusion", "s3"]) { @@ -180,6 +185,39 @@ let syncing = false; const allMaps = Object.values(maps); + function updatePhenocamPhenologyText() { + const el = document.getElementById("phenocamPhenology"); + if (!el) return; + if (!phenocamPhenology || (phenocamPhenology.green_up_50pct_date == null && phenocamPhenology.green_down_50pct_date == null)) { + el.textContent = ""; + return; + } + const up = phenocamPhenology.green_up_50pct_date || "—"; + const dn = phenocamPhenology.green_down_50pct_date || "—"; + el.textContent = `* Green-up: ${up}\n* Green-down: ${dn}`; + } + + const PHENOLOGY_VLINE = "rgba(150, 200, 255, 0.95)"; + + function strokePhenologyVlines(ctx, minDate, maxDate, dateRange, plotW, pad, plotH) { + if (!phenocamPhenology) return; + for (const iso of [phenocamPhenology.green_up_50pct_date, phenocamPhenology.green_down_50pct_date]) { + if (!iso) continue; + const d = new Date(iso); + if (isNaN(d.getTime()) || d < minDate || d > maxDate) continue; + const px = pad + ((d - minDate) / dateRange) * plotW; + ctx.save(); + ctx.setLineDash([4, 3]); + ctx.strokeStyle = PHENOLOGY_VLINE; + ctx.lineWidth = 1.5; + ctx.beginPath(); + ctx.moveTo(px, pad); + ctx.lineTo(px, pad + plotH); + ctx.stroke(); + ctx.restore(); + } + } + const syncMaps = (sourceMap) => { if (syncing) return; syncing = true; @@ -197,18 +235,21 @@ async function loadTimeseries() { metricsData = null; const fusionPath = getFusionPath(); - const [s2, fusion, s3, s2gcc, fusiongcc, s3gcc, phenocam] = await Promise.all([ + const [s2, fusion, s3, s2gcc, fusiongcc, s3gcc, phenocam, phenPhen] = await Promise.all([ fetch(`data/${siteName}/${season}/processed_${strategy}_sigma${sigma}/ndvi/s2/timeseries.json`).then(r => r.json()).catch(() => []), fetch(`data/${siteName}/${season}/${fusionPath}/ndvi/fusion/timeseries.json`).then(r => r.json()).catch(() => []), fetch(`data/${siteName}/${season}/processed_${strategy}_sigma${sigma}/ndvi/s3/timeseries.json`).then(r => r.json()).catch(() => []), fetch(`data/${siteName}/${season}/processed_${strategy}_sigma${sigma}/gcc/s2/timeseries.json`).then(r => r.json()).catch(() => []), fetch(`data/${siteName}/${season}/${fusionPath}/gcc/fusion/timeseries.json`).then(r => r.json()).catch(() => []), fetch(`data/${siteName}/${season}/processed_${strategy}_sigma${sigma}/gcc/s3/timeseries.json`).then(r => r.json()).catch(() => []), - fetch(`data/${siteName}/${season}/raw/phenocam/phenocam_gcc.json`).then(r => r.json()).catch(() => []) + fetch(`data/${siteName}/${season}/raw/phenocam/phenocam_gcc.json`).then(r => r.json()).catch(() => []), + fetch(`data/${siteName}/${season}/raw/phenocam/phenocam_phenology.json`).then(r => r.ok ? r.json() : null).catch(() => null) ]); timeseries = { s2, fusion, s3 }; greennessTimeseries = { s2: s2gcc, fusion: fusiongcc, s3: s3gcc }; phenocamGreennessTimeseries = phenocam; + phenocamPhenology = phenPhen && (phenPhen.green_up_50pct_date != null || phenPhen.green_down_50pct_date != null) ? phenPhen : null; + updatePhenocamPhenologyText(); // Load all scenario GCC timeseries for comparison const scenarios = [ @@ -289,10 +330,11 @@ ctx.fillStyle = "#888"; const axisY = pad + plotH; for (const t of data) ctx.fillRect(x(t.date) - 1, axisY - 1, 2, 2); - + strokePhenologyVlines(ctx, minDate, maxDate, dateRange, plotW, pad, plotH); const xPos = x(currentDate); ctx.strokeStyle = "#f00"; ctx.lineWidth = 2; + ctx.setLineDash([]); ctx.beginPath(); ctx.moveTo(xPos, pad); ctx.lineTo(xPos, pad + plotH); @@ -360,11 +402,12 @@ ctx.fillStyle = "#888"; const axisY = pad + plotH; for (const t of data) ctx.fillRect(x(t.date) - 1, axisY - 1, 2, 2); - + strokePhenologyVlines(ctx, minDate, maxDate, dateRange, plotW, pad, plotH); const currentDate = dateFromDays(parseInt(slider.value)); const xPos = x(currentDate); ctx.strokeStyle = "#f00"; ctx.lineWidth = 2; + ctx.setLineDash([]); ctx.beginPath(); ctx.moveTo(xPos, pad); ctx.lineTo(xPos, pad + plotH); @@ -460,17 +503,17 @@ ctx.fillStyle = "#888"; const axisY = pad + plotH; for (const t of [...s2data, ...fusiondata, ...phenocamdata]) ctx.fillRect(x(t.date) - 1, axisY - 1, 2, 2); - + strokePhenologyVlines(ctx, minDate, maxDate, dateRange, plotW, pad, plotH); const legendX = pad + plotW - 80, legendY = pad + 5; ctx.font = "9px sans-serif"; if (s2data.length) { ctx.strokeStyle = "#ff6600"; ctx.beginPath(); ctx.moveTo(legendX, legendY); ctx.lineTo(legendX + 15, legendY); ctx.stroke(); ctx.fillStyle = "#000"; ctx.fillText("S2", legendX + 18, legendY + 3); } if (fusiondata.length) { ctx.strokeStyle = "#9900cc"; ctx.beginPath(); ctx.moveTo(legendX, legendY + 12); ctx.lineTo(legendX + 15, legendY + 12); ctx.stroke(); ctx.fillStyle = "#000"; ctx.fillText("Fusion", legendX + 18, legendY + 15); } if (phenocamdata.length) { ctx.strokeStyle = "#00aa00"; ctx.beginPath(); ctx.moveTo(legendX, legendY + 24); ctx.lineTo(legendX + 15, legendY + 24); ctx.stroke(); ctx.fillStyle = "#000"; ctx.fillText("PhenoCam", legendX + 18, legendY + 27); } - const currentDate = dateFromDays(parseInt(slider.value)); const xPos = x(currentDate); ctx.strokeStyle = "#f00"; ctx.lineWidth = 2; + ctx.setLineDash([]); ctx.beginPath(); ctx.moveTo(xPos, pad); ctx.lineTo(xPos, pad + plotH); @@ -545,6 +588,7 @@ ctx.fillText(name, legendX + 18, legendY + i * 12 + 3); } }); + strokePhenologyVlines(ctx, minDate, maxDate, dateRange, plotW, pad, plotH); } function drawMetricsTable() { diff --git a/webapp/metrics.html b/webapp/metrics.html index c193e05..9cc23b4 100644 --- a/webapp/metrics.html +++ b/webapp/metrics.html @@ -38,6 +38,7 @@ Fusion Postprocessed Metrics + Phenology

Metrics

diff --git a/webapp/postprocessed.html b/webapp/postprocessed.html index 98d9c61..f38ad70 100644 --- a/webapp/postprocessed.html +++ b/webapp/postprocessed.html @@ -44,6 +44,7 @@ Fusion Postprocessed Metrics + Phenology

Innsbruck

2024

diff --git a/webapp/prepared.html b/webapp/prepared.html index 16fe818..e8336d9 100644 --- a/webapp/prepared.html +++ b/webapp/prepared.html @@ -44,6 +44,7 @@ Fusion Postprocessed Metrics + Phenology

Innsbruck

2024

diff --git a/webapp/preselection.html b/webapp/preselection.html index 8c39edb..fdb171e 100644 --- a/webapp/preselection.html +++ b/webapp/preselection.html @@ -43,6 +43,7 @@ Fusion Postprocessed Metrics + Phenology

Innsbruck

2024