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

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

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