diff --git a/ndvi.py b/ndvi.py index f12b62f..20fe2bc 100644 --- a/ndvi.py +++ b/ndvi.py @@ -94,3 +94,88 @@ def create_ndvi_timeseries(year, site_position, site_name): print( f"[NDVI-{source.upper()}] Saved: {timeseries_file} ({len(timeseries)} entries)" ) + + +def generate_ndvi_fusion(year, site_position, site_name): + input_dir = Path(f"data/{site_name}/{year}/efast/fusion/") + output_dir = Path(f"data/{site_name}/{year}/ndvi/fusion/") + output_dir.mkdir(parents=True, exist_ok=True) + + print(f"[NDVI-FUSION] Processing {input_dir}...") + + geotiff_files = sorted(input_dir.glob("REFL_*.tif")) + if not geotiff_files: + print(f"[NDVI-FUSION] No files found") + return + + for geotiff_file in geotiff_files: + date_str = geotiff_file.stem.split("_")[1] + output_file = output_dir / f"{date_str}_ndvi.geotiff" + + if output_file.exists(): + print(f"[NDVI-FUSION] Skipping {geotiff_file.name} (exists)") + continue + + with rasterio.open(geotiff_file) as src: + red = src.read(3).astype(np.float32) + nir = src.read(4).astype(np.float32) + + mask = (red > 0) & (nir > 0) + ndvi = np.zeros_like(red, dtype=np.float32) + ndvi[mask] = (nir[mask] - red[mask]) / (nir[mask] + red[mask]) + + profile = src.profile.copy() + profile.update( + { + "count": 1, + "dtype": "float32", + "nodata": 0, + "compress": "lzw", + } + ) + + with rasterio.open(output_file, "w", **profile) as dst: + dst.write(ndvi, 1) + dst.set_band_description(1, "NDVI") + + print(f"[NDVI-FUSION] Saved: {output_file}") + + print(f"[NDVI-FUSION] Completed") + + +def create_ndvi_timeseries_fusion(year, site_position, site_name): + output_dir = Path(f"data/{site_name}/{year}/ndvi/fusion/") + + print(f"[NDVI-FUSION] Creating timeseries.json...") + timeseries = [] + + ndvi_files = sorted(output_dir.glob("*.geotiff")) + for ndvi_file in ndvi_files: + filename = ndvi_file.name + date_str = filename.split("_")[0] + try: + date = datetime.strptime(date_str, "%Y%m%d").isoformat() + except ValueError: + date = date_str + + ndvi_value = None + try: + with rasterio.open(ndvi_file) as src: + lon, lat = site_position[1], site_position[0] + x, y = transform_coords("EPSG:4326", src.crs, [lon], [lat]) + samples = list(src.sample([(x[0], y[0])])) + if samples and len(samples) > 0: + value = float(samples[0][0]) + if value != 0 and not np.isnan(value): + ndvi_value = value + except Exception as e: + print(f"[NDVI-FUSION] Warning: Could not sample {filename}: {e}") + + timeseries.append({"date": date, "filename": filename, "ndvi": ndvi_value}) + + timeseries.sort(key=lambda x: x["date"]) + timeseries_file = output_dir / "timeseries.json" + with open(timeseries_file, "w") as f: + json.dump(timeseries, f, indent=2) + + print(f"[NDVI-FUSION] Saved: {timeseries_file} ({len(timeseries)} entries)") diff --git a/run.py b/run.py index ae58b5e..653b053 100644 --- a/run.py +++ b/run.py @@ -1,4 +1,5 @@ from efast import run_efast, prepare_s2, prepare_s3 +from ndvi import generate_ndvi_fusion, create_ndvi_timeseries_fusion year = 2024 site_position = (47.116171, 11.320308) @@ -18,11 +19,16 @@ site_name = "innsbruck" # detect_clouds(year, site_name) # print("Cloud detection completed") -print(f"Preparing data for EFAST fusion for {site_name}, {year}") -prepare_s2(year, site_position, site_name) -prepare_s3(year, site_position, site_name) -print("Data preparation completed") +# print(f"Preparing data for EFAST fusion for {site_name}, {year}") +# prepare_s2(year, site_position, site_name) +# prepare_s3(year, site_position, site_name) +# print("Data preparation completed") -print(f"Running EFAST fusion for {site_name}, {year}") -run_efast(year, site_position, site_name) -print("EFAST fusion completed") +# print(f"Running EFAST fusion for {site_name}, {year}") +# run_efast(year, site_position, site_name) +# print("EFAST fusion completed") + +print(f"Generating NDVI for fusion outputs: {site_name}, {year}") +generate_ndvi_fusion(year, site_position, site_name) +create_ndvi_timeseries_fusion(year, site_position, site_name) +print("Fusion NDVI generation completed")