Added post-processing to align area size.

call_efast.py

@@ -155,13 +155,13 @@ def prepare_s3(season, site_position, site_name, date_range=None):
         raise ValueError(f"No REFL files found in {s2_prepared_dir}")
 
     # Get bounds from REFL file (full coverage, matches S2)
+    # Use integer division to match distance_to_clouds logic exactly
     with rasterio.open(sen2_ref_paths[0]) as s2_ref:
         target_bounds = s2_ref.bounds
         target_crs = s2_ref.crs
-        s2_resolution = abs(s2_ref.transform[0])
-        s3_resolution = s2_resolution * RESOLUTION_RATIO
-        width = int((target_bounds.right - target_bounds.left) / s3_resolution)
-        height = int((target_bounds.top - target_bounds.bottom) / s3_resolution)
+        # Use integer division matching distance_to_clouds: s2_height // ratio, s2_width // ratio
+        width = s2_ref.width // RESOLUTION_RATIO
+        height = s2_ref.height // RESOLUTION_RATIO
         s3_transform = rasterio.transform.from_bounds(
             target_bounds.left,
             target_bounds.bottom,

ndvi.py

@@ -190,36 +190,35 @@ def _fusion_namer(f):
     return f"{date_str}_ndvi.geotiff"
 
 
-def generate_ndvi_prepared(season, site_position, site_name):
+def generate_ndvi_post_process(season, site_position, site_name):
     for source in ["s2", "s3"]:
-        input_dir = Path(f"data/{site_name}/{season}/prepared/{source}/")
-        output_dir = Path(f"data/{site_name}/{season}/prepared/ndvi/{source}/")
-        for pattern in ["*.geotiff", "*.tif"]:
+        input_dir = Path(f"data/{site_name}/{season}/processed/{source}/")
+        output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/{source}/")
         _process_ndvi_files(
             input_dir,
             output_dir,
-                f"PREPARED-{source.upper()}",
-                pattern=pattern,
-                output_namer=_get_output_name_prepared,
+            f"POST-PROCESS-{source.upper()}",
+            pattern="*.geotiff",
+            output_namer=lambda f: f.name.replace(".geotiff", "_ndvi.geotiff"),
         )
 
-    input_dir = Path(f"data/{site_name}/{season}/prepared/fusion/")
-    output_dir = Path(f"data/{site_name}/{season}/prepared/ndvi/fusion/")
+    input_dir = Path(f"data/{site_name}/{season}/processed/fusion/")
+    output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/fusion/")
     _process_ndvi_files(
         input_dir,
         output_dir,
-        "FUSION",
-        pattern="REFL_*.tif",
-        output_namer=_fusion_namer,
+        "POST-PROCESS-FUSION",
+        pattern="*.geotiff",
+        output_namer=lambda f: f.name.replace(".geotiff", "_ndvi.geotiff"),
     )
 
 
-def create_ndvi_timeseries_prepared(season, site_position, site_name):
+def create_ndvi_timeseries_post_process(season, site_position, site_name):
     for source in ["s2", "s3"]:
-        output_dir = Path(f"data/{site_name}/{season}/prepared/ndvi/{source}/")
+        output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/{source}/")
         _create_timeseries_for_dir(
-            output_dir, site_position, f"PREPARED-{source.upper()}"
+            output_dir, site_position, f"POST-PROCESS-{source.upper()}"
         )
 
-    output_dir = Path(f"data/{site_name}/{season}/prepared/ndvi/fusion/")
-    _create_timeseries_for_dir(output_dir, site_position, "FUSION")
+    output_dir = Path(f"data/{site_name}/{season}/processed/ndvi/fusion/")
+    _create_timeseries_for_dir(output_dir, site_position, "POST-PROCESS-FUSION")

post_process.py

@@ -0,0 +1,126 @@
+from pathlib import Path
+import numpy as np
+import rasterio
+from rasterio import windows
+
+
+def _crop_to_bounds(src_file, bounds, output_file, row_based_height=None):
+    """Crop a raster file to given bounds and save."""
+    crop_left, crop_bottom, crop_right, crop_top, crop_crs = bounds
+    
+    with rasterio.open(src_file) as src:
+        # Calculate window from bounds
+        window = windows.from_bounds(crop_left, crop_bottom, crop_right, crop_top, src.transform)
+        
+        # Use row-based height if provided (for fusion), otherwise calculate from bounds
+        if row_based_height is not None:
+            col_off = int(round(window.col_off))
+            window = windows.Window(col_off, 0, src.width, row_based_height)
+            # Calculate bottom Y from row index
+            bottom_y = src.transform[5] + row_based_height * src.transform[4]
+        else:
+            pixel_size = abs(src.transform[0])
+            width = int(round((crop_right - crop_left) / pixel_size))
+            height = int(round((crop_top - crop_bottom) / pixel_size))
+            window = windows.Window(
+                int(round(window.col_off)), int(round(window.row_off)), width, height
+            )
+            bottom_y = crop_bottom
+        
+        # Clip window to source bounds
+        src_window = windows.Window(0, 0, src.width, src.height)
+        window = window.intersection(src_window)
+        if not window or window.height <= 0 or window.width <= 0:
+            return False
+        
+        data = src.read(window=window)
+        transform = rasterio.transform.from_bounds(
+            crop_left, bottom_y, crop_right, crop_top, window.width, window.height
+        )
+        
+        profile = src.profile.copy()
+        profile.update({
+            "height": window.height,
+            "width": window.width,
+            "transform": transform,
+            "crs": crop_crs,
+        })
+    
+    with rasterio.open(output_file, "w", **profile) as dst:
+        dst.write(data)
+    return True
+
+
+def process_cropped(season, site_position, site_name):
+    """Crop prepared S2, S3, and fusion files to fusion valid data bounds."""
+    base = Path(f"data/{site_name}/{season}")
+    prepared = base / "prepared"
+    processed = base / "processed"
+    
+    s2_prep = prepared / "s2"
+    s3_prep = prepared / "s3"
+    fusion_prep = prepared / "fusion"
+    
+    for output_dir in [processed / "s2", processed / "s3", processed / "fusion"]:
+        output_dir.mkdir(parents=True, exist_ok=True)
+    
+    print(f"[PROCESS] Cropping files to fusion valid data bounds: {site_name}, {season}")
+    
+    # Determine valid bounds for each fusion file
+    fusion_bounds = {}
+    fusion_rows = {}
+    
+    for fusion_file in fusion_prep.glob("REFL_*.tif"):
+        date_str = fusion_file.stem.split("_")[1]
+        dist_cloud = s2_prep / f"S2A_MSIL2A_{date_str}_DIST_CLOUD.tif"
+        if not dist_cloud.exists():
+            continue
+        
+        with rasterio.open(dist_cloud) as dist_src:
+            dist_bounds = dist_src.bounds
+            dist_crs = dist_src.crs
+        
+        # Find first valid row from bottom in fusion file
+        with rasterio.open(fusion_file) as fusion_src:
+            data = fusion_src.read()
+            height = data.shape[1]
+            
+            first_valid_row = height
+            for row_idx in range(height - 1, -1, -1):
+                if np.any(~np.isnan(data[:, row_idx, :]) & (data[:, row_idx, :] > 0.001)):
+                    first_valid_row = row_idx
+                    break
+            
+            valid_bottom_y = (fusion_src.transform * (0, first_valid_row + 1))[1]
+            crop_bottom = max(dist_bounds.bottom, valid_bottom_y)
+            
+            fusion_bounds[date_str] = (
+                dist_bounds.left, crop_bottom, dist_bounds.right, dist_bounds.top, dist_crs
+            )
+            fusion_rows[date_str] = first_valid_row
+    
+    # Process S2 files
+    for refl_file in s2_prep.glob("*REFL.tif"):
+        date_str = refl_file.stem.split("_")[2]
+        if date_str in fusion_bounds:
+            output_file = processed / "s2" / f"{date_str}_0.geotiff"
+            if _crop_to_bounds(refl_file, fusion_bounds[date_str], output_file):
+                print(f"[PROCESS] Saved: {output_file}")
+    
+    # Process S3 files
+    for s3_file in s3_prep.glob("composite_*.tif"):
+        date_str = s3_file.stem.split("_")[1]
+        if date_str in fusion_bounds:
+            output_file = processed / "s3" / f"{date_str}_0.geotiff"
+            if _crop_to_bounds(s3_file, fusion_bounds[date_str], output_file):
+                print(f"[PROCESS] Saved: {output_file}")
+    
+    # Process fusion files (use row-based cropping)
+    for date_str, bounds in fusion_bounds.items():
+        fusion_file = fusion_prep / f"REFL_{date_str}.tif"
+        if fusion_file.exists():
+            output_file = processed / "fusion" / f"{date_str}_0.geotiff"
+            if _crop_to_bounds(fusion_file, bounds, output_file, row_based_height=fusion_rows[date_str] + 1):
+                print(f"[PROCESS] Saved: {output_file}")
+    
+    print("[PROCESS] Completed")

run.py

@@ -1,9 +1,10 @@
 from call_efast import run_efast, prepare_s2, prepare_s3
+from post_process import process_cropped
 from ndvi import (
     generate_ndvi_raw,
     create_ndvi_timeseries_raw,
-    generate_ndvi_prepared,
-    create_ndvi_timeseries_prepared,
+    generate_ndvi_post_process,
+    create_ndvi_timeseries_post_process,
 )
 from download_s2 import download_s2
 from download_s3 import download_s3
@@ -14,25 +15,27 @@ def run_pipeline(season, site_position, site_name):
     try:
         # print(f"Downloading data for {site_name}, {season}")
         # download_s2(season, site_position, site_name)
-        download_s3(season, site_position, site_name)
+        # download_s3(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)
+        # generate_ndvi_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)
+        # detect_clouds(season, 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"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)
+        # run_efast(season, site_position, site_name)
 
-        # print(f"Generating NDVI for prepared outputs: {site_name}, {season}")
-        generate_ndvi_prepared(season, site_position, site_name)
-        create_ndvi_timeseries_prepared(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)
+        create_ndvi_timeseries_post_process(season, site_position, site_name)
 
     except Exception as e:
         print(f"Error: {e}")

webapp/index.html

@@ -76,25 +76,34 @@
         const osmUrl = "https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png";
         const osmOpts = { attribution: "OpenStreetMap", opacity: 0.4 };
         const mapOpts = { zoomControl: false };
+        const sitePosition = [47.116171, 11.320308];
         const maps = {
-            s2: L.map("s2map", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
-            fusion: L.map("fusionmap", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
-            s3: L.map("s3map", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts))
+            s2: L.map("s2map", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
+            fusion: L.map("fusionmap", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
+            s3: L.map("s3map", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts))
         };
         const ndvimaps = {
-            s2: L.map("s2ndvimap", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
-            fusion: L.map("fusionndvimap", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
-            s3: L.map("s3ndvimap", mapOpts).setView([47.27, 11.39], 12).addLayer(L.tileLayer(osmUrl, osmOpts))
+            s2: L.map("s2ndvimap", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
+            fusion: L.map("fusionndvimap", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts)),
+            s3: L.map("s3ndvimap", mapOpts).setView(sitePosition, 12).addLayer(L.tileLayer(osmUrl, osmOpts))
         };
         const overlays = { s2: null, fusion: null, s3: null };
         const ndviOverlays = { s2: null, fusion: null, s3: null };
+        const markers = { s2: null, fusion: null, s3: null };
+        const ndviMarkers = { s2: null, fusion: null, s3: null };
         let timeseries = { s2: [], fusion: [], s3: [] };
         
+        // Add site marker to all maps
+        for (const source of ["s2", "fusion", "s3"]) {
+            markers[source] = L.marker(sitePosition, { icon: L.divIcon({ className: "site-marker", html: "<div style='width:5px;height:5px;background:red;border:1px solid white;border-radius:50%;box-shadow:0 0 1px rgba(0,0,0,0.5);'></div>", iconSize: [5, 5] }) }).addTo(maps[source]);
+            ndviMarkers[source] = L.marker(sitePosition, { icon: L.divIcon({ className: "site-marker", html: "<div style='width:5px;height:5px;background:red;border:1px solid white;border-radius:50%;box-shadow:0 0 1px rgba(0,0,0,0.5);'></div>", iconSize: [5, 5] }) }).addTo(ndvimaps[source]);
+        }
+        
         async function loadTimeseries() {
             const [s2, fusion, s3] = await Promise.all([
-                fetch("../data/innsbruck/2024/prepared/ndvi/s2/timeseries.json").then(r => r.json()),
-                fetch("../data/innsbruck/2024/prepared/ndvi/fusion/timeseries.json").then(r => r.json()).catch(() => []),
-                fetch("../data/innsbruck/2024/prepared/ndvi/s3/timeseries.json").then(r => r.json())
+                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())
             ]);
             timeseries = { s2, fusion, s3 };
             drawTimeseries();
@@ -190,25 +199,12 @@
                     const d = new Date(target);
                     d.setDate(d.getDate() + offset * dir);
                     const date = d.toISOString().split("T")[0].replace(/-/g, "");
-                    if (source === "fusion") {
-                        const filename = `REFL_${date}.tif`;
+                    // Processed files use DATE_0.geotiff format
+                    const filename = `${date}_0.geotiff`;
                     try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/fusion/${filename}`);
+                        const res = await fetch(`../data/innsbruck/2024/processed/${source}/${filename}`);
                         if (res.ok) return filename;
                     } catch {}
-                    } else if (source === "s2") {
-                        const filename = `S2A_MSIL2A_${date}_REFL.tif`;
-                        try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/${source}/${filename}`);
-                            if (res.ok) return filename;
-                        } catch {}
-                    } else if (source === "s3") {
-                        const filename = `composite_${date}.tif`;
-                        try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/${source}/${filename}`);
-                            if (res.ok) return filename;
-                        } catch {}
-                    }
                 }
             }
             return null;
@@ -220,7 +216,7 @@
         }
         
         async function loadGeotiff(source, filename) {
-            const path = `../data/innsbruck/2024/prepared/${source}/${filename}`;
+            const path = `../data/innsbruck/2024/processed/${source}/${filename}`;
             const tiff = await GeoTIFF.fromArrayBuffer(await (await fetch(path)).arrayBuffer());
             const image = await tiff.getImage();
             const rasters = await image.readRasters();
@@ -265,23 +261,14 @@
             overlays[source] = L.imageOverlay(canvas.toDataURL(), bounds, { opacity: 0.95 }).addTo(maps[source]);
             maps[source].fitBounds(bounds);
             
-            let dateStr;
-            if (source === "fusion") {
-                // REFL_20240101.tif -> extract 20240101
-                dateStr = filename.split("_")[1].replace(".tif", "");
-            } else if (source === "s2") {
-                // S2A_MSIL2A_20240101_REFL.tif -> extract 20240101
-                dateStr = filename.split("_")[2];
-            } else if (source === "s3") {
-                // composite_20240101.tif -> extract 20240101
-                dateStr = filename.split("_")[1].replace(".tif", "");
-            }
+            // Processed files use DATE_0.geotiff format: 20240101_0.geotiff -> extract 20240101
+            const dateStr = filename.split("_")[0];
             const date = `${dateStr.slice(0,4)}-${dateStr.slice(4,6)}-${dateStr.slice(6,8)}`;
             document.getElementById(`${source}rgbdate`).textContent = date;
         }
         
         async function loadNDVI(source, filename, dateStr) {
-            const tiff = await GeoTIFF.fromArrayBuffer(await (await fetch(`../data/innsbruck/2024/prepared/ndvi/${source}/${filename}`)).arrayBuffer());
+            const tiff = await GeoTIFF.fromArrayBuffer(await (await fetch(`../data/innsbruck/2024/processed/ndvi/${source}/${filename}`)).arrayBuffer());
             const image = await tiff.getImage();
             const data = Array.from((await image.readRasters())[0]);
             const width = image.getWidth();
@@ -317,15 +304,8 @@
             ndviOverlays[source] = L.imageOverlay(canvas.toDataURL(), bounds, { opacity: 0.95 }).addTo(ndvimaps[source]);
             ndvimaps[source].fitBounds(bounds);
             
-            let extractedDateStr;
-            if (source === "fusion") {
-                extractedDateStr = filename.split("_")[0];
-            } else if (source === "s2") {
-                // S2 NDVI files are now named YYYYMMDD_ndvi.geotiff
-                extractedDateStr = filename.split("_")[0];
-            } else if (source === "s3") {
-                extractedDateStr = filename.split("_")[1].split(".")[0];
-            }
+            // Processed NDVI files use DATE_0_ndvi.geotiff format: 20240101_0_ndvi.geotiff -> extract 20240101
+            const extractedDateStr = filename.split("_")[0];
             const date = `${extractedDateStr.slice(0,4)}-${extractedDateStr.slice(4,6)}-${extractedDateStr.slice(6,8)}`;
             document.getElementById(`${source}ndvidate`).textContent = date;
         }
@@ -344,26 +324,12 @@
                     const d = new Date(target);
                     d.setDate(d.getDate() + offset * dir);
                     const date = d.toISOString().split("T")[0].replace(/-/g, "");
-                    if (source === "s2") {
-                        // S2 NDVI files are now named YYYYMMDD_ndvi.geotiff
-                        const filename = `${date}_ndvi.geotiff`;
+                    // Processed NDVI files use DATE_0_ndvi.geotiff format
+                    const filename = `${date}_0_ndvi.geotiff`;
                     try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/ndvi/s2/${filename}`);
+                        const res = await fetch(`../data/innsbruck/2024/processed/ndvi/${source}/${filename}`);
                         if (res.ok) return filename;
                     } catch {}
-                    } else if (source === "fusion") {
-                        const filename = `${date}_ndvi.geotiff`;
-                        try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/ndvi/fusion/${filename}`);
-                            if (res.ok) return filename;
-                        } catch {}
-                    } else if (source === "s3") {
-                        const filename = `composite_${date}.geotiff`;
-                        try {
-                            const res = await fetch(`../data/innsbruck/2024/prepared/ndvi/s3/${filename}`);
-                            if (res.ok) return filename;
-                        } catch {}
-                    }
                 }
             }
             return null;