Updated calculation endpoint

BackendNew/main.py

@@ -250,7 +250,6 @@ def collection_for_year(year: int):
     Sentinel-2 L2A for 2016+, Landsat 8/9 L2 for earlier.
     """
     if year >= 2016:
-        # Sentinel-2 L2A on Planetary Computer
         return {
             "collection": "sentinel-2-l2a",
             "red": "B04",
@@ -259,7 +258,6 @@ def collection_for_year(year: int):
             "native_res_m": 10.0,
         }
     else:
-        # Landsat 8 Collection 2 Level-2 Surface Reflectance (pre-2016)
         return {
             "collection": "landsat-8-c2-l2",
             "red": "SR_B4",
@@ -296,10 +294,8 @@ def _scale_reflectance(arr: np.ndarray, scale_info) -> np.ndarray:
     kind, params = scale_info
     arr = arr.astype("float32")
     if kind == "sentinel":
-        # Planetary Computer S2 L2A are typically scaled 0..10000
         return np.clip(arr / 10000.0, 0.0, 1.0)
     elif kind == "landsat_sr":
-        # Landsat C2 L2 SR scale/offset
         mult, off = params
         return np.clip(arr * mult + off, 0.0, 1.0)
     else:
@@ -315,15 +311,11 @@ def _compute_ndvi(nir: np.ndarray, red: np.ndarray, nodata_mask: np.ndarray = No
 def _write_float_geotiff(path: str, arr: np.ndarray, profile, nodata=-9999.0):
     """
     Write a single-band float32 GeoTIFF.
-    - Use a finite nodata (e.g. -9999.0) for best GDAL compatibility on Windows.
-    - Sanitize profile to avoid conflicts from copied profiles.
     """
     arr_out = arr.astype("float32", copy=False)
-    # Replace non-finite values with nodata
     if not np.isfinite(arr_out).all():
         arr_out = np.where(np.isfinite(arr_out), arr_out, nodata).astype("float32")
 
-    # Start from a minimal, clean profile
     clean = {
         "driver": "GTiff",
         "height": profile["height"],
@@ -338,9 +330,7 @@ def _write_float_geotiff(path: str, arr: np.ndarray, profile, nodata=-9999.0):
         "nodata": float(nodata),
     }
 
-    # Make sure directory exists
     os.makedirs(os.path.dirname(path), exist_ok=True)
-
     with rasterio.open(path, "w", **clean) as dst:
         dst.write(arr_out, 1)
 
@@ -402,22 +392,18 @@ def _read_band_over_aoi(asset_href: str, aoi_geom_wgs84: dict):
     href = pc.sign(asset_href)
     with GDAL_ENV:
         with rasterio.open(href) as src:
-            # Reproject AOI to raster CRS and compute minimal read window
             aoi_geom = shape(aoi_geom_wgs84)
             aoi = gpd.GeoSeries([aoi_geom], crs="EPSG:4326").to_crs(src.crs)
             geom = aoi.iloc[0]
 
-            # Compute the read window; if AOI is completely outside, bail early
             try:
                 win = geometry_window(src, [mapping(geom)], pad_x=0, pad_y=0, north_up=True, pixel_precision=3)
             except ValueError:
-                # No overlap
                 return np.full((0, 0), np.nan, dtype="float32"), {
                     "height": 0, "width": 0, "count": 1, "dtype": "float32",
                     "crs": src.crs, "transform": src.transform
                 }
 
-            # Windowed read as float32 masked array; fill mask with NaN
             m = src.read(1, window=win, boundless=False, masked=True, out_dtype="float32")
             arr = m.filled(np.nan)
 
@@ -432,7 +418,10 @@ def _read_band_over_aoi(asset_href: str, aoi_geom_wgs84: dict):
             }
             return arr, prof
 
-def _ndvi_metric_from_arrays(nir_arr, nir_prof, red_arr, red_prof, scale_info, metric: str) -> float:
+def _ndvi_metrics_from_arrays(nir_arr, nir_prof, red_arr, red_prof, scale_info) -> dict:
+    """
+    Compute NDVI and return all summary metrics we care about: mean, median, p05, p95.
+    """
     # scale to reflectance
     nir = _scale_reflectance(nir_arr, scale_info)
     red = _scale_reflectance(red_arr, scale_info)
@@ -450,24 +439,17 @@ def _ndvi_metric_from_arrays(nir_arr, nir_prof, red_arr, red_prof, scale_info, m
 
     valid = np.isfinite(ndvi)
     if not np.any(valid):
-        return float("nan")
+        return {"mean": float("nan"), "median": float("nan"), "p05": float("nan"), "p95": float("nan")}
     v = ndvi[valid]
-    metric = (metric or "mean").lower()
-    if metric == "mean":
-        return float(np.nanmean(v))
-    if metric == "median":
-        return float(np.nanmedian(v))
-    if metric.startswith("p"):
-        try:
-            q = float(metric[1:])
-        except Exception:
-            q = 50.0
-        return float(np.nanpercentile(v, q))
-    # default to mean
-    return float(np.nanmean(v))
-
-def _ndvi_metric_for_year(year: int, geometry_wgs84: dict, max_cloud: int, metric: str) -> Tuple[float, dict]:
-    """Fast, in-memory NDVI metric for a year over AOI."""
+    return {
+        "mean": float(np.nanmean(v)),
+        "median": float(np.nanmedian(v)),
+        "p05": float(np.nanpercentile(v, 5)),
+        "p95": float(np.nanpercentile(v, 95)),
+    }
+
+def _ndvi_metrics_for_year(year: int, geometry_wgs84: dict, max_cloud: int) -> Tuple[dict, dict]:
+    """Fast, in-memory NDVI metrics (mean/median/p05/p95) for a year over AOI."""
     cfg = collection_for_year(year)
     date_range = year_to_range(year)
     item = _fetch_least_cloudy_item(cfg["collection"], geometry_wgs84, date_range, max_cloud)
@@ -481,14 +463,14 @@ def _ndvi_metric_for_year(year: int, geometry_wgs84: dict, max_cloud: int, metri
     red_arr, red_prof = _read_band_over_aoi(red_href, geometry_wgs84)
     nir_arr, nir_prof = _read_band_over_aoi(nir_href, geometry_wgs84)
 
-    value = _ndvi_metric_from_arrays(nir_arr, nir_prof, red_arr, red_prof, cfg["scale"], metric)
+    metrics = _ndvi_metrics_from_arrays(nir_arr, nir_prof, red_arr, red_prof, cfg["scale"])
     meta = {
         "collection": cfg["collection"],
         "native_res_m": cfg["native_res_m"],
         "cloud_cover": item.properties.get("eo:cloud_cover"),
         "item_id": item.id
     }
-    return value, meta
+    return metrics, meta
 
 # --- API setup ----------------------------------------------------------------
 
@@ -499,7 +481,7 @@ An API that returns a satellite image file for a given year and computes NDVI &
 - Swagger UI: `/docs`
 - ReDoc: `/redoc`
 """,
-    version="1.4.0"
+    version="1.5.0"
 )
 
 # Defaults; override via env or query params
@@ -552,11 +534,9 @@ def download_imagery(
     if not os.path.exists(AOI_PATH):
         raise HTTPException(status_code=500, detail=f"AOI file not found at {AOI_PATH}")
 
-    # Prepare temp workspace
     workdir = tempfile.mkdtemp(prefix=f"imagery_{year}_")
     aoi_path = AOI_PATH
 
-    # Optional buffering (result saved as WGS84 GeoJSON for STAC search)
     if buffer_meters and buffer_meters > 0:
         aoi_path = os.path.join(workdir, "aoi_buffered.geojson")
         buffer_aoi(AOI_PATH, buffer_meters, aoi_path)
@@ -567,7 +547,6 @@ def download_imagery(
     if not out_path:
         raise HTTPException(status_code=404, detail=f"No imagery found for year={year}, range={date_range}, max_cloud={max_cloud}")
 
-    # Optional: clip to AOI (must reproject AOI to raster CRS before masking)
     final_path = out_path
     if clip_to_aoi:
         with GDAL_ENV:
@@ -592,17 +571,14 @@ def download_imagery(
                 with rasterio.open(clipped_tif, "w", **out_meta) as dest:
                     dest.write(out_image)
 
-        # Reproject to a clean target resolution (e.g., 10 m) for crisp pixels
         reproj_tif = os.path.join(workdir, f"imagery_{year}_clipped_{int(target_res_m)}m.tif")
         reproject_to_resolution(clipped_tif, reproj_tif, target_res_m=target_res_m)
 
-        # Add internal overviews for sharp display at various zoom levels
         if build_overviews:
             add_overviews_inplace(reproj_tif)
 
         final_path = reproj_tif
 
-    # Serve the file
     return FileResponse(
         final_path,
         media_type="image/tiff",
@@ -653,30 +629,25 @@ def ndvi_year(
     if item is None:
         raise HTTPException(status_code=404, detail="No imagery found")
 
-    # Read bands
     red_href = item.assets[cfg["red"]].href
     nir_href = item.assets[cfg["nir"]].href
     red, red_prof = _read_band_to_array_signed(red_href)
     nir, nir_prof = _read_band_to_array_signed(nir_href)
 
-    # Scale reflectance to 0..1
     red = _scale_reflectance(red, cfg["scale"])
     nir = _scale_reflectance(nir, cfg["scale"])
 
-    # Ensure bands are on identical grid (usually true, but safeguard)
     same_grid = (red_prof["transform"] == nir_prof["transform"]) and (red_prof["crs"] == nir_prof["crs"]) and (red.shape == nir.shape)
     if not same_grid:
         nir = _resample_match(nir, nir_prof, red_prof, resampling=Resampling.bilinear)
 
     ndvi = _compute_ndvi(nir, red)
 
-    # Base profile from red band
     with GDAL_ENV:
         with rasterio.open(pc.sign(red_href)) as ref:
             base_profile = ref.profile.copy()
     base_profile.update(count=1, dtype="float32")
 
-    # Write temp NDVI before clipping
     tmp_ndvi = os.path.join(workdir, f"ndvi_{year}_raw.tif")
     _write_float_geotiff(tmp_ndvi, ndvi, base_profile)
 
@@ -694,7 +665,6 @@ def ndvi_year(
                     dst.write(out_image)
                 final_ndvi = clipped
 
-    # Optionally resample to a standard comparison grid (e.g., 30 m)
     if res_for_compare_m and res_for_compare_m > 0:
         with GDAL_ENV:
             with rasterio.open(final_ndvi) as src:
@@ -714,7 +684,6 @@ def ndvi_year(
                     )
                 final_ndvi = repro
 
-    # Compute stats
     with rasterio.open(final_ndvi) as ds:
         arr = ds.read(1)
         resx = abs(ds.transform.a)
@@ -742,7 +711,6 @@ def ndvi_change(
     if end_year < start_year:
         raise HTTPException(status_code=422, detail="end_year must be >= start_year")
 
-    # Reuse the /ndvi logic twice, then difference
     start = ndvi_year(
         year=start_year,
         max_cloud=max_cloud,
@@ -761,10 +729,8 @@ def ndvi_change(
     ndvi_start_path = os.path.join(start["workdir"], start["ndvi_tif"])
     ndvi_end_path   = os.path.join(end["workdir"], end["ndvi_tif"])
 
-    # Align and difference
     with GDAL_ENV:
         with rasterio.open(ndvi_end_path) as end_ds, rasterio.open(ndvi_start_path) as start_ds:
-            # If shapes differ (shouldn't if same res), resample start to end grid
             if (start_ds.transform != end_ds.transform) or (start_ds.width != end_ds.width) or (start_ds.height != end_ds.height):
                 arr_s = _resample_match(start_ds.read(1), start_ds.profile, end_ds.profile, resampling=Resampling.bilinear)
                 arr_e = end_ds.read(1)
@@ -776,7 +742,6 @@ def ndvi_change(
 
             change = (arr_e - arr_s).astype("float32")
 
-            # Save change raster
             workdir = tempfile.mkdtemp(prefix=f"ndvi_change_{start_year}_{end_year}_")
             change_tif = os.path.join(workdir, f"ndvi_change_{start_year}_{end_year}.tif")
             _write_float_geotiff(change_tif, change, profile)
@@ -785,7 +750,6 @@ def ndvi_change(
             stats_start = _summarize_ndvi(arr_s, resx)
             stats_end = _summarize_ndvi(arr_e, resx)
 
-            # Change stats: summarize positive/negative change
             valid = np.isfinite(change)
             ch = change[valid]
             change_stats = {
@@ -812,21 +776,23 @@ def ndvi_change(
 
 # --- FAST NDVI change VALUE endpoint (in-memory, no file writes) --------------
 
-@app.get("/ndvi-change/value", summary="Fast NDVI change value (in-memory, no file writes)")
+@app.get("/ndvi-change/value", summary="Fast NDVI change values (mean/median/p05/p95)")
 def ndvi_change_value(
     start_year: int = Query(..., ge=2013),
     end_year:   int = Query(..., ge=2013),
-    metric: str = Query("mean", description="NDVI summary to compare: mean | median | p05 | p95"),
     max_cloud: int = Query(MAX_CLOUD_DEFAULT, ge=0, le=100),
     buffer_meters: float = Query(BUFFER_METERS_DEFAULT, ge=0, description="Optional AOI buffer (m)")
 ):
+    """
+    Returns summary NDVI metrics for start & end years (mean, median, p05, p95)
+    and their deltas (end - start). Always computes all metrics.
+    """
     if end_year < start_year:
         raise HTTPException(status_code=422, detail="end_year must be >= start_year")
     if not os.path.exists(AOI_PATH):
         raise HTTPException(status_code=500, detail=f"AOI file not found at {AOI_PATH}")
 
     t0 = time.time()
-    # prepare (buffered) AOI once
     with tempfile.TemporaryDirectory(prefix="aoi_") as tmp:
         aoi_path = AOI_PATH
         if buffer_meters and buffer_meters > 0:
@@ -834,18 +800,28 @@ def ndvi_change_value(
             buffer_aoi(AOI_PATH, buffer_meters, aoi_path)
         geom_wgs84 = read_aoi_geom_wgs84(aoi_path)
 
-        v_start, meta_start = _ndvi_metric_for_year(start_year, geom_wgs84, max_cloud, metric)
-        v_end,   meta_end   = _ndvi_metric_for_year(end_year,   geom_wgs84, max_cloud, metric)
+        start_metrics, meta_start = _ndvi_metrics_for_year(start_year, geom_wgs84, max_cloud)
+        end_metrics,   meta_end   = _ndvi_metrics_for_year(end_year,   geom_wgs84, max_cloud)
 
     elapsed_s = time.time() - t0
-    change_value = None if (np.isnan(v_start) or np.isnan(v_end)) else float(v_end - v_start)
-    return {
-        "metric": metric,
-        "start": {"year": start_year, "value": v_start, **meta_start},
-        "end":   {"year": end_year,   "value": v_end,   **meta_end},
-        "change_value": change_value,
-        "elapsed_seconds": round(elapsed_s, 3)
+
+    def delta(key):
+        sv = start_metrics.get(key)
+        ev = end_metrics.get(key)
+        if sv is None or ev is None or np.isnan(sv) or np.isnan(ev):
+            return None
+        return float(ev - sv)
+
+    change = {
+        "mean": delta("mean"),
+        "median": delta("median"),
+        "p05": delta("p05"),
+        "p95": delta("p95"),
     }
 
-# Run:
-# uvicorn main:app --host 0.0.0.0 --port 8000
+    return {
+        "start": {"year": start_year, "metrics": start_metrics, **meta_start},
+        "end":   {"year": end_year,   "metrics": end_metrics,   **meta_end},
+        "change": change,
+        "elapsed_seconds": round(elapsed_s, 3)
+    }
\ No newline at end of file