Skip to content
GitLab
Commit 36353266 authored by Heisenberg5124's avatar Heisenberg5124
Browse files

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

# Conflicts:
#	BackendNew/main.py
parents fabaa42f 5cbbcb10
Hide whitespace changes
Inline Side-by-side
BackendNew/main.py
View file @ 36353266
......@@ -6,6 +6,7 @@
import os
import json
import tempfile
import time
from typing import Optional, Tuple
from fastapi import FastAPI, Query, HTTPException
......@@ -18,13 +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 fastapi.middleware.cors import CORSMiddleware
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 --------------------------------------------------------
......@@ -40,7 +53,6 @@ def read_aoi_geom_wgs84(aoi_path: str) -> dict:
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),
......@@ -56,7 +68,6 @@ def buffer_aoi(aoi_path: str, buffer_meters: float, output_path: str) -> str:
aoi_buffered_gdf.to_file(output_path, driver="GeoJSON")
return output_path
# --- High-quality asset selection / composition -------------------------------
def _pick_fullres_asset_or_none(item):
......@@ -75,7 +86,6 @@ def _pick_fullres_asset_or_none(item):
return a
return None
def _scale_to_uint8(arr: np.ndarray) -> np.ndarray:
"""
Simple 2–98 percentile contrast stretch to uint8.
......@@ -90,7 +100,6 @@ def _scale_to_uint8(arr: np.ndarray) -> np.ndarray:
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).
......@@ -104,38 +113,38 @@ 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 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",
)
with GDAL_ENV:
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")
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)
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],
......@@ -181,15 +190,13 @@ def download_best_image(aoi_geojson: str,
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:
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.
......@@ -227,7 +234,6 @@ def reproject_to_resolution(src_path: str, dst_path: str, target_res_m: float, f
)
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.
......@@ -236,7 +242,6 @@ def add_overviews_inplace(tif_path: str, levels=(2, 4, 8, 16)):
ds.build_overviews(levels, Resampling.average)
ds.update_tags(ns="rio_overview", resampling="average")
# --- NDVI helpers -------------------------------------------------------------
def collection_for_year(year: int):
......@@ -263,7 +268,6 @@ def collection_for_year(year: int):
"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(
......@@ -278,15 +282,15 @@ def _fetch_least_cloudy_item(collection: str, geometry: dict, date_range: tuple,
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()
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
def _scale_reflectance(arr: np.ndarray, scale_info) -> np.ndarray:
"""Return reflectance in 0..1 range where possible."""
kind, params = scale_info
......@@ -301,7 +305,6 @@ def _scale_reflectance(arr: np.ndarray, scale_info) -> np.ndarray:
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")
......@@ -309,7 +312,6 @@ def _compute_ndvi(nir: np.ndarray, red: np.ndarray, nodata_mask: np.ndarray = No
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.
......@@ -342,24 +344,23 @@ 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:
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"]
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"]
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),
......@@ -372,7 +373,6 @@ def _resample_match(src_arr, src_profile, ref_profile, resampling=Resampling.bil
)
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):
......@@ -392,6 +392,103 @@ 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 ----------------------------------------------------------------
......@@ -402,18 +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"
)
app.add_middleware(
CORSMiddleware,
allow_origins=[
"http://localhost:5173",
"http://127.0.0.1:5173",
], # use ["*"] for local dev only
allow_credentials=True,
allow_methods=["*"],
allow_headers=["*"],
version="1.4.0"
)
# Defaults; override via env or query params
......@@ -422,28 +508,21 @@ 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)")
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")
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",
......@@ -457,15 +536,14 @@ class DownloadParams(BaseModel):
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")
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`.
......@@ -487,32 +565,32 @@ def download_imagery(
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}")
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)
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]
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")
......@@ -531,7 +609,6 @@ def download_imagery(
filename=os.path.basename(final_path)
)
@app.post(
"/download",
response_class=FileResponse,
......@@ -549,17 +626,15 @@ 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(
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"),
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}")
......@@ -589,16 +664,16 @@ def ndvi_year(
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)
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()
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
......@@ -607,35 +682,37 @@ def ndvi_year(
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
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]
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
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
)
final_ndvi = repro
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:
......@@ -653,15 +730,14 @@ def ndvi_year(
"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"),
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")
......@@ -683,43 +759,43 @@ 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"])
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,
}
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)
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,
......@@ -734,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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Dies ist die Gitlab-Instanz des Transferportals der Hochschule für Technik Stuttgart. Hier geht es zurück zum Portal