diff --git a/call_efast.py b/call_efast.py index 20543b5..acb7a6c 100644 --- a/call_efast.py +++ b/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, diff --git a/ndvi.py b/ndvi.py index 15ebb5e..2877a3e 100644 --- a/ndvi.py +++ b/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"]: - _process_ndvi_files( - input_dir, - output_dir, - f"PREPARED-{source.upper()}", - pattern=pattern, - output_namer=_get_output_name_prepared, - ) + 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"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") diff --git a/post_process.py b/post_process.py new file mode 100644 index 0000000..b207d7c --- /dev/null +++ b/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") diff --git a/run.py b/run.py index 35f789b..01c4b13 100644 --- a/run.py +++ b/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}") diff --git a/webapp/index.html b/webapp/index.html index 4e92d11..b375ee2 100644 --- a/webapp/index.html +++ b/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: "
", iconSize: [5, 5] }) }).addTo(maps[source]); + ndviMarkers[source] = L.marker(sitePosition, { icon: L.divIcon({ className: "site-marker", html: "
", 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`; - try { - const res = await fetch(`../data/innsbruck/2024/prepared/fusion/${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 {} - } + // Processed files use DATE_0.geotiff format + const filename = `${date}_0.geotiff`; + try { + const res = await fetch(`../data/innsbruck/2024/processed/${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`; - try { - const res = await fetch(`../data/innsbruck/2024/prepared/ndvi/s2/${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 {} - } + // Processed NDVI files use DATE_0_ndvi.geotiff format + const filename = `${date}_0_ndvi.geotiff`; + try { + const res = await fetch(`../data/innsbruck/2024/processed/ndvi/${source}/${filename}`); + if (res.ok) return filename; + } catch {} } } return null;