diff --git a/call_efast.py b/call_efast.py index a9aec22..54ade92 100644 --- a/call_efast.py +++ b/call_efast.py @@ -9,18 +9,21 @@ from rasterio.warp import Resampling from rasterio.vrt import WarpedVRT from rasterio import shutil as rio_shutil + def _import_efast(): """Lazy import of efast to avoid import errors when not using efast functions.""" try: import efast from efast.s2_processing import distance_to_clouds from efast.s3_processing import reproject_and_crop_s3 + return efast, distance_to_clouds, reproject_and_crop_s3 except ImportError: raise ImportError( "efast package not found. Install with: pip install git+https://github.com/DHI-GRAS/efast.git" ) + RESOLUTION_RATIO = 21 @@ -33,11 +36,22 @@ def _load_clouds(clouds_file): return clouds -def _reproject_raster_to_target(src_path, dst_path, target_bounds, target_crs, width, height, resampling=Resampling.cubic): +def _reproject_raster_to_target( + src_path, + dst_path, + target_bounds, + target_crs, + width, + height, + resampling=Resampling.cubic, +): dst_transform = rasterio.transform.from_bounds( - target_bounds.left, target_bounds.bottom, - target_bounds.right, target_bounds.top, - width, height + target_bounds.left, + target_bounds.bottom, + target_bounds.right, + target_bounds.top, + width, + height, ) with rasterio.open(src_path) as src: vrt_options = { @@ -88,7 +102,9 @@ def prepare_s2(season, site_position, site_name, date_range=None): with rasterio.open(temp_normalized, "w", **profile) as dst: dst.write(data) - _reproject_raster_to_target(temp_normalized, refl_dst, target_bounds, target_crs, s2_width, s2_height) + _reproject_raster_to_target( + temp_normalized, refl_dst, target_bounds, target_crs, s2_width, s2_height + ) temp_normalized.unlink() _, distance_to_clouds, _ = _import_efast() @@ -149,6 +165,8 @@ def run_efast(season, site_position, site_name, date_range=None): fusion_output_dir.mkdir(parents=True, exist_ok=True) print(f"[EFAST] Starting fusion: {site_name} ({lat:.6f}, {lon:.6f}), {season}") + efast, _, _ = _import_efast() + start_str, end_str = datetime_range.split("/") start_date = datetime.strptime(start_str, "%Y-%m-%d") end_date = datetime.strptime(end_str, "%Y-%m-%d") @@ -157,18 +175,25 @@ def run_efast(season, site_position, site_name, date_range=None): while current_date <= end_date: date_str = current_date.strftime("%Y%m%d") output_file = fusion_output_dir / f"REFL_{date_str}.tif" - if output_file.exists(): - print(f"[EFAST] Skipping {date_str} (exists)") - else: - 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, - ) - print(f"[EFAST] Saved: {output_file}" if output_file.exists() else f"[EFAST] No output for {date_str} (insufficient nearby data)") - except Exception as e: - print(f"[EFAST] Error processing {date_str}: {e}") + 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, + ) + print( + f"[EFAST] Saved: {output_file}" + if output_file.exists() + else f"[EFAST] No output for {date_str} (insufficient nearby data)" + ) + except Exception as e: + print(f"[EFAST] Error processing {date_str}: {e}") current_date += timedelta(days=1) print("[EFAST] Completed") diff --git a/ndvi.py b/ndvi.py index 84fd348..bd3b588 100644 --- a/ndvi.py +++ b/ndvi.py @@ -38,14 +38,26 @@ def _get_ndvi_value(ndvi_file, site_position): 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]) + + # Check if point is within bounds + if not ( + src.bounds.left <= x[0] <= src.bounds.right + and src.bounds.bottom <= y[0] <= src.bounds.top + ): + return None # Point is outside raster bounds + samples = list(src.sample([(x[0], y[0])])) if samples: value = float(samples[0][0]) - if value != 0 and not np.isnan(value): - return value - # Return the raw value even if 0 or NaN for diagnostic purposes + # Check if it's actually nodata (using raster's nodata value) + if src.nodata is not None and value == src.nodata: + return None # This is nodata, not a valid 0 value + if np.isnan(value): + return None # NaN is invalid + # 0 is a valid NDVI value (no vegetation), so return it return value - except Exception: + except Exception as e: + print(f"Error sampling {ndvi_file.name}: {e}") pass return None @@ -56,21 +68,41 @@ def _create_timeseries_for_dir(output_dir, site_position, source_name): for ndvi_file in sorted(output_dir.glob("*.geotiff")): filename = ndvi_file.name - date_str = filename.split("_")[0] - try: - date = datetime.strptime(date_str, "%Y%m%d").isoformat() - except ValueError: + # Extract date from filename + # Format examples: + # - YYYYMMDD_ndvi.geotiff -> date is at [0] + # - YYYYMMDD_0.geotiff -> date is at [0] + # - composite_YYYYMMDD.geotiff -> date is at [1] + parts = filename.replace(".geotiff", "").split("_") + date_str = None + + # Try to find a date pattern (8 digits) + for part in parts: + if len(part) == 8 and part.isdigit(): + date_str = part + break + + if date_str: + try: + date = datetime.strptime(date_str, "%Y%m%d").isoformat() + except ValueError: + date = date_str + else: + # Fallback: use first part (for old MSIL2A_ndvi.geotiff files) + date_str = parts[0] date = date_str + print( + f"[NDVI-{source_name}] Warning: Could not extract date from {filename}, using '{date_str}'" + ) ndvi_value = _get_ndvi_value(ndvi_file, site_position) if ndvi_value is None: - print(f"[NDVI-{source_name}] Warning: Could not sample {filename}") - elif ndvi_value == 0: - print(f"[NDVI-{source_name}] Warning: Could not sample {filename} (NoData)") - ndvi_value = None # Set to None for timeseries - elif np.isnan(ndvi_value): - print(f"[NDVI-{source_name}] Warning: Could not sample {filename} (NaN)") - ndvi_value = None # Set to None for timeseries + print( + f"[NDVI-{source_name}] Warning: Could not sample {filename} (outside bounds or nodata)" + ) + # Note: 0 is a valid NDVI value (no vegetation), so we keep it + # The _get_ndvi_value function now properly distinguishes between + # valid 0 values and nodata values timeseries.append({"date": date, "filename": filename, "ndvi": ndvi_value}) @@ -98,20 +130,20 @@ def _process_ndvi_files( try: with rasterio.open(geotiff_file) as src: if src.count < 4: - print(f"[NDVI-{source_name}] Skipping {geotiff_file.name} (only {src.count} band(s), need 4+)") + print( + f"[NDVI-{source_name}] Skipping {geotiff_file.name} (only {src.count} band(s), need 4+)" + ) continue except Exception as e: - print(f"[NDVI-{source_name}] Skipping {geotiff_file.name} (error reading: {e})") + print( + f"[NDVI-{source_name}] Skipping {geotiff_file.name} (error reading: {e})" + ) continue output_file = output_dir / ( output_namer(geotiff_file) if output_namer else geotiff_file.name ) - if output_file.exists(): - print(f"[NDVI-{source_name}] Skipping {geotiff_file.name} (exists)") - continue - _calculate_and_write_ndvi(geotiff_file, output_file) print(f"[NDVI-{source_name}] Saved: {output_file}") @@ -132,7 +164,18 @@ def create_ndvi_timeseries_raw(season, site_position, site_name): def _get_output_name_prepared(geotiff_file): if geotiff_file.suffix == ".tif": if "REFL" in geotiff_file.stem: - date_str = geotiff_file.stem.split("_")[1] + # For S2: S2A_MSIL2A_20240101_REFL -> date is at index [2] + # For S3: composite_20240101.tif -> date is at index [1] after removing .tif + parts = geotiff_file.stem.split("_") + if len(parts) >= 3 and parts[0].startswith("S2"): + # S2 format: S2A_MSIL2A_YYYYMMDD_REFL + date_str = parts[2] + elif len(parts) >= 2 and parts[0] == "composite": + # S3 format: composite_YYYYMMDD + date_str = parts[1] + else: + # Fallback: try index [1] for other formats + date_str = parts[1] if len(parts) > 1 else parts[0] return f"{date_str}_ndvi.geotiff" return geotiff_file.name.replace(".tif", ".geotiff") return geotiff_file.name diff --git a/run.py b/run.py index 929c072..16149c2 100644 --- a/run.py +++ b/run.py @@ -24,11 +24,11 @@ def run_pipeline(season, site_position, 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) + 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) diff --git a/webapp/cloudy.html b/webapp/cloudy.html index 00dae3a..b45e660 100644 --- a/webapp/cloudy.html +++ b/webapp/cloudy.html @@ -18,6 +18,7 @@ .timeseries-label { font-size: 12px; margin-bottom: 5px; color: #666; } .timeseries { width: 100%; height: 120px; border: 1px solid #ccc; margin-bottom: 10px; } .map-label { font-size: 12px; margin-bottom: 5px; color: #666; } + .map-date { font-size: 11px; margin-top: 5px; color: #999; } .map { height: 500px; border: 1px solid #ccc; } .leaflet-image-layer { image-rendering: pixelated; } .leaflet-control-attribution { display: none; } @@ -28,6 +29,10 @@