Harmon/insar-management-system-v2
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Integrate GF3 SARscape flood workflow

Browse Source
This commit is contained in:
Harmon
2026-06-02 01:25:42 +08:00
parent 16cac0e292
commit 79e08b3a47
24 changed files with 4220 additions and 84 deletions
Download Patch File Download Diff File Expand all files Collapse all files
backend/app/services/water_detect_service.py
+159 -21
View File
@@ -10,6 +10,7 @@ import logging
import math
import os
import struct
from pathlib import Path
from typing import Any, Dict, Optional, Tuple
import numpy as np
@@ -21,6 +22,65 @@ logger = logging.getLogger(__name__)
# SRTM HGT helpers
# ---------------------------------------------------------------------------
def _otsu_threshold(values: np.ndarray, bins: int = 512) -> float:
"""Compute Otsu's threshold without depending on scikit-image."""
finite_values = values[np.isfinite(values)]
if finite_values.size == 0:
raise ValueError("No finite pixels available for Otsu threshold")
min_value = float(np.min(finite_values))
max_value = float(np.max(finite_values))
if min_value == max_value:
return min_value
counts, edges = np.histogram(finite_values, bins=bins, range=(min_value, max_value))
centers = (edges[:-1] + edges[1:]) / 2.0
total = float(counts.sum())
if total <= 0:
return min_value
weight_background = np.cumsum(counts).astype(np.float64)
weight_foreground = total - weight_background
cumulative_mean = np.cumsum(counts * centers)
total_mean = cumulative_mean[-1]
valid = (weight_background > 0) & (weight_foreground > 0)
between = np.zeros_like(centers, dtype=np.float64)
mean_background = np.zeros_like(centers, dtype=np.float64)
mean_foreground = np.zeros_like(centers, dtype=np.float64)
mean_background[valid] = cumulative_mean[valid] / weight_background[valid]
mean_foreground[valid] = (total_mean - cumulative_mean[valid]) / weight_foreground[valid]
between[valid] = (
weight_background[valid]
* weight_foreground[valid]
* (mean_background[valid] - mean_foreground[valid]) ** 2
)
return float(centers[int(np.argmax(between))])
def _disk_structure(radius: int) -> np.ndarray:
y, x = np.ogrid[-radius: radius + 1, -radius: radius + 1]
return (x * x + y * y) <= radius * radius
def _prepare_detection_image(img: np.ndarray, valid: np.ndarray) -> tuple[np.ndarray, np.ndarray, str]:
"""Normalize SAR values for water thresholding and keep the original mask."""
working = img.astype(np.float32, copy=True)
working[~valid] = np.nan
valid_values = working[valid]
if valid_values.size == 0:
return working, valid, "raw"
min_value = float(np.nanmin(valid_values))
p99_value = float(np.nanpercentile(valid_values, 99))
max_value = float(np.nanmax(valid_values))
if min_value >= 0.0 and (p99_value > 1.0 or max_value > 5.0):
positive_valid = valid & (working > 0)
converted = np.full_like(working, np.nan, dtype=np.float32)
converted[positive_valid] = 10.0 * np.log10(np.maximum(working[positive_valid], 1e-12))
return converted, positive_valid, "linear_to_db"
return working, valid, "raw"
def _read_hgt(filepath: str) -> np.ndarray:
"""Read a single SRTM .hgt file. Auto-detect SRTM1 (3601) vs SRTM3 (1201)."""
file_size = os.path.getsize(filepath)
@@ -96,6 +156,70 @@ def _load_srtm3_dem(
return mosaic, dem_bounds
def _candidate_dem_paths(dem_path: str) -> list[str]:
text = str(dem_path or "").strip()
if not text:
return []
path = Path(text)
candidates: list[Path] = []
for suffix in (".vrt", ".wgs84.vrt", ".tif", ".tiff", ".img", ".wgs84"):
candidate = Path(text + suffix)
if candidate.exists() and candidate not in candidates:
candidates.append(candidate)
if path.is_file() and path not in candidates:
candidates.append(path)
if path.is_dir():
for pattern in ("*.vrt", "*.tif", "*.tiff", "*.img"):
for candidate in path.glob(pattern):
if candidate not in candidates:
candidates.append(candidate)
return [str(candidate) for candidate in candidates]
def _load_raster_dem(
bounds: Tuple[float, float, float, float],
dem_path: str,
out_shape: tuple[int, int],
) -> Optional[np.ndarray]:
"""Read a DEM raster subset and resample it to the SAR image grid size."""
import rasterio
from rasterio.enums import Resampling
from rasterio.windows import from_bounds
height, width = out_shape
for candidate in _candidate_dem_paths(dem_path):
try:
with rasterio.open(candidate) as src:
if src.crs and not src.crs.is_geographic:
continue
dem_bounds = src.bounds
min_lon, min_lat, max_lon, max_lat = bounds
if (
max_lon <= dem_bounds.left
or min_lon >= dem_bounds.right
or max_lat <= dem_bounds.bottom
or min_lat >= dem_bounds.top
):
continue
window = from_bounds(min_lon, min_lat, max_lon, max_lat, transform=src.transform)
data = src.read(
1,
window=window,
out_shape=(height, width),
boundless=True,
fill_value=np.nan,
resampling=Resampling.bilinear,
).astype(np.float32)
nodata = src.nodata
if nodata is not None and np.isfinite(nodata):
data[data == np.float32(nodata)] = np.nan
logger.info("[WaterDetect] DEM raster loaded: %s", candidate)
return data
except Exception as exc:
logger.warning("[WaterDetect] DEM raster candidate skipped: %s (%s)", candidate, exc)
return None
# ---------------------------------------------------------------------------
# Core detection
# ---------------------------------------------------------------------------
@@ -119,14 +243,12 @@ def run_water_detection(
Returns dict with keys: ok, output_path, water_area_km2, water_pixel_count, otsu_threshold_db
"""
import rasterio
from scipy import ndimage
from scipy.ndimage import median_filter, gaussian_filter, label, zoom
from skimage.filters import threshold_otsu
from skimage.morphology import disk, binary_dilation, binary_erosion
from skimage.measure import regionprops
from ..config import settings
dem_dir = settings.SRTM_DEM_DIR
dem_path = settings.GF3_SARSCAPE_DEM_PATH or settings.GF3_GEO_DEM_PATH or settings.SRTM_DEM_DIR
os.makedirs(output_dir, exist_ok=True)
logger.info("[WaterDetect] Reading input: %s", geo_tiff_path)
@@ -136,6 +258,7 @@ def run_water_detection(
img = src.read(1).astype(np.float32)
transform = src.transform
crs = src.crs
nodata = src.nodata
height, width = img.shape
pixel_size_x = transform.a # degrees per pixel (x)
pixel_size_y = transform.e # degrees per pixel (y, negative)
@@ -155,17 +278,23 @@ def run_water_detection(
# Step 5: Valid mask
valid = np.isfinite(img) & (img != 0)
if nodata is not None and np.isfinite(nodata):
valid &= img != np.float32(nodata)
if np.count_nonzero(valid) < 100:
return {"ok": False, "error": "Too few valid pixels in input image"}
detection_img, valid, value_transform = _prepare_detection_image(img, valid)
logger.info("[WaterDetect] Value transform: %s", value_transform)
# Step 6: Otsu threshold on valid pixels
valid_pixels = img[valid]
thresh = threshold_otsu(valid_pixels)
valid_pixels = detection_img[valid]
thresh = _otsu_threshold(valid_pixels)
logger.info("[WaterDetect] Otsu threshold: %.4f", thresh)
# Step 7-8: Median + Gaussian filtering
filtered = median_filter(img, size=3)
filtered_input = np.where(valid, detection_img, np.nanmedian(valid_pixels)).astype(np.float32)
filtered = median_filter(filtered_input, size=3)
filtered = gaussian_filter(filtered, sigma=1.0)
# Step 9: Initial water mask
@@ -173,8 +302,11 @@ def run_water_detection(
# Step 3-4: Load and resample DEM (if available)
dem_applied = False
if dem_dir and os.path.isdir(dem_dir):
dem, dem_bounds = _load_srtm3_dem(bounds, dem_dir)
if dem_path and (os.path.isdir(dem_path) or os.path.isfile(dem_path)):
dem_resampled = _load_raster_dem(bounds, dem_path, (height, width))
dem, dem_bounds = (None, None)
if dem_resampled is None and os.path.isdir(dem_path):
dem, dem_bounds = _load_srtm3_dem(bounds, dem_path)
if dem is not None and dem_bounds is not None:
# Resample DEM to image resolution
zoom_y = height / dem.shape[0]
@@ -183,6 +315,7 @@ def run_water_detection(
# Clip to match image shape exactly
dem_resampled = dem_resampled[:height, :width]
if dem_resampled is not None:
# Step 10: DEM height constraint (0m <= DEM <= 1000m)
dem_valid = np.isfinite(dem_resampled)
height_mask = dem_valid & (dem_resampled >= 0) & (dem_resampled <= 1000)
@@ -204,14 +337,14 @@ def run_water_detection(
else:
logger.warning("[WaterDetect] DEM not available, skipping DEM constraints")
else:
logger.warning("[WaterDetect] SRTM_DEM_DIR not configured, skipping DEM constraints")
logger.warning("[WaterDetect] DEM path not configured, skipping DEM constraints")
# Step 12: Morphological processing — disk(5) dilate→erode→dilate→erode
selem = disk(5)
water = binary_dilation(water, selem)
water = binary_erosion(water, selem)
water = binary_dilation(water, selem)
water = binary_erosion(water, selem)
selem = _disk_structure(5)
water = ndimage.binary_dilation(water, structure=selem)
water = ndimage.binary_erosion(water, structure=selem)
water = ndimage.binary_dilation(water, structure=selem)
water = ndimage.binary_erosion(water, structure=selem)
# Step 13: Connected component filtering
labeled, num_features = label(water)
@@ -220,18 +353,22 @@ def run_water_detection(
pixel_area_m2 = abs(pixel_size_x) * 111320 * abs(pixel_size_y) * 111320
min_pixels_by_area = max(1, int(3000 / max(pixel_area_m2, 1)))
props = regionprops(labeled)
areas = [p.area for p in props]
if areas:
areas = ndimage.sum(
np.ones_like(labeled, dtype=np.uint8),
labeled,
index=np.arange(1, num_features + 1),
)
areas = np.asarray(areas, dtype=np.float64)
if areas.size:
median_area = float(np.median(areas))
min_area = max(min_pixels_by_area, int(median_area))
else:
min_area = min_pixels_by_area
# Remove small components
for prop in props:
if prop.area < min_area:
water[labeled == prop.label] = False
small_labels = np.where(areas < min_area)[0] + 1
if small_labels.size:
water[np.isin(labeled, small_labels)] = False
logger.info("[WaterDetect] Connected component filter: min_area=%d pixels, kept %d/%d components",
min_area, np.count_nonzero(np.unique(labeled[water])), num_features)
@@ -265,4 +402,5 @@ def run_water_detection(
"water_area_km2": round(water_area, 4),
"water_pixel_count": water_pixel_count,
"otsu_threshold_db": round(float(thresh), 4),
"value_transform": value_transform,
}