Merge remote-tracking branch 'origin/backend' into frontend

BackendNew/main.py

+# main.py
+# Requirements:
+#   pip install fastapi uvicorn geopandas rasterio pystac-client planetary-computer requests pydantic shapely numpy
+# (plus system deps needed by rasterio/GEOS/GDAL depending on your env)
+
+import os
+import json
+import tempfile
+from typing import Optional, Tuple
+
+from fastapi import FastAPI, Query, HTTPException
+from fastapi.responses import FileResponse
+from pydantic import BaseModel, Field
+
+import geopandas as gpd
+import rasterio
+from rasterio.mask import mask
+from rasterio.enums import Resampling
+from rasterio.warp import calculate_default_transform, reproject
+from rasterio.io import MemoryFile
+from pystac_client import Client
+import planetary_computer as pc
+import requests
+import numpy as np
+from shapely.geometry import mapping
+
+# --- CRS / AOI helpers --------------------------------------------------------
+
+def read_aoi_geom_wgs84(aoi_path: str) -> dict:
+    """
+    Load the AOI file and return a single GeoJSON geometry in EPSG:4326.
+    If there are multiple features, union them.
+    """
+    gdf = gpd.read_file(aoi_path)
+    if gdf.empty:
+        raise ValueError(f"AOI at {aoi_path} contains no features.")
+    gdf = gdf.to_crs(epsg=4326)
+    geom = gdf.unary_union
+    return mapping(geom)  # GeoJSON geometry dict
+
+def buffer_aoi(aoi_path: str, buffer_meters: float, output_path: str) -> str:
+    """
+    Buffer the AOI by N meters (using EPSG:3857 for meter-based buffering),
+    and save the result as a GeoJSON in EPSG:4326 for STAC search.
+    """
+    aoi = gpd.read_file(aoi_path)
+    if aoi.empty:
+        raise ValueError(f"AOI at {aoi_path} contains no features.")
+    aoi_merc = aoi.to_crs(epsg=3857)  # meter units
+    aoi_buffered = aoi_merc.buffer(buffer_meters)
+    aoi_buffered_gdf = gpd.GeoDataFrame(geometry=aoi_buffered, crs="EPSG:3857")
+    aoi_buffered_gdf = aoi_buffered_gdf.to_crs(epsg=4326)  # lon/lat for STAC
+    aoi_buffered_gdf.to_file(output_path, driver="GeoJSON")
+    return output_path
+
+# --- High-quality asset selection / composition -------------------------------
+
+def _pick_fullres_asset_or_none(item):
+    """
+    Return a full-resolution 'visual' GeoTIFF/COG asset if present.
+    Avoid thumbnails/overviews (PNG/JPEG renders).
+    """
+    a = item.assets.get("visual")
+    if not a:
+        return None
+    mt = (a.media_type or "").lower()
+    roles = [r.lower() for r in (a.roles or [])]
+    if any(r in roles for r in ["thumbnail", "overview", "rendered_preview"]):
+        return None
+    if ("tiff" in mt) or ("geotiff" in mt) or ("cog" in mt):
+        return a
+    return None
+
+def _scale_to_uint8(arr: np.ndarray) -> np.ndarray:
+    """
+    Simple 2–98 percentile contrast stretch to uint8.
+    Ignores zeros when computing percentiles (common for nodata).
+    """
+    mask_pos = arr > 0
+    if not np.any(mask_pos):
+        return np.zeros_like(arr, dtype="uint8")
+    p2, p98 = np.percentile(arr[mask_pos], (2, 98))
+    if p98 <= p2:
+        return np.clip(arr, 0, 255).astype("uint8")
+    scaled = (arr.astype("float32") - p2) * (255.0 / (p98 - p2))
+    return np.clip(scaled, 0, 255).astype("uint8")
+
+def _build_true_color_from_bands(item, out_path: str, rgb_uint8: bool) -> str:
+    """
+    Build a 10 m true-color GeoTIFF from Sentinel-2 L2A bands (B04,B03,B02).
+    - Keeps native dtype (usually uint16) unless rgb_uint8=True (applies 2–98 stretch).
+    - Writes a tiled, compressed GeoTIFF.
+    """
+    bands = {}
+    for bn in ["B04", "B03", "B02"]:
+        a = item.assets.get(bn)
+        if not a:
+            raise RuntimeError(f"Item lacks required band {bn} to build RGB.")
+        bands[bn] = pc.sign(a.href)
+
+    with rasterio.open(bands["B04"]) as rsrc, \
+         rasterio.open(bands["B03"]) as gsrc, \
+         rasterio.open(bands["B02"]) as bsrc:
+
+        red = rsrc.read(1)
+        grn = gsrc.read(1)
+        blu = bsrc.read(1)
+
+        profile = rsrc.profile.copy()
+        profile.update(
+            driver="GTiff",
+            count=3,
+            tiled=True,
+            compress="lzw",
+            predictor=2,
+            photometric="RGB",
+        )
+
+        if rgb_uint8:
+            red = _scale_to_uint8(red)
+            grn = _scale_to_uint8(grn)
+            blu = _scale_to_uint8(blu)
+            profile.update(dtype="uint8")
+
+        with rasterio.open(out_path, "w", **profile) as dst:
+            dst.write(red, 1)
+            dst.write(grn, 2)
+            dst.write(blu, 3)
+
+    return out_path
+
+def download_best_image(aoi_geojson: str,
+                        collection: str,
+                        date_range: Tuple[str, str],
+                        max_cloud: int,
+                        out_path: str,
+                        rgb_uint8: bool) -> Optional[str]:
+    """
+    Search STAC with a WGS84 AOI, pick a full-res asset:
+    - Prefer the 'visual' GeoTIFF (10 m), avoiding thumbnails.
+    - If not available, compose RGB from B04/B03/B02 (10 m).
+    """
+    geometry = read_aoi_geom_wgs84(aoi_geojson)
+
+    catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
+    search = catalog.search(
+        collections=[collection],
+        intersects=geometry,
+        datetime=f"{date_range[0]}/{date_range[1]}",
+        query={"eo:cloud_cover": {"lt": max_cloud}}
+    )
+    items = list(search.get_items())
+    if not items:
+        return None
+
+    # Least cloudy first
+    items.sort(key=lambda x: x.properties.get("eo:cloud_cover", 100))
+    item = items[0]
+
+    # 1) Try a proper full-res 'visual' GeoTIFF/COG
+    visual_asset = _pick_fullres_asset_or_none(item)
+    if visual_asset is not None:
+        signed_href = pc.sign(visual_asset.href)
+        r = requests.get(signed_href, stream=True, timeout=180)
+        r.raise_for_status()
+        with open(out_path, "wb") as f:
+            for chunk in r.iter_content(chunk_size=1 << 20):
+                f.write(chunk)
+        return out_path
+
+    # 2) Fall back to composing RGB from 10 m bands
+    try:
+        return _build_true_color_from_bands(item, out_path, rgb_uint8=rgb_uint8)
+    except Exception:
+        return None
+
+# --- Post-processing: resampling & overviews ----------------------------------
+
+def _pixel_size_from_transform(transform) -> Tuple[float, float]:
+    # (pixel width, pixel height) in CRS units (usually meters)
+    return abs(transform.a), abs(transform.e)
+
+def reproject_to_resolution(src_path: str, dst_path: str, target_res_m: float, force_upsample: bool=False) -> str:
+    """
+    Reproject/resample to a clean target pixel size (meters).
+    If target_res is finer than native and force_upsample=False, clamp to native.
+    """
+    with rasterio.open(src_path) as src:
+        native_x = abs(src.transform.a)
+        native_y = abs(src.transform.e)
+        native = (native_x + native_y) / 2.0
+
+        # Prevent unintentional upsampling (blur)
+        out_res = target_res_m
+        if (target_res_m < native) and not force_upsample:
+            out_res = native
+
+        dst_transform, width, height = calculate_default_transform(
+            src.crs, src.crs, src.width, src.height, *src.bounds, resolution=out_res
+        )
+        profile = src.profile.copy()
+        profile.update(transform=dst_transform, width=width, height=height)
+
+        # Choose resampling: nearest for upsampling (keeps edges), cubic for downsampling
+        upsampling = out_res < native
+        resamp = Resampling.nearest if upsampling else Resampling.cubic
+
+        with rasterio.open(dst_path, "w", **profile) as dst:
+            for b in range(1, src.count + 1):
+                reproject(
+                    source=rasterio.band(src, b),
+                    destination=rasterio.band(dst, b),
+                    src_transform=src.transform,
+                    src_crs=src.crs,
+                    dst_transform=dst_transform,
+                    dst_crs=src.crs,
+                    resampling=resamp
+                )
+    return dst_path
+
+def add_overviews_inplace(tif_path: str, levels=(2, 4, 8, 16)):
+    """
+    Build internal overviews to improve on-screen clarity at multiple zooms.
+    """
+    with rasterio.open(tif_path, "r+") as ds:
+        ds.build_overviews(levels, Resampling.average)
+        ds.update_tags(ns="rio_overview", resampling="average")
+
+# --- NDVI helpers -------------------------------------------------------------
+
+def collection_for_year(year: int):
+    """
+    Pick dataset + band names by year.
+    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",
+            "nir": "B08",
+            "scale": ("sentinel", None),  # 0..10000 -> reflectance 0..1
+            "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",
+            "nir": "SR_B5",
+            "scale": ("landsat_sr", (0.0000275, -0.2)),  # reflectance = DN*0.0000275 - 0.2
+            "native_res_m": 30.0,
+        }
+
+def _fetch_least_cloudy_item(collection: str, geometry: dict, date_range: tuple, max_cloud: int):
+    cat = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
+    search = cat.search(
+        collections=[collection],
+        intersects=geometry,
+        datetime=f"{date_range[0]}/{date_range[1]}",
+        query={"eo:cloud_cover": {"lt": max_cloud}}
+    )
+    items = list(search.get_items())
+    if not items:
+        return None
+    items.sort(key=lambda x: x.properties.get("eo:cloud_cover", 100))
+    return items[0]
+
+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:
+        arr = src.read(1)
+        profile = src.profile.copy()
+    return arr, profile
+
+def _scale_reflectance(arr: np.ndarray, scale_info) -> np.ndarray:
+    """Return reflectance in 0..1 range where possible."""
+    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:
+        return arr
+
+def _compute_ndvi(nir: np.ndarray, red: np.ndarray, nodata_mask: np.ndarray = None) -> np.ndarray:
+    denom = (nir + red)
+    ndvi = np.where(denom != 0, (nir - red) / denom, np.nan).astype("float32")
+    if nodata_mask is not None:
+        ndvi = np.where(nodata_mask, np.nan, ndvi)
+    return ndvi
+
+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"],
+        "width": profile["width"],
+        "count": 1,
+        "dtype": "float32",
+        "crs": profile.get("crs"),
+        "transform": profile.get("transform"),
+        "tiled": True,
+        "compress": "lzw",
+        "predictor": 2,
+        "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)
+
+
+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:
+        with rasterio.open(
+            mem_src, "w", driver="GTiff",
+            height=src_profile["height"], width=src_profile["width"],
+            count=1, dtype="float32",
+            crs=src_profile["crs"], transform=src_profile["transform"]
+        ) as src_ds:
+            src_ds.write(src_arr.astype("float32"), 1)
+
+            with MemoryFile() as mem_dst:
+                with rasterio.open(
+                    mem_dst, "w", driver="GTiff",
+                    height=ref_profile["height"], width=ref_profile["width"],
+                    count=1, dtype="float32",
+                    crs=ref_profile["crs"], transform=ref_profile["transform"]
+                ) as dst_ds:
+                    reproject(
+                        source=rasterio.band(src_ds, 1),
+                        destination=rasterio.band(dst_ds, 1),
+                        src_transform=src_profile["transform"],
+                        src_crs=src_profile["crs"],
+                        dst_transform=ref_profile["transform"],
+                        dst_crs=ref_profile["crs"],
+                        resampling=resampling
+                    )
+                    return dst_ds.read(1)
+
+def _summarize_ndvi(ndvi: np.ndarray, pixel_size_m: float) -> dict:
+    valid = np.isfinite(ndvi)
+    if not np.any(valid):
+        return {"count": 0}
+    v = ndvi[valid]
+    area_per_pixel_m2 = pixel_size_m * pixel_size_m
+    return {
+        "count": int(v.size),
+        "mean": float(np.nanmean(v)),
+        "median": float(np.nanmedian(v)),
+        "p05": float(np.nanpercentile(v, 5)),
+        "p95": float(np.nanpercentile(v, 95)),
+        "frac_gt_0_2": float(np.mean(v > 0.2)),
+        "frac_gt_0_4": float(np.mean(v > 0.4)),
+        "frac_gt_0_6": float(np.mean(v > 0.6)),
+        "area_gt_0_4_m2": float(np.sum(v > 0.4) * area_per_pixel_m2),
+        "area_gt_0_6_m2": float(np.sum(v > 0.6) * area_per_pixel_m2),
+    }
+
+# --- API setup ----------------------------------------------------------------
+
+app = FastAPI(
+    title="AOI Imagery API",
+    description="""
+An API that returns a satellite image file for a given year and computes NDVI & change.
+- Swagger UI: `/docs`
+- ReDoc: `/redoc`
+""",
+    version="1.3.0"
+)
+
+# Defaults; override via env or query params
+AOI_PATH = os.environ.get("AOI_PATH", "map.geojson")
+BUFFER_METERS_DEFAULT = float(os.environ.get("BUFFER_METERS", "200"))
+COLLECTION_DEFAULT = os.environ.get("COLLECTION", "sentinel-2-l2a")
+MAX_CLOUD_DEFAULT = int(os.environ.get("MAX_CLOUD", "10"))
+
+# Map a year -> date window (customize as needed; here June–Aug)
+def year_to_range(year: int) -> Tuple[str, str]:
+    return (f"{year}-06-01", f"{year}-08-31")
+
+# --- Download Endpoint --------------------------------------------------------
+
+class DownloadParams(BaseModel):
+    year: int = Field(..., ge=2015, le=2100, description="Year to fetch imagery for")
+    collection: Optional[str] = Field(default=COLLECTION_DEFAULT, description="STAC collection (e.g., sentinel-2-l2a)")
+    max_cloud: Optional[int] = Field(default=MAX_CLOUD_DEFAULT, ge=0, le=100, description="Max cloud cover percent")
+    buffer_meters: Optional[float] = Field(default=BUFFER_METERS_DEFAULT, ge=0, description="Buffer to apply to AOI (meters)")
+    clip_to_aoi: Optional[bool] = Field(default=False, description="If true, clip the image to buffered AOI")
+    rgb_uint8: Optional[bool] = Field(default=False, description="If true, write 8-bit RGB with a simple stretch (preview-friendly)")
+    target_res_m: Optional[float] = Field(default=10.0, ge=0, description="Reproject/resample clipped output to this pixel size in meters")
+    build_overviews: Optional[bool] = Field(default=True, description="If true, add internal overviews to the output GeoTIFF")
+
+@app.get(
+    "/download",
+    response_class=FileResponse,
+    responses={
+        200: {"content": {"image/tiff": {}}, "description": "GeoTIFF imagery file"},
+        404: {"description": "No imagery found for given inputs"},
+        422: {"description": "Validation error"},
+        500: {"description": "Server error"}
+    },
+    summary="Download imagery for a given year (GeoTIFF)"
+)
+def download_imagery(
+    year: int = Query(..., ge=2015, le=2100, description="Year to fetch imagery for"),
+    collection: str = Query(COLLECTION_DEFAULT, description="STAC collection (e.g., sentinel-2-l2a)"),
+    max_cloud: int = Query(MAX_CLOUD_DEFAULT, ge=0, le=100, description="Max cloud cover percent"),
+    buffer_meters: float = Query(BUFFER_METERS_DEFAULT, ge=0, description="Buffer in meters to apply to AOI"),
+    clip_to_aoi: bool = Query(False, description="If true, clip the image to buffered AOI"),
+    rgb_uint8: bool = Query(False, description="If true, write 8-bit RGB with a simple stretch (preview-friendly)"),
+    target_res_m: float = Query(10.0, ge=0, description="Reproject/resample clipped output to this pixel size in meters"),
+    build_overviews: bool = Query(True, description="If true, add internal overviews to the output GeoTIFF")
+):
+    """
+    Returns a GeoTIFF image file for the requested `year`.
+    By default returns the best (least cloudy) June–August image intersecting your AOI.
+    """
+    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)
+
+    date_range = year_to_range(year)
+    raw_tif = os.path.join(workdir, f"imagery_{year}.tif")
+    out_path = download_best_image(aoi_path, collection, date_range, max_cloud, raw_tif, rgb_uint8=rgb_uint8)
+    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 rasterio.open(out_path) as src:
+            geoms = gpd.read_file(aoi_path).to_crs(src.crs)
+            shapes = [mapping(geom) for geom in geoms.geometry]
+
+            out_image, out_transform = mask(src, shapes, crop=True)
+            out_meta = src.meta.copy()
+            out_meta.update({
+                "driver": "GTiff",
+                "height": out_image.shape[1],
+                "width": out_image.shape[2],
+                "transform": out_transform,
+                "photometric": "RGB",
+                "tiled": True,
+                "compress": "lzw",
+                "predictor": 2,
+            })
+
+            clipped_tif = os.path.join(workdir, f"imagery_{year}_clipped.tif")
+            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",
+        filename=os.path.basename(final_path)
+    )
+
+@app.post(
+    "/download",
+    response_class=FileResponse,
+    summary="Download imagery (POST body)"
+)
+def download_imagery_post(params: DownloadParams):
+    return download_imagery(
+        year=params.year,
+        collection=params.collection,
+        max_cloud=params.max_cloud,
+        buffer_meters=params.buffer_meters,
+        clip_to_aoi=params.clip_to_aoi,
+        rgb_uint8=params.rgb_uint8,
+        target_res_m=params.target_res_m,
+        build_overviews=params.build_overviews,
+    )
+
+# --- NDVI Endpoints -----------------------------------------------------------
+
+@app.get("/ndvi", summary="Compute NDVI for a given year; returns stats and writes a GeoTIFF")
+def ndvi_year(
+    year: int = Query(..., ge=2013, le=2100),
+    max_cloud: int = Query(MAX_CLOUD_DEFAULT, ge=0, le=100),
+    buffer_meters: float = Query(BUFFER_METERS_DEFAULT, ge=0),
+    clip_to_aoi: bool = Query(True, description="Clip NDVI to AOI"),
+    res_for_compare_m: float = Query(30.0, ge=0, description="Optional resampling resolution (m) for cross-year/sensor comparability"),
+):
+    if not os.path.exists(AOI_PATH):
+        raise HTTPException(status_code=500, detail=f"AOI file not found at {AOI_PATH}")
+
+    workdir = tempfile.mkdtemp(prefix=f"ndvi_{year}_")
+    aoi_path = AOI_PATH
+    if buffer_meters and buffer_meters > 0:
+        aoi_path = os.path.join(workdir, "aoi_buffered.geojson")
+        buffer_aoi(AOI_PATH, buffer_meters, aoi_path)
+
+    geometry = read_aoi_geom_wgs84(aoi_path)
+    cfg = collection_for_year(year)
+    date_range = year_to_range(year)
+
+    item = _fetch_least_cloudy_item(cfg["collection"], geometry, date_range, max_cloud)
+    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 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)
+
+    final_ndvi = tmp_ndvi
+    if clip_to_aoi:
+        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]
+            out_image, out_transform = mask(src, shapes, crop=True)
+            out_meta = src.meta.copy()
+            out_meta.update(transform=out_transform, height=out_image.shape[1], width=out_image.shape[2])
+            clipped = os.path.join(workdir, f"ndvi_{year}_clip.tif")
+            with rasterio.open(clipped, "w", **out_meta) as dst:
+                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 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
+            )
+            prof = src.profile.copy()
+            prof.update(transform=dst_transform, width=width, height=height)
+            repro = os.path.join(workdir, f"ndvi_{year}_{int(res_for_compare_m)}m.tif")
+            with rasterio.open(repro, "w", **prof) as dst:
+                reproject(
+                    source=rasterio.band(src, 1),
+                    destination=rasterio.band(dst, 1),
+                    src_transform=src.transform, src_crs=src.crs,
+                    dst_transform=dst_transform, dst_crs=src.crs,
+                    resampling=Resampling.bilinear
+                )
+            final_ndvi = repro
+
+    # Compute stats
+    with rasterio.open(final_ndvi) as ds:
+        arr = ds.read(1)
+        resx = abs(ds.transform.a)
+        stats = _summarize_ndvi(arr, resx)
+
+    return {
+        "year": year,
+        "collection": cfg["collection"],
+        "native_res_m": cfg["native_res_m"],
+        "compare_res_m": res_for_compare_m,
+        "ndvi_tif": os.path.basename(final_ndvi),
+        "workdir": workdir,
+        "stats": stats
+    }
+
+@app.get("/ndvi-change", summary="Compute NDVI change (end - start); returns stats and paths to GeoTIFFs")
+def ndvi_change(
+    start_year: int = Query(..., ge=2013),
+    end_year:   int = Query(..., ge=2013),
+    max_cloud: int = Query(MAX_CLOUD_DEFAULT, ge=0, le=100),
+    buffer_meters: float = Query(BUFFER_METERS_DEFAULT, ge=0),
+    clip_to_aoi: bool = Query(True),
+    compare_res_m: float = Query(30.0, ge=0, description="Common resolution (m) for both years"),
+):
+    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,
+        buffer_meters=buffer_meters,
+        clip_to_aoi=clip_to_aoi,
+        res_for_compare_m=compare_res_m
+    )
+    end = ndvi_year(
+        year=end_year,
+        max_cloud=max_cloud,
+        buffer_meters=buffer_meters,
+        clip_to_aoi=clip_to_aoi,
+        res_for_compare_m=compare_res_m
+    )
+
+    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 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)
+            profile = end_ds.profile.copy()
+        else:
+            arr_s = start_ds.read(1)
+            arr_e = end_ds.read(1)
+            profile = end_ds.profile.copy()
+
+        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)
+
+        resx = abs(profile["transform"].a)
+        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 = {
+            "mean_change": float(np.nanmean(ch)) if ch.size else None,
+            "median_change": float(np.nanmedian(ch)) if ch.size else None,
+            "p05_change": float(np.nanpercentile(ch, 5)) if ch.size else None,
+            "p95_change": float(np.nanpercentile(ch, 95)) if ch.size else None,
+            "frac_change_gt_0_1": float(np.mean(ch > 0.1)) if ch.size else None,
+            "frac_change_lt_-0_1": float(np.mean(ch < -0.1)) if ch.size else None,
+        }
+
+    return {
+        "start_year": start_year,
+        "end_year": end_year,
+        "compare_res_m": compare_res_m,
+        "ndvi_start_tif": os.path.basename(ndvi_start_path),
+        "ndvi_end_tif": os.path.basename(ndvi_end_path),
+        "ndvi_change_tif": os.path.basename(change_tif),
+        "outputs_dir": workdir,
+        "start_stats": stats_start,
+        "end_stats": stats_end,
+        "change_stats": change_stats
+    }
+
+# Run:
+# uvicorn main:app --host 0.0.0.0 --port 8000

BackendNew/map.geojson

+{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"coordinates":[[[9.16145374315991,48.78558930650067],[9.16145374315991,48.77109004651268],[9.18585591426276,48.77109004651268],[9.18585591426276,48.78558930650067],[9.16145374315991,48.78558930650067]]],"type":"Polygon"}}]}
\ No newline at end of file

BackendNew/preprocess_aoi.py

+# Preprocess AOI and Download Satellite Imagery from Microsoft Planetary Computer
+# Requirements: pip install geopandas rasterio pystac-client planetary-computer requests
+
+import geopandas as gpd
+import rasterio
+from rasterio.mask import mask
+import json
+import os
+from pystac_client import Client
+import planetary_computer as pc
+import requests
+import matplotlib.pyplot as plt
+
+# --- CONFIG ---
+AOI_PATH = "map.geojson"  # Your AOI file (exported from geojson.io or Google Earth)
+BUFFER_METERS = 200           # Buffer distance in meters
+OUTPUT_DIR = "processed"     # Output directory for processed files
+COLLECTION = "sentinel-2-l2a"  # or "landsat-8-c2-l2"
+DATE_RANGE = [
+    ("2016-06-01", "2016-08-31"),  # Past
+    ("2025-06-01", "2025-08-31")   # Present
+]
+MAX_CLOUD = 10  # %
+
+# --- STEP 1: Buffer AOI ---
+def buffer_aoi(aoi_path, buffer_meters, output_path):
+    aoi = gpd.read_file(aoi_path)
+    aoi_buffered = aoi.to_crs(epsg=3857).buffer(buffer_meters).to_crs(aoi.crs)
+    aoi_buffered_gdf = gpd.GeoDataFrame(geometry=aoi_buffered)
+    aoi_buffered_gdf.to_file(output_path, driver="GeoJSON")
+    print(f"Buffered AOI saved to {output_path}")
+    return output_path
+
+# --- STEP 2: Download Image from MPC ---
+def download_best_image(aoi_geojson, collection, date_range, max_cloud, out_path):
+    # Load AOI geometry
+    with open(aoi_geojson) as f:
+        geojson = json.load(f)
+    geometry = geojson['features'][0]['geometry']
+    # Search MPC STAC
+    catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")
+    search = catalog.search(
+        collections=[collection],
+        intersects=geometry,
+        datetime=f"{date_range[0]}/{date_range[1]}",
+        query={"eo:cloud_cover": {"lt": max_cloud}}
+    )
+    items = list(search.get_items())
+    if not items:
+        print(f"No images found for {date_range}")
+        return None
+    # Pick the least cloudy image
+    items.sort(key=lambda x: x.properties.get("eo:cloud_cover", 100))
+    item = items[0]
+    asset = item.assets["visual"] if "visual" in item.assets else list(item.assets.values())[0]
+    signed_href = pc.sign(asset.href)
+    # Download the image
+    print(f"Downloading {signed_href} ...")
+    r = requests.get(signed_href, stream=True)
+    with open(out_path, "wb") as f:
+        for chunk in r.iter_content(chunk_size=8192):
+            f.write(chunk)
+    print(f"Downloaded image to {out_path}")
+    return out_path
+
+# --- STEP 3: Display Images Side by Side ---
+def show_images_side_by_side(image_paths):
+    fig, axes = plt.subplots(1, 2, figsize=(16, 8))
+    for i, img_path in enumerate(image_paths):
+        with rasterio.open(img_path) as src:
+            img = src.read([1, 2, 3])  # RGB bands
+            img = img.transpose(1, 2, 0)
+            # Normalize for display
+            img = (img - img.min()) / (img.max() - img.min())
+            axes[i].imshow(img)
+            axes[i].set_title(f"Image {i+1}")
+            axes[i].axis('off')
+    plt.tight_layout()
+    plt.show()
+
+if __name__ == "__main__":
+    os.makedirs(OUTPUT_DIR, exist_ok=True)
+    # Download images only, skip buffer and clipping
+    image_paths = []
+    for i, drange in enumerate(DATE_RANGE):
+        out_img = os.path.join(OUTPUT_DIR, f"raw_{i+1}.tif")
+        img_path = download_best_image(AOI_PATH, COLLECTION, drange, MAX_CLOUD, out_img)
+        if img_path:
+            image_paths.append(img_path)
+    if len(image_paths) == 2:
+        show_images_side_by_side(image_paths)
+    print("\nAll done! You now have raw rasters in the 'processed' folder and can view them side by side.")

BackendNew/processed/map_buffered.geojson

+{
+"type": "FeatureCollection",
+"name": "map_buffered",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 9.166994209422283, 48.784091968852827 ], [ 9.178154201230427, 48.784091968852827 ], [ 9.178330301820953, 48.784086268546858 ], [ 9.178504706466526, 48.78406922252217 ], [ 9.178675735555073, 48.784040994929839 ], [ 9.178841741982971, 48.784001857598447 ], [ 9.179001127017576, 48.783952187416872 ], [ 9.179152355693873, 48.783892462705673 ], [ 9.179293971597041, 48.783823258611697 ], [ 9.179424610888514, 48.783745241570486 ], [ 9.179543015440506, 48.783659162889563 ], [ 9.179648044952451, 48.783565851514489 ], [ 9.179738687932733, 48.783466206047081 ], [ 9.179814071439905, 48.783361186092904 ], [ 9.179873469489584, 48.783251803021074 ], [ 9.179916310046078, 48.783139110225356 ], [ 9.179942180531409, 48.783024192980392 ], [ 9.179950831798665, 48.782908157990647 ], [ 9.179950831798665, 48.778299983780748 ], [ 9.179942180531409, 48.778183937869152 ], [ 9.179916310046078, 48.77806900927871 ], [ 9.179873469489584, 48.777956304846001 ], [ 9.179814071439905, 48.777846909994913 ], [ 9.179738687932733, 48.777741878282676 ], [ 9.179648044952451, 48.777642221252542 ], [ 9.179543015440506, 48.777548898690974 ], [ 9.179424610888514, 48.777462809382932 ], [ 9.179293971597041, 48.777384782454767 ], [ 9.179152355693873, 48.77731556938766 ], [ 9.179001127017576, 48.777255836778998 ], [ 9.178841741982971, 48.777206159921221 ], [ 9.178675735555073, 48.777167017260076 ], [ 9.178504706466526, 48.777138785785802 ], [ 9.178330301820953, 48.777121737401508 ], [ 9.178154201230427, 48.777116036303902 ], [ 9.166994209422283, 48.777116036303902 ], [ 9.166818108831755, 48.777121737401508 ], [ 9.166643704186182, 48.777138785785802 ], [ 9.166472675097637, 48.777167017260076 ], [ 9.166306668669737, 48.777206159921221 ], [ 9.166147283635134, 48.777255836778998 ], [ 9.165996054958837, 48.77731556938766 ], [ 9.16585443905567, 48.777384782454767 ], [ 9.165723799764194, 48.777462809382932 ], [ 9.165605395212204, 48.777548898690974 ], [ 9.165500365700257, 48.777642221252542 ], [ 9.165409722719977, 48.777741878282676 ], [ 9.165334339212803, 48.777846909994913 ], [ 9.165274941163126, 48.777956304846001 ], [ 9.165232100606632, 48.77806900927871 ], [ 9.165206230121299, 48.778183937869152 ], [ 9.165197578854045, 48.778299983780748 ], [ 9.165197578854045, 48.782908157990647 ], [ 9.165206230121299, 48.783024192980392 ], [ 9.165232100606632, 48.783139110225356 ], [ 9.165274941163126, 48.783251803021074 ], [ 9.165334339212803, 48.783361186092904 ], [ 9.165409722719977, 48.783466206047081 ], [ 9.165500365700257, 48.783565851514489 ], [ 9.165605395212204, 48.783659162889563 ], [ 9.165723799764194, 48.783745241570486 ], [ 9.16585443905567, 48.783823258611697 ], [ 9.165996054958837, 48.783892462705673 ], [ 9.166147283635134, 48.783952187416872 ], [ 9.166306668669737, 48.784001857598447 ], [ 9.166472675097637, 48.784040994929839 ], [ 9.166643704186182, 48.78406922252217 ], [ 9.166818108831755, 48.784086268546858 ], [ 9.166994209422283, 48.784091968852827 ] ] ] } }
+]
+}