From 903bdc2598a4b6ad14f1e4053e17756dec97f628 Mon Sep 17 00:00:00 2001
From: Felix Delattre <felix@delattre.de>
Date: Sun, 8 Feb 2026 20:59:34 +0100
Subject: [PATCH] four scenarios.

---
 call_efast.py       | 57 +++++++++++++++++++++++++++----------------
 clouds.py           | 12 ++++-----
 generate_indexes.py | 44 ++++++++++++++++++---------------
 post_process.py     | 18 ++++++++++----
 run.py              | 29 ++++++++--------------
 webapp/index.html   | 59 ++++++++++++++++++++++++++++++++++++---------
 6 files changed, 136 insertions(+), 83 deletions(-)

diff --git a/call_efast.py b/call_efast.py
index acb7a6c..022fd27 100644
--- a/call_efast.py
+++ b/call_efast.py
@@ -36,6 +36,10 @@ def _load_clouds(clouds_file):
     return clouds
 
 
+def _get_base_dir(season, site_name, cleaning_strategy):
+    return Path(f"data/{site_name}/{season}/prepared_{cleaning_strategy}/")
+
+
 def _reproject_raster_to_target(
     src_path,
     dst_path,
@@ -67,11 +71,11 @@ def _reproject_raster_to_target(
             rio_shutil.copy(vrt, dst_path, **profile)
 
 
-def prepare_s2(season, site_position, site_name, date_range=None):
+def prepare_s2(season, site_position, site_name, cleaning_strategy="aggressive", date_range=None):
     s2_dir = Path(f"data/{site_name}/{season}/raw/s2/")
     s3_dir = Path(f"data/{site_name}/{season}/raw/s3/")
-    s2_output_dir = Path(f"data/{site_name}/{season}/prepared/s2/")
-    clouds_file = Path(f"data/{site_name}/{season}/clouds.json")
+    s2_output_dir = _get_base_dir(season, site_name, cleaning_strategy) / "s2"
+    clouds_file = Path(f"data/{site_name}/{season}/clouds_{cleaning_strategy}.json")
 
     clouds = _load_clouds(clouds_file)
     s2_output_dir.mkdir(parents=True, exist_ok=True)
@@ -111,11 +115,12 @@ def prepare_s2(season, site_position, site_name, date_range=None):
     distance_to_clouds(s2_output_dir, ratio=RESOLUTION_RATIO)
 
 
-def prepare_s3(season, site_position, site_name, date_range=None):
+def prepare_s3(season, site_position, site_name, cleaning_strategy="aggressive", date_range=None):
     s3_dir = Path(f"data/{site_name}/{season}/raw/s3/")
-    s2_prepared_dir = Path(f"data/{site_name}/{season}/prepared/s2/")
-    s3_preprocessed_dir = Path(f"data/{site_name}/{season}/prepared/s3/")
-    clouds_file = Path(f"data/{site_name}/{season}/clouds.json")
+    base_dir = _get_base_dir(season, site_name, cleaning_strategy)
+    s2_prepared_dir = base_dir / "s2"
+    s3_preprocessed_dir = base_dir / "s3"
+    clouds_file = Path(f"data/{site_name}/{season}/clouds_{cleaning_strategy}.json")
 
     clouds = _load_clouds(clouds_file)
     s3_preprocessed_dir.mkdir(parents=True, exist_ok=True)
@@ -192,14 +197,14 @@ def prepare_s3(season, site_position, site_name, date_range=None):
     shutil.rmtree(temp_composite_dir)
 
 
-def run_efast(season, site_position, site_name, date_range=None):
+def run_efast(season, site_position, site_name, cleaning_strategy="aggressive", sigma=None, date_range=None):
     lat, lon = site_position
     datetime_range = date_range or f"{season}-01-01/{season}-12-31"
 
-    efast_base_dir = Path(f"data/{site_name}/{season}/prepared/")
+    efast_base_dir = _get_base_dir(season, site_name, cleaning_strategy)
     s2_output_dir = efast_base_dir / "s2"
     s3_output_dir = efast_base_dir / "s3"
-    fusion_output_dir = efast_base_dir / "fusion"
+    fusion_output_dir = efast_base_dir / (f"fusion_sigma{sigma}" if sigma else "fusion")
 
     fusion_output_dir.mkdir(parents=True, exist_ok=True)
     print(f"[EFAST] Starting fusion: {site_name} ({lat:.6f}, {lon:.6f}), {season}")
@@ -215,17 +220,16 @@ def run_efast(season, site_position, site_name, date_range=None):
         date_str = current_date.strftime("%Y%m%d")
         output_file = fusion_output_dir / f"REFL_{date_str}.tif"
         try:
-            efast.fusion(
-                current_date,
-                s3_output_dir,
-                s2_output_dir,
-                fusion_output_dir,
-                product="REFL",
-                max_days=30,
-                date_position=2,
-                minimum_acquisition_importance=0.0,
-                ratio=RESOLUTION_RATIO,
-            )
+            kwargs = {
+                "product": "REFL",
+                "max_days": 30,
+                "date_position": 2,
+                "minimum_acquisition_importance": 0.0,
+                "ratio": RESOLUTION_RATIO,
+            }
+            if sigma is not None:
+                kwargs["sigma"] = sigma
+            efast.fusion(current_date, s3_output_dir, s2_output_dir, fusion_output_dir, **kwargs)
             print(
                 f"[EFAST] Saved: {output_file}"
                 if output_file.exists()
@@ -236,3 +240,14 @@ def run_efast(season, site_position, site_name, date_range=None):
         current_date += timedelta(days=1)
 
     print("[EFAST] Completed")
+
+
+def run_all_efast_scenarios(season, site_position, site_name, sigma_value=30, date_range=None):
+    from clouds import detect_clouds
+    
+    for strategy in ["aggressive", "nonaggressive"]:
+        detect_clouds(season, site_name, cleaning_strategy=strategy)
+        prepare_s2(season, site_position, site_name, cleaning_strategy=strategy, date_range=date_range)
+        prepare_s3(season, site_position, site_name, cleaning_strategy=strategy, date_range=date_range)
+        run_efast(season, site_position, site_name, cleaning_strategy=strategy, sigma=None, date_range=date_range)
+        run_efast(season, site_position, site_name, cleaning_strategy=strategy, sigma=sigma_value, date_range=date_range)
diff --git a/clouds.py b/clouds.py
index 61f1a78..ec6b85c 100644
--- a/clouds.py
+++ b/clouds.py
@@ -3,14 +3,14 @@ from pathlib import Path
 from datetime import datetime
 
 WINDOW_DAYS = 14
-NDVI_THRESHOLD = 0.3
-NDVI_DELTA = 0.15
 MIN_WINDOW_SIZE = 3
+THRESHOLDS = {"aggressive": {"threshold": 0.3, "delta": 0.15}, "nonaggressive": {"threshold": 0.2, "delta": 0.25}}
 
 
-def detect_clouds(season, site_name):
-    output_file = Path(f"data/{site_name}/{season}/clouds.json")
+def detect_clouds(season, site_name, cleaning_strategy="aggressive"):
+    output_file = Path(f"data/{site_name}/{season}/clouds_{cleaning_strategy}.json")
     clouds = {"s2": [], "s3": []}
+    thresholds = THRESHOLDS[cleaning_strategy]
 
     for source in ["s2", "s3"]:
         timeseries_file = Path(
@@ -47,9 +47,9 @@ def detect_clouds(season, site_name):
                 continue
 
             max_ndvi = max(window_ndvi)
-            threshold = max_ndvi - NDVI_DELTA
+            threshold = max_ndvi - thresholds["delta"]
 
-            if entry["ndvi"] < threshold and entry["ndvi"] < NDVI_THRESHOLD:
+            if entry["ndvi"] < threshold and entry["ndvi"] < thresholds["threshold"]:
                 clouds[source].append(entry["filename"])
 
         print(
diff --git a/generate_indexes.py b/generate_indexes.py
index 42b8b34..5bc32fb 100644
--- a/generate_indexes.py
+++ b/generate_indexes.py
@@ -235,16 +235,18 @@ def generate_ndvi_post_process(season, site_position, site_name):
 
 
 def create_ndvi_timeseries_post_process(season, site_position, site_name):
-    for source in ["s2", "s3"]:
-        input_dir = Path(f"data/{site_name}/{season}/processed/{source}/")
-        output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/{source}/")
-        _create_timeseries_for_dir(
-            input_dir, output_dir, site_position, f"POST-PROCESS-{source.upper()}"
-        )
-
-    input_dir = Path(f"data/{site_name}/{season}/processed/fusion/")
-    output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/fusion/")
-    _create_timeseries_for_dir(input_dir, output_dir, site_position, "POST-PROCESS-FUSION")
+    for strategy in ["aggressive", "nonaggressive"]:
+        for sigma in [20, 30]:
+            processed_dir = f"processed_{strategy}_sigma{sigma}"
+            for source in ["s2", "s3"]:
+                input_dir = Path(f"data/{site_name}/{season}/{processed_dir}/{source}/")
+                output_dir = Path(f"data/{site_name}/{season}/{processed_dir}/ndvi/{source}/")
+                _create_timeseries_for_dir(
+                    input_dir, output_dir, site_position, f"POST-PROCESS-{source.upper()}-{strategy}-σ{sigma}"
+                )
+            input_dir = Path(f"data/{site_name}/{season}/{processed_dir}/fusion/")
+            output_dir = Path(f"data/{site_name}/{season}/{processed_dir}/ndvi/fusion/")
+            _create_timeseries_for_dir(input_dir, output_dir, site_position, f"POST-PROCESS-FUSION-{strategy}-σ{sigma}")
 
 
 def _calculate_and_write_gcc(input_file, output_file):
@@ -434,13 +436,15 @@ def generate_gcc_post_process(season, site_position, site_name):
 
 
 def create_gcc_timeseries_post_process(season, site_position, site_name):
-    for source in ["s2", "s3"]:
-        input_dir = Path(f"data/{site_name}/{season}/processed/{source}/")
-        output_dir = Path(f"data/{site_name}/{season}/processed/gcc/{source}/")
-        _create_gcc_timeseries_for_dir(
-            input_dir, output_dir, site_position, f"POST-PROCESS-{source.upper()}"
-        )
-
-    input_dir = Path(f"data/{site_name}/{season}/processed/fusion/")
-    output_dir = Path(f"data/{site_name}/{season}/processed/gcc/fusion/")
-    _create_gcc_timeseries_for_dir(input_dir, output_dir, site_position, "POST-PROCESS-FUSION")
+    for strategy in ["aggressive", "nonaggressive"]:
+        for sigma in [20, 30]:
+            processed_dir = f"processed_{strategy}_sigma{sigma}"
+            for source in ["s2", "s3"]:
+                input_dir = Path(f"data/{site_name}/{season}/{processed_dir}/{source}/")
+                output_dir = Path(f"data/{site_name}/{season}/{processed_dir}/gcc/{source}/")
+                _create_gcc_timeseries_for_dir(
+                    input_dir, output_dir, site_position, f"POST-PROCESS-{source.upper()}-{strategy}-σ{sigma}"
+                )
+            input_dir = Path(f"data/{site_name}/{season}/{processed_dir}/fusion/")
+            output_dir = Path(f"data/{site_name}/{season}/{processed_dir}/gcc/fusion/")
+            _create_gcc_timeseries_for_dir(input_dir, output_dir, site_position, f"POST-PROCESS-FUSION-{strategy}-σ{sigma}")
diff --git a/post_process.py b/post_process.py
index bf7d2af..f2bcc43 100644
--- a/post_process.py
+++ b/post_process.py
@@ -6,20 +6,21 @@ from rasterio.warp import reproject, Resampling
 from rasterio.io import MemoryFile
 
 
-def process_cropped(season, site_position, site_name):
+def process_cropped(season, site_position, site_name, cleaning_strategy="aggressive", sigma=None):
     """Crop fusion to valid data, then crop S2/S3 to match."""
     base = Path(f"data/{site_name}/{season}")
-    prepared = base / "prepared"
-    processed = base / "processed"
+    prepared = base / f"prepared_{cleaning_strategy}"
+    processed_dir = f"processed_{cleaning_strategy}_sigma{sigma}" if sigma else f"processed_{cleaning_strategy}_sigma20"
+    processed = base / processed_dir
     
     s2_prep = prepared / "s2"
     s3_prep = prepared / "s3"
-    fusion_prep = prepared / "fusion"
+    fusion_prep = prepared / (f"fusion_sigma{sigma}" if sigma else "fusion")
     
     for output_dir in [processed / "s2", processed / "s3", processed / "fusion"]:
         output_dir.mkdir(parents=True, exist_ok=True)
     
-    print(f"[PROCESS] Processing files: {site_name}, {season}")
+    print(f"[PROCESS] Processing files: {site_name}, {season}, {cleaning_strategy}, sigma={sigma or 20}")
     
     # Crop fusion to valid data and get dimensions
     fusion_dims = {}
@@ -87,3 +88,10 @@ def process_cropped(season, site_position, site_name):
             print(f"[PROCESS] Cropped: {output_file}")
     
     print("[PROCESS] Completed")
+
+
+def process_all_scenarios(season, site_position, site_name):
+    """Process all 4 EFAST scenarios."""
+    for strategy in ["aggressive", "nonaggressive"]:
+        for sigma in [None, 30]:
+            process_cropped(season, site_position, site_name, cleaning_strategy=strategy, sigma=sigma)
diff --git a/run.py b/run.py
index c34a63f..da8df42 100644
--- a/run.py
+++ b/run.py
@@ -1,5 +1,5 @@
-from call_efast import run_efast, prepare_s2, prepare_s3
-from post_process import process_cropped
+from call_efast import run_all_efast_scenarios
+from post_process import process_all_scenarios
 from generate_indexes import (
     generate_ndvi_raw,
     create_ndvi_timeseries_raw,
@@ -23,26 +23,17 @@ def run_pipeline(season, site_position, site_name):
         #download_phenocam_greenness(season, site_position, site_name)
 
         # print(f"Generating NDVI for raw data: {site_name}, {season}")
-        # generate_ndvi_raw(season, site_position, site_name)
-        create_ndvi_timeseries_raw(season, site_position, site_name)
+        # create_ndvi_timeseries_raw(season, site_position, site_name)
 
-        # print(f"Detecting clouds for {site_name}, {season}")
-        #detect_clouds(season, site_name)
+        # print(f"Running EFAST fusion for all scenarios: {site_name}, {season}")
+        # run_all_efast_scenarios(season, site_position, site_name)
 
-        #print(f"Preparing data for EFAST fusion for {site_name}, {season}")
-        #prepare_s2(season, site_position, site_name)
-        #prepare_s3(season, site_position, site_name)
-
-        #print(f"Running EFAST fusion for {site_name}, {season}")
-        #run_efast(season, site_position, site_name)
-
-        #print(f"Post-processing data: {site_name}, {season}")
-        #process_cropped(season, site_position, site_name)
-        #print(f"Generating NDVI for final outputs: {site_name}, {season}")
-        #generate_ndvi_post_process(season, site_position, site_name)
+        print(f"Post-processing data: {site_name}, {season}")
+        process_all_scenarios(season, site_position, site_name)
+        print(f"Generating NDVI for final outputs: {site_name}, {season}")
         create_ndvi_timeseries_post_process(season, site_position, site_name)
-        #print(f"Generating GCC for final outputs: {site_name}, {season}")
-        generate_gcc_post_process(season, site_position, site_name)
+        print(f"Generating GCC for final outputs: {site_name}, {season}")
+        # generate_gcc_post_process(season, site_position, site_name)  # No-op function
         create_gcc_timeseries_post_process(season, site_position, site_name)
 
     except Exception as e:
diff --git a/webapp/index.html b/webapp/index.html
index 918cc1d..c30da18 100644
--- a/webapp/index.html
+++ b/webapp/index.html
@@ -9,6 +9,8 @@
     <style>
         body { margin: 0; font-family: sans-serif; }
         .slider-container { position: sticky; top: 0; background: white; padding: 20px; z-index: 1000; border-bottom: 1px solid #ccc; }
+        .scenario-selector { margin-bottom: 10px; }
+        .scenario-selector select { padding: 5px 10px; font-size: 14px; }
         .container { max-width: 1400px; margin: 0 auto; padding: 20px; }
         .header { display: flex; gap: 20px; margin-bottom: 20px; border-bottom: 1px solid #ccc; padding-top: 10px;padding-bottom: 20px;}
         .header-col { flex: 1; }
@@ -58,6 +60,15 @@
                 <div id="sitemap" class="sitemap"></div>
             </div>
         </div>
+        <div class="scenario-selector">
+            <label for="scenarioSelect">Scenario: </label>
+            <select id="scenarioSelect">
+                <option value="aggressive_20">Aggressive (σ=20)</option>
+                <option value="aggressive_30">Aggressive (σ=30)</option>
+                <option value="nonaggressive_20">Non-aggressive (σ=20)</option>
+                <option value="nonaggressive_30">Non-aggressive (σ=30)</option>
+            </select>
+        </div>
         <div class="maps">
             <div class="map-container">
                 <h3>S2</h3>
@@ -106,6 +117,16 @@
         const osmOpts = { attribution: "OpenStreetMap", opacity: 0.4 };
         const mapOpts = { zoomControl: false };
         const sitePosition = [47.116171, 11.320308];
+        const siteName = "innsbruck";
+        const season = "2024";
+        
+        const urlParams = new URLSearchParams(location.search);
+        const strategy = urlParams.get("strategy") || "aggressive";
+        const sigma = urlParams.get("sigma") || "20";
+        
+        function getFusionPath() {
+            return `processed_${strategy}_sigma${sigma}`;
+        }
         const siteMap = L.map("sitemap", { zoomControl: false }).setView(sitePosition, 4).addLayer(L.tileLayer(osmUrl, { attribution: "OpenStreetMap", opacity: 1 }));
         L.marker(sitePosition, { icon: L.divIcon({ className: "site-marker", html: "<div style='width:8px;height:8px;background:red;border:2px solid white;border-radius:50%;box-shadow:0 0 2px rgba(0,0,0,0.5);'></div>", iconSize: [8, 8] }) }).addTo(siteMap);
         const maps = {
@@ -141,14 +162,15 @@
         });
         
         async function loadTimeseries() {
+            const fusionPath = getFusionPath();
             const [s2, fusion, s3, s2gcc, fusiongcc, s3gcc, phenocam] = await Promise.all([
-                fetch("../data/innsbruck/2024/processed/ndvi/s2/timeseries.json").then(r => r.json()),
-                fetch("../data/innsbruck/2024/processed/ndvi/fusion/timeseries.json").then(r => r.json()).catch(() => []),
-                fetch("../data/innsbruck/2024/processed/ndvi/s3/timeseries.json").then(r => r.json()),
-                fetch("../data/innsbruck/2024/processed/gcc/s2/timeseries.json").then(r => r.json()).catch(() => []),
-                fetch("../data/innsbruck/2024/processed/gcc/fusion/timeseries.json").then(r => r.json()).catch(() => []),
-                fetch("../data/innsbruck/2024/processed/gcc/s3/timeseries.json").then(r => r.json()).catch(() => []),
-                fetch("../data/innsbruck/2024/raw/phenocam/timeseries.json").then(r => r.json()).catch(() => [])
+                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/timeseries.json`).then(r => r.json()).catch(() => [])
             ]);
             timeseries = { s2, fusion, s3 };
             greennessTimeseries = { s2: s2gcc, fusion: fusiongcc, s3: s3gcc };
@@ -487,6 +509,7 @@
         
         async function findFile(dateStr, source) {
             const target = new Date(dateStr);
+            const basePath = source === "fusion" ? getFusionPath() : `processed_${strategy}_sigma${sigma}`;
             // Search outward from target date (0, ±1, ±2, ±3, ...) until we find the closest file
             // Check dates in order: exact, then -1, +1, then -2, +2, etc.
             // Limit to ±365 days to avoid infinite search
@@ -496,7 +519,7 @@
                     const date = target.toISOString().split("T")[0].replace(/-/g, "");
                     const filename = `${date}_0.geotiff`;
                     try {
-                        const res = await fetch(`../data/innsbruck/2024/processed/${source}/${filename}`, { method: 'HEAD' });
+                        const res = await fetch(`../data/${siteName}/${season}/${basePath}/${source}/${filename}`, { method: 'HEAD' });
                         if (res.ok) return filename;
                     } catch {}
                 } else {
@@ -507,7 +530,7 @@
                         const date = d.toISOString().split("T")[0].replace(/-/g, "");
                         const filename = `${date}_0.geotiff`;
                         try {
-                            const res = await fetch(`../data/innsbruck/2024/processed/${source}/${filename}`, { method: 'HEAD' });
+                            const res = await fetch(`../data/${siteName}/${season}/${basePath}/${source}/${filename}`, { method: 'HEAD' });
                             if (res.ok) return filename;
                         } catch {}
                     }
@@ -522,7 +545,8 @@
         }
         
         async function loadGeotiff(source, filename) {
-            const path = `../data/innsbruck/2024/processed/${source}/${filename}`;
+            const basePath = source === "fusion" ? getFusionPath() : `processed_${strategy}_sigma${sigma}`;
+            const path = `../data/${siteName}/${season}/${basePath}/${source}/${filename}`;
             const tiff = await GeoTIFF.fromArrayBuffer(await (await fetch(path)).arrayBuffer());
             const image = await tiff.getImage();
             const rasters = await image.readRasters();
@@ -587,7 +611,7 @@
                     const d = new Date(target);
                     d.setDate(d.getDate() + offset * dir);
                     const date = d.toISOString().split("T")[0].replace(/-/g, "");
-                    const url = `../data/innsbruck/2024/raw/phenocam/${date}.jpg`;
+                    const url = `../data/${siteName}/${season}/raw/phenocam/${date}.jpg`;
                     try {
                         const res = await fetch(url, { method: 'HEAD' });
                         if (res.ok) {
@@ -607,6 +631,8 @@
             dateDisplay.textContent = date;
             const params = new URLSearchParams();
             params.set("date", date);
+            params.set("strategy", strategy);
+            if (sigma !== "20") params.set("sigma", sigma);
             history.replaceState({}, "", `?${params}`);
             drawTimeseries();
             drawGreennessTimeseries();
@@ -625,7 +651,16 @@
             await loadPhenoCam(date);
         }
         
-        const urlParams = new URLSearchParams(location.search);
+        const scenarioSelect = document.getElementById("scenarioSelect");
+        scenarioSelect.value = `${strategy}_${sigma}`;
+        scenarioSelect.addEventListener("change", function() {
+            const [newStrategy, newSigma] = this.value.split("_");
+            const params = new URLSearchParams(location.search);
+            params.set("strategy", newStrategy);
+            params.set("sigma", newSigma);
+            window.location.search = params.toString();
+        });
+        
         const urlDate = urlParams.get("date");
         if (urlDate) slider.value = daysFromDate(urlDate);
         slider.addEventListener("input", updateImages);