four scenarios.

This commit is contained in:
Felix Delattre 2026-02-08 20:59:34 +01:00
parent 46df3be8e7
commit 903bdc2598
6 changed files with 136 additions and 83 deletions

View file

@ -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)

View file

@ -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(

View file

@ -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}")

View file

@ -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)

29
run.py
View file

@ -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:

View file

@ -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);