Added new endpoint to calculate change

BackendNew/main.py

@@ -6,6 +6,7 @@
 import os
 import json
 import tempfile
+import time
 from typing import Optional, Tuple
 
 from fastapi import FastAPI, Query, HTTPException
@@ -18,11 +19,25 @@ from rasterio.mask import mask
 from rasterio.enums import Resampling
 from rasterio.warp import calculate_default_transform, reproject
 from rasterio.io import MemoryFile
+from rasterio.features import geometry_window
+from rasterio.env import Env
+
 from pystac_client import Client
 import planetary_computer as pc
 import requests
 import numpy as np
-from shapely.geometry import mapping
+from shapely.geometry import mapping, shape
+
+# ------------------------------------------------------------------------------
+# Optional GDAL/Rasterio HTTP tuning for faster COG range requests
+# ------------------------------------------------------------------------------
+GDAL_ENV = Env(
+    GDAL_DISABLE_READDIR_ON_OPEN='EMPTY_DIR',
+    GDAL_HTTP_MULTIRANGE='YES',
+    CPL_VSIL_CURL_ALLOWED_EXTENSIONS='.tif',
+    VSI_CACHE='TRUE',
+    CPL_VSIL_CURL_NON_CACHED=''
+)
 
 # --- CRS / AOI helpers --------------------------------------------------------
 
@@ -98,6 +113,7 @@ def _build_true_color_from_bands(item, out_path: str, rgb_uint8: bool) -> str:
             raise RuntimeError(f"Item lacks required band {bn} to build RGB.")
         bands[bn] = pc.sign(a.href)
 
+    with GDAL_ENV:
         with rasterio.open(bands["B04"]) as rsrc, \
              rasterio.open(bands["B03"]) as gsrc, \
              rasterio.open(bands["B02"]) as bsrc:
@@ -268,7 +284,9 @@ def _fetch_least_cloudy_item(collection: str, geometry: dict, date_range: tuple,
 
 def _read_band_to_array_signed(asset_href: str) -> tuple:
     """Open a single-band raster and return (array, profile)."""
-    with rasterio.open(pc.sign(asset_href)) as src:
+    href = pc.sign(asset_href)
+    with GDAL_ENV:
+        with rasterio.open(href) as src:
             arr = src.read(1)
             profile = src.profile.copy()
     return arr, profile
@@ -326,7 +344,6 @@ def _write_float_geotiff(path: str, arr: np.ndarray, profile, nodata=-9999.0):
     with rasterio.open(path, "w", **clean) as dst:
         dst.write(arr_out, 1)
 
-
 def _resample_match(src_arr, src_profile, ref_profile, resampling=Resampling.bilinear):
     """Resample a single-band array to the reference profile's grid."""
     with MemoryFile() as mem_src:
@@ -375,6 +392,104 @@ def _summarize_ndvi(ndvi: np.ndarray, pixel_size_m: float) -> dict:
         "area_gt_0_6_m2": float(np.sum(v > 0.6) * area_per_pixel_m2),
     }
 
+# --- FAST in-memory NDVI metric helpers --------------------------------------
+
+def _read_band_over_aoi(asset_href: str, aoi_geom_wgs84: dict):
+    """
+    Read a single-band raster only over the AOI window, returning (array, profile).
+    Uses HTTP range requests against COGs and returns float32 with NaNs for nodata.
+    """
+    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)
+
+            transform = src.window_transform(win)
+            prof = {
+                "height": arr.shape[0],
+                "width": arr.shape[1],
+                "count": 1,
+                "dtype": "float32",
+                "crs": src.crs,
+                "transform": transform
+            }
+            return arr, prof
+
+def _ndvi_metric_from_arrays(nir_arr, nir_prof, red_arr, red_prof, scale_info, metric: str) -> float:
+    # scale to reflectance
+    nir = _scale_reflectance(nir_arr, scale_info)
+    red = _scale_reflectance(red_arr, scale_info)
+
+    # align grids if needed
+    same_grid = (
+        red_prof["transform"] == nir_prof["transform"] and
+        red_prof["crs"] == nir_prof["crs"] and
+        red_arr.shape == nir_arr.shape
+    )
+    if not same_grid:
+        nir = _resample_match(nir, nir_prof, red_prof, resampling=Resampling.bilinear)
+
+    ndvi = _compute_ndvi(nir, red)
+
+    valid = np.isfinite(ndvi)
+    if not np.any(valid):
+        return 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."""
+    cfg = collection_for_year(year)
+    date_range = year_to_range(year)
+    item = _fetch_least_cloudy_item(cfg["collection"], geometry_wgs84, date_range, max_cloud)
+    if item is None:
+        raise HTTPException(status_code=404, detail=f"No imagery for year={year}")
+
+    red_href = item.assets[cfg["red"]].href
+    nir_href = item.assets[cfg["nir"]].href
+
+    # Read only AOI window for each band
+    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)
+    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
+
 # --- API setup ----------------------------------------------------------------
 
 app = FastAPI(
@@ -384,7 +499,7 @@ An API that returns a satellite image file for a given year and computes NDVI &
 - Swagger UI: `/docs`
 - ReDoc: `/redoc`
 """,
-    version="1.3.0"
+    version="1.4.0"
 )
 
 # Defaults; override via env or query params
@@ -455,6 +570,7 @@ def download_imagery(
     # Optional: clip to AOI (must reproject AOI to raster CRS before masking)
     final_path = out_path
     if clip_to_aoi:
+        with GDAL_ENV:
             with rasterio.open(out_path) as src:
                 geoms = gpd.read_file(aoi_path).to_crs(src.crs)
                 shapes = [mapping(geom) for geom in geoms.geometry]
@@ -510,7 +626,7 @@ def download_imagery_post(params: DownloadParams):
         build_overviews=params.build_overviews,
     )
 
-# --- NDVI Endpoints -----------------------------------------------------------
+# --- NDVI Endpoints (original, file-based) -----------------------------------
 
 @app.get("/ndvi", summary="Compute NDVI for a given year; returns stats and writes a GeoTIFF")
 def ndvi_year(
@@ -555,6 +671,7 @@ def ndvi_year(
     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")
@@ -565,6 +682,7 @@ def ndvi_year(
 
     final_ndvi = tmp_ndvi
     if clip_to_aoi:
+        with GDAL_ENV:
             with rasterio.open(tmp_ndvi) as src:
                 geoms = gpd.read_file(aoi_path).to_crs(src.crs)
                 shapes = [mapping(geom) for geom in geoms.geometry]
@@ -578,6 +696,7 @@ def ndvi_year(
 
     # 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:
                 dst_transform, width, height = calculate_default_transform(
                     src.crs, src.crs, src.width, src.height, *src.bounds, resolution=res_for_compare_m
@@ -643,6 +762,7 @@ def ndvi_change(
     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):
@@ -690,5 +810,42 @@ def ndvi_change(
         "change_stats": change_stats
     }
 
+# --- 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)")
+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)")
+):
+    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:
+            aoi_path = os.path.join(tmp, "aoi_buffered.geojson")
+            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)
+
+    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)
+    }
+
 # Run:
 # uvicorn main:app --host 0.0.0.0 --port 8000