diff --git a/.env.example b/.env.example index d1bf001..de07c7e 100644 --- a/.env.example +++ b/.env.example @@ -119,6 +119,10 @@ IDL_DINSAR_DEM_BASE_FILE=D:\SRTM30m\SRTMDEM_RSP_SARscape SRTM_DEM_DIR=D:\SRTM30m GF3_GEO_DEM_PATH=D:\DEM\gf3_dem.jp2 WATER_RESULTS_DIR=D:\WaterResult +GF3_WATER_DLTB_CACHE_DIR=D:\production_assets\gf3_water\priors\dltb_cache\heilongjiang +GF3_WATER_DEM_PATH= +GF3_WATER_DEFAULT_CARTOGRAPHIC=true +GF3_WATER_DEFAULT_OUT_VECTOR=true SAR_ANALYSIS_READY_ROOT=D:\production_results\sar_analysis_ready SAR_ANALYSIS_WORK_ROOT=D:\production_runtime\sar_analysis_work SAR_ANALYSIS_NODATA_VALUE=-9999 diff --git a/backend/app/config.py b/backend/app/config.py index 958e3d6..020be19 100644 --- a/backend/app/config.py +++ b/backend/app/config.py @@ -215,6 +215,10 @@ class Settings(BaseSettings): RADAR_PREVIEW_BUILD_ON_DEMAND: bool = True WATER_RESULTS_DIR: str = "" + GF3_WATER_DLTB_CACHE_DIR: str = "" + GF3_WATER_DEM_PATH: str = "" + GF3_WATER_DEFAULT_CARTOGRAPHIC: bool = True + GF3_WATER_DEFAULT_OUT_VECTOR: bool = True SAR_ANALYSIS_READY_ROOT: str = "" SAR_ANALYSIS_WORK_ROOT: str = "" SAR_ANALYSIS_NODATA_VALUE: float = -9999.0 diff --git a/backend/app/processors/__init__.py b/backend/app/processors/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/backend/app/processors/gf3_water/__init__.py b/backend/app/processors/gf3_water/__init__.py new file mode 100644 index 0000000..04948b3 --- /dev/null +++ b/backend/app/processors/gf3_water/__init__.py @@ -0,0 +1,8 @@ +"""GF-3 HH/HV water extraction package.""" + +from .api import run_water_extraction +from .cli import main +from .config import WaterExtractionConfig +from .pipeline import run_from_args + +__all__ = ["WaterExtractionConfig", "main", "run_from_args", "run_water_extraction"] diff --git a/backend/app/processors/gf3_water/api.py b/backend/app/processors/gf3_water/api.py new file mode 100644 index 0000000..3996015 --- /dev/null +++ b/backend/app/processors/gf3_water/api.py @@ -0,0 +1,25 @@ +"""Programmatic API for GF-3 water extraction.""" + +from __future__ import annotations + +from argparse import Namespace +from dataclasses import asdict +from pathlib import Path + +from .config import WaterExtractionConfig +from .pipeline import run_from_args + + +def run_water_extraction(config: WaterExtractionConfig) -> int: + """Run the processor from a typed config object. + + The return value matches the CLI process exit code. Outputs are written under + ``config.out_dir`` and optional vector output paths. + """ + data = asdict(config) + for key, value in list(data.items()): + if isinstance(value, Path): + data[key] = value + elif isinstance(value, list): + data[key] = [Path(item) for item in value] + return run_from_args(Namespace(**data)) diff --git a/backend/app/processors/gf3_water/cli.py b/backend/app/processors/gf3_water/cli.py new file mode 100644 index 0000000..07afa33 --- /dev/null +++ b/backend/app/processors/gf3_water/cli.py @@ -0,0 +1,64 @@ +"""Command-line interface for GF-3 water extraction.""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +from .pipeline import run_from_args + + +def build_arg_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser(description="Extract water from GF-3 HH/HV ENVI images with a non-DL baseline.") + parser.add_argument("--hh", required=True, type=Path) + parser.add_argument("--hv", required=True, type=Path) + parser.add_argument("--out-dir", required=True, type=Path) + parser.add_argument("--dem", type=Path, default=None) + parser.add_argument("--dltb-gdb", type=Path, default=None, help="DLTB FileGDB path used as soft land-use prior.") + parser.add_argument("--dltb-cache-dir", type=Path, default=None, help="Directory containing water_prior.shp, paddy.shp, strict_review.shp DLTB cache layers.") + parser.add_argument("--dltb-layer", default="DLTB") + parser.add_argument("--dltb-field", default="DLMC") + parser.add_argument("--dltb-mode", choices=["soft", "strict", "off"], default="soft") + parser.add_argument("--dltb-max-features", type=int, default=None, help="Safety limit for DLTB features read from the source.") + parser.add_argument("--water-vector", type=Path, action="append", default=[], help="Known river/lake shapefile. Can be passed multiple times.") + parser.add_argument("--paddy-vector", type=Path, action="append", default=[], help="Paddy field/farmland water-sensitive vector. Can be passed multiple times.") + parser.add_argument("--river-buffer-meters", type=float, default=120.0, help="Buffer width for line river vectors.") + parser.add_argument("--paddy-buffer-meters", type=float, default=0.0, help="Buffer width for line paddy vectors, if any.") + parser.add_argument("--threshold-method", choices=["otsu", "percentile"], default="otsu") + parser.add_argument("--score-percentile", type=float, default=92.0) + parser.add_argument("--hv-percentile", type=float, default=35.0) + parser.add_argument("--prior-score-percentile", type=float, default=85.0) + parser.add_argument("--prior-hv-percentile", type=float, default=45.0) + parser.add_argument("--paddy-score-percentile", type=float, default=88.0) + parser.add_argument("--paddy-hv-percentile", type=float, default=45.0) + parser.add_argument("--candidate-score-percentile", type=float, default=90.0) + parser.add_argument("--candidate-hv-percentile", type=float, default=50.0) + parser.add_argument("--slope-max", type=float, default=8.0) + parser.add_argument("--close-pixels", type=int, default=0, help="Binary closing radius for final water mask before vectorization.") + parser.add_argument("--open-pixels", type=int, default=0, help="Binary opening radius for final water mask after hole filling.") + parser.add_argument("--paddy-close-pixels", type=int, default=0, help="Binary closing radius for paddy water-like mask.") + parser.add_argument("--paddy-open-pixels", type=int, default=0, help="Binary opening radius for paddy water-like mask.") + parser.add_argument("--candidate-close-pixels", type=int, default=0, help="Binary closing radius for low-confidence candidate mask.") + parser.add_argument("--candidate-open-pixels", type=int, default=0, help="Binary opening radius for low-confidence candidate mask.") + parser.add_argument("--cartographic-water", action="store_true", help="Export a map-production layer that merges high-confidence and review candidate water.") + parser.add_argument("--cartographic-include-paddy", action="store_true", help="Include paddy water-like pixels in cartographic water.") + parser.add_argument("--cartographic-close-pixels", type=int, default=3) + parser.add_argument("--cartographic-open-pixels", type=int, default=0) + parser.add_argument("--cartographic-fill-hole-pixels", type=int, default=4096) + parser.add_argument("--cartographic-min-component-pixels", type=int, default=512) + parser.add_argument("--min-component-pixels", type=int, default=128) + parser.add_argument("--fill-hole-pixels", type=int, default=512) + parser.add_argument("--no-morphology", action="store_true") + parser.add_argument("--out-vector-gpkg", type=Path, default=None, help="Write cartographic vector products to this GeoPackage.") + parser.add_argument("--out-vector-shp-dir", type=Path, default=None, help="Write one ESRI Shapefile per cartographic vector layer.") + parser.add_argument("--min-polygon-area-m2", type=float, default=1000.0) + parser.add_argument("--simplify-meters", type=float, default=0.0) + parser.add_argument("--smooth-meters", type=float, default=0.0, help="Vector boundary smoothing distance using buffer(+d)/buffer(-d).") + parser.add_argument("--min-hole-area-m2", type=float, default=0.0, help="Remove polygon interior holes smaller than this area.") + return parser + + +def main(argv: list[str] | None = None) -> int: + parser = build_arg_parser() + return run_from_args(parser.parse_args(argv)) + diff --git a/backend/app/processors/gf3_water/config.py b/backend/app/processors/gf3_water/config.py new file mode 100644 index 0000000..0e5d7af --- /dev/null +++ b/backend/app/processors/gf3_water/config.py @@ -0,0 +1,55 @@ +"""Configuration objects for GF-3 water extraction.""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path + + +@dataclass +class WaterExtractionConfig: + hh: Path + hv: Path + out_dir: Path + dem: Path | None = None + dltb_gdb: Path | None = None + dltb_cache_dir: Path | None = None + dltb_layer: str = "DLTB" + dltb_field: str = "DLMC" + dltb_mode: str = "soft" + dltb_max_features: int | None = None + water_vector: list[Path] = field(default_factory=list) + paddy_vector: list[Path] = field(default_factory=list) + river_buffer_meters: float = 120.0 + paddy_buffer_meters: float = 0.0 + threshold_method: str = "otsu" + score_percentile: float = 92.0 + hv_percentile: float = 35.0 + prior_score_percentile: float = 85.0 + prior_hv_percentile: float = 45.0 + paddy_score_percentile: float = 88.0 + paddy_hv_percentile: float = 45.0 + candidate_score_percentile: float = 90.0 + candidate_hv_percentile: float = 50.0 + slope_max: float = 8.0 + close_pixels: int = 0 + open_pixels: int = 0 + paddy_close_pixels: int = 0 + paddy_open_pixels: int = 0 + candidate_close_pixels: int = 0 + candidate_open_pixels: int = 0 + cartographic_water: bool = False + cartographic_include_paddy: bool = False + cartographic_close_pixels: int = 3 + cartographic_open_pixels: int = 0 + cartographic_fill_hole_pixels: int = 4096 + cartographic_min_component_pixels: int = 512 + min_component_pixels: int = 128 + fill_hole_pixels: int = 512 + no_morphology: bool = False + out_vector_gpkg: Path | None = None + out_vector_shp_dir: Path | None = None + min_polygon_area_m2: float = 1000.0 + simplify_meters: float = 0.0 + smooth_meters: float = 0.0 + min_hole_area_m2: float = 0.0 diff --git a/backend/app/processors/gf3_water/constants.py b/backend/app/processors/gf3_water/constants.py new file mode 100644 index 0000000..fdc14a9 --- /dev/null +++ b/backend/app/processors/gf3_water/constants.py @@ -0,0 +1,20 @@ +"""Class ids and product labels used by GF-3 water extraction.""" + +CLASS_NON_WATER = 0 +CLASS_HIGH_CONFIDENCE_WATER = 1 +CLASS_KNOWN_WATER = 2 +CLASS_PADDY_WATER_LIKE = 3 +CLASS_LOW_CONFIDENCE_WATER = 4 +CLASS_CARTOGRAPHIC_WATER = 5 +CLASS_INVALID = 255 + +CLASS_NAMES = { + CLASS_NON_WATER: "non_water", + CLASS_HIGH_CONFIDENCE_WATER: "high_confidence_water", + CLASS_KNOWN_WATER: "known_river_lake_water", + CLASS_PADDY_WATER_LIKE: "paddy_water_like", + CLASS_LOW_CONFIDENCE_WATER: "low_confidence_water", + CLASS_CARTOGRAPHIC_WATER: "cartographic_water", + CLASS_INVALID: "invalid", +} + diff --git a/backend/app/processors/gf3_water/dltb.py b/backend/app/processors/gf3_water/dltb.py new file mode 100644 index 0000000..bf512b4 --- /dev/null +++ b/backend/app/processors/gf3_water/dltb.py @@ -0,0 +1,146 @@ +"""DLTB land-use prior definitions. + +DLTB should be treated as a soft background layer for flood-period mapping: +it can change interpretation and confidence, but should not hard-exclude strong +water evidence outside normal water classes. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field as dc_field +from pathlib import Path + +import numpy as np +from shapely.geometry import box, shape +from shapely.ops import transform as shapely_transform + +try: + import fiona +except Exception: # pragma: no cover - optional at runtime + fiona = None + +try: + from pyproj import Transformer +except Exception: # pragma: no cover - optional at runtime + Transformer = None + + +@dataclass(frozen=True) +class DltbConfig: + gdb: Path + layer: str = "DLTB" + field: str = "DLMC" + mode: str = "soft" + max_features: int | None = None + water_names: tuple[str, ...] = dc_field(default_factory=lambda: ("河流水面", "湖泊水面", "水库水面", "坑塘水面", "沟渠")) + paddy_names: tuple[str, ...] = dc_field(default_factory=lambda: ("水田",)) + strict_names: tuple[str, ...] = dc_field(default_factory=lambda: ("城镇村道路用地", "公路用地", "农村道路", "设施农用地")) + + +@dataclass +class DltbSceneZones: + water_geoms: list + paddy_geoms: list + strict_geoms: list + normal_geoms: list + feature_count: int + class_counts: dict[str, int] + dlmc_values: list[str] + crs: str | None + + +def classify_dlmc(name: str, config: DltbConfig) -> str: + """Classify one DLMC value into a soft policy zone.""" + value = (name or "").strip() + if value in config.water_names: + return "water_prior" + if value in config.paddy_names: + return "paddy" + if value in config.strict_names: + return "strict_review" + return "normal" + + +def _transform_bounds(bounds: tuple[float, float, float, float], src_crs: str, dst_crs) -> tuple[float, float, float, float]: + if Transformer is None: + raise RuntimeError("pyproj is required to use DLTB priors with CRS transformation") + transformer = Transformer.from_crs(src_crs, dst_crs, always_xy=True) + left, bottom, right, top = bounds + xs = [left, left, right, right] + ys = [bottom, top, bottom, top] + tx, ty = transformer.transform(xs, ys) + return min(tx), min(ty), max(tx), max(ty) + + +def _transform_geom(geom, src_crs, dst_crs: str): + if Transformer is None: + raise RuntimeError("pyproj is required to transform DLTB geometries") + transformer = Transformer.from_crs(src_crs, dst_crs, always_xy=True) + return shapely_transform(lambda x, y, z=None: transformer.transform(np.asarray(x), np.asarray(y)), geom) + + +def load_dltb_scene_zones(config: DltbConfig, sar_bounds_wgs84: tuple[float, float, float, float]) -> DltbSceneZones: + """Load only DLTB features intersecting one SAR scene. + + The input bounds are WGS84 lon/lat. Returned geometries are transformed to + WGS84 and clipped to the SAR bounds. + """ + if fiona is None: + raise RuntimeError("fiona is required for DLTB geodatabase access") + if not config.gdb.exists(): + raise FileNotFoundError(f"DLTB geodatabase not found: {config.gdb}") + + roi_wgs84 = box(*sar_bounds_wgs84) + water_geoms = [] + paddy_geoms = [] + strict_geoms = [] + normal_geoms = [] + class_counts = {"water_prior": 0, "paddy": 0, "strict_review": 0, "normal": 0} + dlmc_values = set() + feature_count = 0 + + with fiona.open(config.gdb, layer=config.layer) as src: + src_crs = src.crs_wkt or src.crs + bbox = sar_bounds_wgs84 + if src_crs: + bbox = _transform_bounds(sar_bounds_wgs84, "EPSG:4326", src_crs) + + for feat in src.filter(bbox=bbox): + if config.max_features is not None and feature_count >= config.max_features: + break + geom_data = feat.get("geometry") + if not geom_data: + continue + geom = shape(geom_data) + if geom.is_empty: + continue + if src_crs: + geom = _transform_geom(geom, src_crs, "EPSG:4326") + geom = geom.intersection(roi_wgs84) + if geom.is_empty: + continue + + dlmc = str(feat.get("properties", {}).get(config.field, "") or "").strip() + zone = classify_dlmc(dlmc, config) + dlmc_values.add(dlmc) + class_counts[zone] += 1 + feature_count += 1 + if zone == "water_prior": + water_geoms.append(geom) + elif zone == "paddy": + paddy_geoms.append(geom) + elif zone == "strict_review": + strict_geoms.append(geom) + else: + normal_geoms.append(geom) + + return DltbSceneZones( + water_geoms=water_geoms, + paddy_geoms=paddy_geoms, + strict_geoms=strict_geoms, + normal_geoms=normal_geoms, + feature_count=feature_count, + class_counts=class_counts, + dlmc_values=sorted(v for v in dlmc_values if v), + crs=str(src_crs) if "src_crs" in locals() else None, + ) diff --git a/backend/app/processors/gf3_water/envi.py b/backend/app/processors/gf3_water/envi.py new file mode 100644 index 0000000..adb7b64 --- /dev/null +++ b/backend/app/processors/gf3_water/envi.py @@ -0,0 +1,188 @@ +"""ENVI raster parsing and GeoTIFF writing.""" + +from __future__ import annotations + +import math +import re +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +try: + import rasterio + from rasterio.crs import CRS + from rasterio.transform import Affine +except Exception: # pragma: no cover - optional at runtime + rasterio = None + CRS = None + Affine = None + + +ENVI_DTYPES = { + 1: np.uint8, + 2: np.int16, + 3: np.int32, + 4: np.float32, + 5: np.float64, + 12: np.uint16, + 13: np.uint32, + 14: np.int64, + 15: np.uint64, +} + + +@dataclass(frozen=True) +class EnviInfo: + path: Path + hdr_path: Path + samples: int + lines: int + bands: int + header_offset: int + dtype: np.dtype + byte_order: int + interleave: str + x0: float + y0: float + dx: float + dy: float + crs_wkt: str | None + + @property + def bounds(self) -> tuple[float, float, float, float]: + left = self.x0 + top = self.y0 + right = left + self.samples * self.dx + bottom = top - self.lines * self.dy + return left, bottom, right, top + + @property + def transform(self): + if Affine is None: + return None + return Affine(self.dx, 0.0, self.x0, 0.0, -self.dy, self.y0) + + +def read_hdr_text(data_path: Path) -> tuple[Path, str]: + hdr_path = data_path.with_suffix(data_path.suffix + ".hdr") if data_path.suffix else Path(str(data_path) + ".hdr") + if not hdr_path.exists(): + alt = data_path.with_suffix(".hdr") + if alt.exists(): + hdr_path = alt + if not hdr_path.exists(): + raise FileNotFoundError(f"ENVI header not found for {data_path}") + return hdr_path, hdr_path.read_text(encoding="utf-8", errors="ignore") + + +def hdr_value(text: str, key: str, default: str | None = None) -> str: + match = re.search(rf"(?im)^\s*{re.escape(key)}\s*=\s*(.+?)\s*$", text) + if match: + return match.group(1).strip() + if default is not None: + return default + raise ValueError(f"Missing ENVI header key: {key}") + + +def parse_map_info(text: str) -> tuple[float, float, float, float]: + match = re.search(r"(?is)map info\s*=\s*\{(.+?)\}", text) + if not match: + raise ValueError("Missing map info in ENVI header") + parts = [p.strip() for p in match.group(1).replace("\n", " ").split(",")] + if len(parts) < 7: + raise ValueError(f"Unexpected map info: {match.group(0)}") + return float(parts[3]), float(parts[4]), abs(float(parts[5])), abs(float(parts[6])) + + +def parse_crs_wkt(text: str) -> str | None: + match = re.search(r"(?is)coordinate system string\s*=\s*\{(.+?)\}", text) + return match.group(1).strip() if match else None + + +def parse_envi(path: Path) -> EnviInfo: + hdr_path, text = read_hdr_text(path) + dtype_code = int(hdr_value(text, "data type")) + if dtype_code not in ENVI_DTYPES: + raise ValueError(f"Unsupported ENVI data type: {dtype_code}") + x0, y0, dx, dy = parse_map_info(text) + dtype = np.dtype(ENVI_DTYPES[dtype_code]) + byte_order = int(hdr_value(text, "byte order", "0")) + if byte_order == 1: + dtype = dtype.newbyteorder(">") + else: + dtype = dtype.newbyteorder("<") + return EnviInfo( + path=path, + hdr_path=hdr_path, + samples=int(hdr_value(text, "samples")), + lines=int(hdr_value(text, "lines")), + bands=int(hdr_value(text, "bands", "1")), + header_offset=int(hdr_value(text, "header offset", "0")), + dtype=dtype, + byte_order=byte_order, + interleave=hdr_value(text, "interleave", "bsq").lower(), + x0=x0, + y0=y0, + dx=dx, + dy=dy, + crs_wkt=parse_crs_wkt(text), + ) + + +def read_envi_band(info: EnviInfo) -> np.ndarray: + if info.bands != 1 or info.interleave != "bsq": + raise ValueError("This baseline expects one-band BSQ ENVI inputs") + count = info.lines * info.samples + data = np.memmap(info.path, dtype=info.dtype, mode="r", offset=info.header_offset, shape=(count,)) + arr = np.asarray(data.reshape(info.lines, info.samples), dtype=np.float32) + arr = arr.copy() + arr[~np.isfinite(arr)] = np.nan + return arr + + +def read_dem_for_sar(dem_info: EnviInfo, sar_info: EnviInfo) -> np.ndarray: + left, bottom, right, top = sar_info.bounds + pad = 2 + col0 = max(0, int(math.floor((left - dem_info.x0) / dem_info.dx)) - pad) + col1 = min(dem_info.samples, int(math.ceil((right - dem_info.x0) / dem_info.dx)) + pad) + row0 = max(0, int(math.floor((dem_info.y0 - top) / dem_info.dy)) - pad) + row1 = min(dem_info.lines, int(math.ceil((dem_info.y0 - bottom) / dem_info.dy)) + pad) + if col1 <= col0 or row1 <= row0: + raise ValueError("SAR image does not overlap DEM") + + mm = np.memmap(dem_info.path, dtype=dem_info.dtype, mode="r", offset=dem_info.header_offset, shape=(dem_info.lines, dem_info.samples)) + dem_window = np.asarray(mm[row0:row1, col0:col1], dtype=np.float32).copy() + dem_window[~np.isfinite(dem_window)] = np.nan + + x = sar_info.x0 + (np.arange(sar_info.samples) + 0.5) * sar_info.dx + y = sar_info.y0 - (np.arange(sar_info.lines) + 0.5) * sar_info.dy + dem_cols = np.clip(np.rint((x - dem_info.x0) / dem_info.dx - 0.5).astype(np.int64) - col0, 0, dem_window.shape[1] - 1) + dem_rows = np.clip(np.rint((dem_info.y0 - y) / dem_info.dy - 0.5).astype(np.int64) - row0, 0, dem_window.shape[0] - 1) + return dem_window[dem_rows[:, None], dem_cols[None, :]] + + +def write_tif(path: Path, arr: np.ndarray, info: EnviInfo, dtype: str, nodata=None) -> None: + if rasterio is None: + return + crs = None + if CRS is not None: + try: + crs = CRS.from_epsg(4326) + except Exception: + crs = None + profile = { + "driver": "GTiff", + "height": arr.shape[0], + "width": arr.shape[1], + "count": 1, + "dtype": dtype, + "compress": "deflate", + "predictor": 2 if dtype.startswith("float") else 1, + "transform": info.transform, + "nodata": nodata, + } + if crs is not None: + profile["crs"] = crs + with rasterio.open(path, "w", **profile) as dst: + dst.write(arr.astype(dtype), 1) + diff --git a/backend/app/processors/gf3_water/geo.py b/backend/app/processors/gf3_water/geo.py new file mode 100644 index 0000000..3116717 --- /dev/null +++ b/backend/app/processors/gf3_water/geo.py @@ -0,0 +1,25 @@ +"""Small geospatial math helpers for lon/lat products.""" + +from __future__ import annotations + +import math + + +def meters_to_degrees(meters: float, lat: float) -> tuple[float, float]: + deg_lat = meters / 111_320.0 + deg_lon = meters / max(111_320.0 * math.cos(math.radians(lat)), 1.0) + return deg_lon, deg_lat + + +def pixel_area_m2(info) -> float: + center_lat = info.y0 - info.lines * info.dy * 0.5 + meters_per_deg_lat = 111_320.0 + meters_per_deg_lon = 111_320.0 * math.cos(math.radians(center_lat)) + return abs(info.dx * meters_per_deg_lon * info.dy * meters_per_deg_lat) + + +def polygon_area_m2(geom, center_lat: float) -> float: + meters_per_deg_lat = 111_320.0 + meters_per_deg_lon = 111_320.0 * math.cos(math.radians(center_lat)) + return float(abs(geom.area) * meters_per_deg_lon * meters_per_deg_lat) + diff --git a/backend/app/processors/gf3_water/pipeline.py b/backend/app/processors/gf3_water/pipeline.py new file mode 100644 index 0000000..de501cb --- /dev/null +++ b/backend/app/processors/gf3_water/pipeline.py @@ -0,0 +1,404 @@ +"""High-level GF-3 HH/HV water extraction pipeline.""" + +from __future__ import annotations + +import json +from argparse import Namespace + +import numpy as np + +from .constants import ( + CLASS_HIGH_CONFIDENCE_WATER, + CLASS_INVALID, + CLASS_KNOWN_WATER, + CLASS_LOW_CONFIDENCE_WATER, + CLASS_NAMES, + CLASS_NON_WATER, + CLASS_PADDY_WATER_LIKE, +) +from .dltb import DltbConfig, load_dltb_scene_zones +from .envi import parse_envi, read_dem_for_sar, read_envi_band, write_tif +from .previews import save_class_preview, save_dltb_zone_preview, save_gray_png, save_mask_png, save_preview +from .raster_ops import close_mask, fill_small_holes, open_mask, otsu_threshold, remove_small_components, robust_normalize, slope_degrees, to_db +from .vector_io import rasterize_geometries, rasterize_vector_mask, rasterize_water_prior +from .vector_products import write_classified_vectors + + +def _load_dltb_masks(args: Namespace, info) -> tuple[dict | None, np.ndarray, np.ndarray, np.ndarray]: + water_mask = np.zeros((info.lines, info.samples), dtype=bool) + paddy_mask = np.zeros((info.lines, info.samples), dtype=bool) + strict_mask = np.zeros((info.lines, info.samples), dtype=bool) + + if args.dltb_cache_dir is not None and args.dltb_mode != "off": + water_cache = args.dltb_cache_dir / "water_prior.shp" + paddy_cache = args.dltb_cache_dir / "paddy.shp" + strict_cache = args.dltb_cache_dir / "strict_review.shp" + if water_cache.exists(): + water_mask = rasterize_vector_mask([water_cache], info, 0.0, "--dltb-cache-dir water_prior") + if paddy_cache.exists(): + paddy_mask = rasterize_vector_mask([paddy_cache], info, 0.0, "--dltb-cache-dir paddy") + if strict_cache.exists(): + strict_mask = rasterize_vector_mask([strict_cache], info, 0.0, "--dltb-cache-dir strict_review") + stats = { + "cache_dir": str(args.dltb_cache_dir), + "mode": args.dltb_mode, + "source": "cache", + "water_prior_path": str(water_cache) if water_cache.exists() else None, + "paddy_path": str(paddy_cache) if paddy_cache.exists() else None, + "strict_review_path": str(strict_cache) if strict_cache.exists() else None, + } + return stats, water_mask, paddy_mask, strict_mask + + if args.dltb_gdb is not None and args.dltb_mode != "off": + zones = load_dltb_scene_zones( + DltbConfig(gdb=args.dltb_gdb, layer=args.dltb_layer, field=args.dltb_field, mode=args.dltb_mode, max_features=args.dltb_max_features), + info.bounds, + ) + water_mask = rasterize_geometries(zones.water_geoms, info) + paddy_mask = rasterize_geometries(zones.paddy_geoms, info) + strict_mask = rasterize_geometries(zones.strict_geoms, info) + stats = { + "gdb": str(args.dltb_gdb), + "source": "gdb", + "layer": args.dltb_layer, + "field": args.dltb_field, + "mode": args.dltb_mode, + "source_crs": zones.crs, + "features_in_scene": zones.feature_count, + "class_counts": zones.class_counts, + "dlmc_values_in_scene": zones.dlmc_values, + } + return stats, water_mask, paddy_mask, strict_mask + + return None, water_mask, paddy_mask, strict_mask + + +def _ensure_matching_grids(hh_info, hv_info) -> None: + hh_grid = (hh_info.samples, hh_info.lines, hh_info.x0, hh_info.y0, hh_info.dx, hh_info.dy) + hv_grid = (hv_info.samples, hv_info.lines, hv_info.x0, hv_info.y0, hv_info.dx, hv_info.dy) + if hh_grid != hv_grid: + raise ValueError("HH and HV grids do not match") + + +def _classify_products(mask: np.ndarray, low_confidence_candidate: np.ndarray, paddy_candidate: np.ndarray, prior_mask: np.ndarray, valid: np.ndarray) -> np.ndarray: + high_confidence_mask = mask & ~paddy_candidate + known_water_mask = mask & prior_mask + classified = np.full(mask.shape, CLASS_INVALID, dtype=np.uint8) + classified[valid] = CLASS_NON_WATER + classified[low_confidence_candidate & valid & ~mask & ~paddy_candidate] = CLASS_LOW_CONFIDENCE_WATER + classified[high_confidence_mask & valid] = CLASS_HIGH_CONFIDENCE_WATER + classified[paddy_candidate & valid] = CLASS_PADDY_WATER_LIKE + classified[known_water_mask & valid] = CLASS_KNOWN_WATER + return classified + + +def _write_rasters( + args: Namespace, + info, + probability: np.ndarray, + valid: np.ndarray, + mask: np.ndarray, + raw_mask: np.ndarray, + classified: np.ndarray, + cartographic_water: np.ndarray | None, + prior_mask: np.ndarray, + paddy_mask: np.ndarray, + paddy_candidate: np.ndarray, + dltb_enabled: bool, + dltb_water_mask: np.ndarray, + dltb_paddy_mask: np.ndarray, + dltb_strict_mask: np.ndarray, +) -> None: + mask_u8 = np.where(valid, mask.astype(np.uint8), 255).astype(np.uint8) + raw_mask_u8 = np.where(valid, raw_mask.astype(np.uint8), 255).astype(np.uint8) + + write_tif(args.out_dir / "water_score.tif", np.where(np.isfinite(probability), probability, -9999.0), info, "float32", nodata=-9999.0) + write_tif(args.out_dir / "water_mask.tif", mask_u8, info, "uint8", nodata=255) + write_tif(args.out_dir / "water_mask_raw.tif", raw_mask_u8, info, "uint8", nodata=255) + write_tif(args.out_dir / "classified_water.tif", classified, info, "uint8", nodata=CLASS_INVALID) + if cartographic_water is not None: + write_tif(args.out_dir / "cartographic_water.tif", np.where(valid, cartographic_water.astype(np.uint8), 255).astype(np.uint8), info, "uint8", nodata=255) + if args.water_vector: + write_tif(args.out_dir / "known_water_prior.tif", prior_mask.astype(np.uint8), info, "uint8", nodata=0) + if args.paddy_vector or dltb_enabled: + write_tif(args.out_dir / "paddy_prior.tif", paddy_mask.astype(np.uint8), info, "uint8", nodata=0) + write_tif(args.out_dir / "paddy_water_like.tif", np.where(valid, paddy_candidate.astype(np.uint8), 255).astype(np.uint8), info, "uint8", nodata=255) + if dltb_enabled: + write_tif(args.out_dir / "dltb_water_prior.tif", dltb_water_mask.astype(np.uint8), info, "uint8", nodata=0) + write_tif(args.out_dir / "dltb_paddy_prior.tif", dltb_paddy_mask.astype(np.uint8), info, "uint8", nodata=0) + write_tif(args.out_dir / "dltb_strict_zone.tif", dltb_strict_mask.astype(np.uint8), info, "uint8", nodata=0) + + +def _write_previews( + args: Namespace, + probability: np.ndarray, + valid: np.ndarray, + mask: np.ndarray, + raw_mask: np.ndarray, + hh_norm: np.ndarray, + hv_norm: np.ndarray, + classified: np.ndarray, + prior_mask: np.ndarray, + paddy_candidate: np.ndarray, + dltb_enabled: bool, + dltb_water_mask: np.ndarray, + dltb_paddy_mask: np.ndarray, + dltb_strict_mask: np.ndarray, + cartographic_water: np.ndarray | None, +) -> None: + save_gray_png(args.out_dir / "water_score.png", probability, valid) + save_mask_png(args.out_dir / "water_mask.png", mask) + save_mask_png(args.out_dir / "water_mask_raw.png", raw_mask) + if args.water_vector: + save_mask_png(args.out_dir / "known_water_prior.png", prior_mask) + if args.paddy_vector or dltb_enabled: + save_mask_png(args.out_dir / "paddy_water_like.png", paddy_candidate) + if dltb_enabled: + save_dltb_zone_preview(args.out_dir / "dltb_zone_preview.png", dltb_water_mask, dltb_paddy_mask, dltb_strict_mask) + if cartographic_water is not None: + save_mask_png(args.out_dir / "cartographic_water.png", cartographic_water) + save_preview(args.out_dir / "preview_overlay.png", hh_norm, hv_norm, mask, valid) + save_class_preview(args.out_dir / "classified_preview.png", classified, hh_norm, hv_norm, valid) + + +def run_from_args(args: Namespace) -> int: + args.out_dir.mkdir(parents=True, exist_ok=True) + + hh_info = parse_envi(args.hh) + hv_info = parse_envi(args.hv) + _ensure_matching_grids(hh_info, hv_info) + + dltb_stats, dltb_water_mask, dltb_paddy_mask, dltb_strict_mask = _load_dltb_masks(args, hh_info) + dltb_enabled = dltb_stats is not None + + hh = read_envi_band(hh_info) + hv = read_envi_band(hv_info) + valid = np.isfinite(hh) & np.isfinite(hv) & (hh > 0) & (hv > 0) + + hh_db = to_db(hh) + hv_db = to_db(hv) + valid &= np.isfinite(hh_db) & np.isfinite(hv_db) + hh_norm, hh_lo, hh_hi = robust_normalize(hh_db, valid) + hv_norm, hv_lo, hv_hi = robust_normalize(hv_db, valid) + + low_backscatter = 1.0 - 0.5 * (hh_norm + hv_norm) + low_backscatter[~valid] = np.nan + + values = low_backscatter[valid] + if args.threshold_method == "otsu": + score_threshold = otsu_threshold(values) + else: + score_threshold = float(np.nanpercentile(values, args.score_percentile)) + hv_dark_threshold = float(np.nanpercentile(hv_norm[valid], args.hv_percentile)) + + mask = (low_backscatter >= score_threshold) & (hv_norm <= hv_dark_threshold) & valid + prior_mask = np.zeros(mask.shape, dtype=bool) + prior_candidate_pixels = 0 + prior_score_threshold = None + prior_hv_threshold = None + if args.water_vector: + prior_mask = rasterize_water_prior(args.water_vector, hh_info, args.river_buffer_meters) & valid + prior_score_threshold = float(np.nanpercentile(values, args.prior_score_percentile)) + prior_hv_threshold = float(np.nanpercentile(hv_norm[valid], args.prior_hv_percentile)) + prior_candidate = prior_mask & (low_backscatter >= prior_score_threshold) & (hv_norm <= prior_hv_threshold) + prior_candidate_pixels = int(prior_candidate.sum()) + mask |= prior_candidate + if dltb_enabled: + prior_mask |= dltb_water_mask & valid + prior_score_threshold = prior_score_threshold if prior_score_threshold is not None else float(np.nanpercentile(values, args.prior_score_percentile)) + prior_hv_threshold = prior_hv_threshold if prior_hv_threshold is not None else float(np.nanpercentile(hv_norm[valid], args.prior_hv_percentile)) + dltb_prior_candidate = prior_mask & (low_backscatter >= prior_score_threshold) & (hv_norm <= prior_hv_threshold) + prior_candidate_pixels += int((dltb_water_mask & dltb_prior_candidate).sum()) + mask |= dltb_prior_candidate + + paddy_mask = np.zeros(mask.shape, dtype=bool) + paddy_candidate = np.zeros(mask.shape, dtype=bool) + paddy_score_threshold = None + paddy_hv_threshold = None + if args.paddy_vector: + paddy_mask = rasterize_vector_mask(args.paddy_vector, hh_info, args.paddy_buffer_meters, "--paddy-vector") & valid + paddy_score_threshold = float(np.nanpercentile(values, args.paddy_score_percentile)) + paddy_hv_threshold = float(np.nanpercentile(hv_norm[valid], args.paddy_hv_percentile)) + paddy_candidate = paddy_mask & (low_backscatter >= paddy_score_threshold) & (hv_norm <= paddy_hv_threshold) + if dltb_enabled: + paddy_mask |= dltb_paddy_mask & valid + paddy_score_threshold = paddy_score_threshold if paddy_score_threshold is not None else float(np.nanpercentile(values, args.paddy_score_percentile)) + paddy_hv_threshold = paddy_hv_threshold if paddy_hv_threshold is not None else float(np.nanpercentile(hv_norm[valid], args.paddy_hv_percentile)) + paddy_candidate |= paddy_mask & (low_backscatter >= paddy_score_threshold) & (hv_norm <= paddy_hv_threshold) + + candidate_score_threshold = float(np.nanpercentile(values, args.candidate_score_percentile)) + candidate_hv_threshold = float(np.nanpercentile(hv_norm[valid], args.candidate_hv_percentile)) + low_confidence_candidate = (low_backscatter >= candidate_score_threshold) & (hv_norm <= candidate_hv_threshold) & valid + if dltb_enabled and args.dltb_mode == "soft": + strong_candidate = (low_backscatter >= score_threshold) & (hv_norm <= hv_dark_threshold) & valid + low_confidence_candidate &= ~dltb_strict_mask | strong_candidate + elif dltb_enabled and args.dltb_mode == "strict": + low_confidence_candidate &= ~dltb_strict_mask + + dem_used = False + slope_threshold = None + slope = None + if args.dem is not None: + dem_info = parse_envi(args.dem) + dem = read_dem_for_sar(dem_info, hh_info) + slope = slope_degrees(dem, hh_info.dx, hh_info.dy, center_lat=hh_info.y0 - hh_info.lines * hh_info.dy * 0.5) + slope_threshold = float(args.slope_max) + mask &= np.isfinite(slope) & (slope <= args.slope_max) + paddy_candidate &= np.isfinite(slope) & (slope <= args.slope_max) + low_confidence_candidate &= np.isfinite(slope) & (slope <= args.slope_max) + valid &= np.isfinite(slope) + dem_used = True + save_gray_png(args.out_dir / "slope_preview.png", slope, np.isfinite(slope)) + write_tif(args.out_dir / "slope_degrees.tif", np.where(np.isfinite(slope), slope, -9999.0), hh_info, "float32", nodata=-9999.0) + + raw_mask = mask.copy() + if not args.no_morphology: + mask = close_mask(mask, args.close_pixels, valid) + mask = remove_small_components(mask, args.min_component_pixels) + mask = fill_small_holes(mask, args.fill_hole_pixels) + mask = open_mask(mask, args.open_pixels, valid) + paddy_candidate = close_mask(paddy_candidate, args.paddy_close_pixels, valid) + paddy_candidate = open_mask(paddy_candidate, args.paddy_open_pixels, valid) + low_confidence_candidate = close_mask(low_confidence_candidate, args.candidate_close_pixels, valid) + low_confidence_candidate = open_mask(low_confidence_candidate, args.candidate_open_pixels, valid) + + probability = np.clip(low_backscatter, 0.0, 1.0) + probability[~valid] = np.nan + classified = _classify_products(mask, low_confidence_candidate, paddy_candidate, prior_mask, valid) + + cartographic_water = None + if args.cartographic_water: + cartographic_water = (mask | low_confidence_candidate) & valid + if args.cartographic_include_paddy: + cartographic_water |= paddy_candidate & valid + else: + cartographic_water &= ~paddy_candidate + if not args.no_morphology: + cartographic_water = close_mask(cartographic_water, args.cartographic_close_pixels, valid) + cartographic_water = fill_small_holes(cartographic_water, args.cartographic_fill_hole_pixels) + cartographic_water = remove_small_components(cartographic_water, args.cartographic_min_component_pixels) + cartographic_water = open_mask(cartographic_water, args.cartographic_open_pixels, valid) + + _write_rasters( + args, + hh_info, + probability, + valid, + mask, + raw_mask, + classified, + cartographic_water, + prior_mask, + paddy_mask, + paddy_candidate, + dltb_enabled, + dltb_water_mask, + dltb_paddy_mask, + dltb_strict_mask, + ) + _write_previews( + args, + probability, + valid, + mask, + raw_mask, + hh_norm, + hv_norm, + classified, + prior_mask, + paddy_candidate, + dltb_enabled, + dltb_water_mask, + dltb_paddy_mask, + dltb_strict_mask, + cartographic_water, + ) + + vector_stats = None + if args.out_vector_gpkg is not None or args.out_vector_shp_dir is not None: + vector_stats = write_classified_vectors( + args.out_vector_gpkg, + args.out_vector_shp_dir, + classified, + cartographic_water, + hh_info, + probability, + hh_db, + hv_db, + slope, + prior_mask, + paddy_mask, + args.min_polygon_area_m2, + args.simplify_meters, + args.smooth_meters, + args.min_hole_area_m2, + ) + + valid_count = int(valid.sum()) + stats = { + "hh": str(args.hh), + "hv": str(args.hv), + "dem": str(args.dem) if args.dem else None, + "dltb": dltb_stats, + "water_vectors": [str(p) for p in args.water_vector], + "paddy_vectors": [str(p) for p in args.paddy_vector], + "river_buffer_meters": float(args.river_buffer_meters), + "paddy_buffer_meters": float(args.paddy_buffer_meters), + "shape": [hh_info.lines, hh_info.samples], + "bounds_wgs84": list(hh_info.bounds), + "valid_pixels": valid_count, + "valid_ratio": float(valid_count / valid.size), + "hh_db_percentile_2_98": [hh_lo, hh_hi], + "hv_db_percentile_2_98": [hv_lo, hv_hi], + "threshold_method": args.threshold_method, + "score_threshold": float(score_threshold), + "hv_norm_dark_threshold": float(hv_dark_threshold), + "prior_score_threshold": prior_score_threshold, + "prior_hv_norm_dark_threshold": prior_hv_threshold, + "paddy_score_threshold": paddy_score_threshold, + "paddy_hv_norm_dark_threshold": paddy_hv_threshold, + "candidate_score_threshold": candidate_score_threshold, + "candidate_hv_norm_dark_threshold": candidate_hv_threshold, + "slope_threshold_degrees": slope_threshold, + "dem_used": dem_used, + "known_water_prior_pixels": int(prior_mask.sum()), + "known_water_prior_ratio_valid": float(prior_mask.sum() / max(valid_count, 1)), + "dltb_water_prior_pixels": int(dltb_water_mask.sum()), + "dltb_paddy_prior_pixels": int(dltb_paddy_mask.sum()), + "dltb_strict_zone_pixels": int(dltb_strict_mask.sum()), + "prior_candidate_pixels": prior_candidate_pixels, + "paddy_prior_pixels": int(paddy_mask.sum()), + "paddy_water_like_pixels": int(paddy_candidate.sum()), + "paddy_water_like_ratio_valid": float(paddy_candidate.sum() / max(valid_count, 1)), + "low_confidence_water_pixels": int((classified == CLASS_LOW_CONFIDENCE_WATER).sum()), + "cartographic_water_pixels": int(cartographic_water.sum()) if cartographic_water is not None else 0, + "cartographic_water_ratio_valid": float(cartographic_water.sum() / max(valid_count, 1)) if cartographic_water is not None else 0.0, + "raw_water_pixels": int(raw_mask.sum()), + "raw_water_ratio_valid": float(raw_mask.sum() / max(valid_count, 1)), + "water_pixels": int(mask.sum()), + "water_ratio_valid": float(mask.sum() / max(valid_count, 1)), + "classified_counts": {CLASS_NAMES[class_id]: int((classified == class_id).sum()) for class_id in CLASS_NAMES}, + "min_component_pixels": int(args.min_component_pixels), + "fill_hole_pixels": int(args.fill_hole_pixels), + "close_pixels": int(args.close_pixels), + "open_pixels": int(args.open_pixels), + "paddy_close_pixels": int(args.paddy_close_pixels), + "paddy_open_pixels": int(args.paddy_open_pixels), + "candidate_close_pixels": int(args.candidate_close_pixels), + "candidate_open_pixels": int(args.candidate_open_pixels), + "cartographic_water_enabled": bool(args.cartographic_water), + "cartographic_include_paddy": bool(args.cartographic_include_paddy), + "cartographic_close_pixels": int(args.cartographic_close_pixels), + "cartographic_open_pixels": int(args.cartographic_open_pixels), + "cartographic_fill_hole_pixels": int(args.cartographic_fill_hole_pixels), + "cartographic_min_component_pixels": int(args.cartographic_min_component_pixels), + "min_polygon_area_m2": float(args.min_polygon_area_m2), + "simplify_meters": float(args.simplify_meters), + "smooth_meters": float(args.smooth_meters), + "min_hole_area_m2": float(args.min_hole_area_m2), + "vector_output": vector_stats, + } + (args.out_dir / "metadata.json").write_text(json.dumps(stats, indent=2), encoding="utf-8") + print(json.dumps(stats, indent=2)) + return 0 + diff --git a/backend/app/processors/gf3_water/previews.py b/backend/app/processors/gf3_water/previews.py new file mode 100644 index 0000000..69c551b --- /dev/null +++ b/backend/app/processors/gf3_water/previews.py @@ -0,0 +1,67 @@ +"""PNG preview writers for extraction products.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from PIL import Image + +from .constants import ( + CLASS_HIGH_CONFIDENCE_WATER, + CLASS_KNOWN_WATER, + CLASS_LOW_CONFIDENCE_WATER, + CLASS_PADDY_WATER_LIKE, +) + + +def save_gray_png(path: Path, arr: np.ndarray, valid: np.ndarray) -> None: + out = np.zeros(arr.shape, dtype=np.uint8) + vals = arr[valid & np.isfinite(arr)] + if vals.size: + lo, hi = np.nanpercentile(vals, [2, 98]) + scaled = np.clip((arr - lo) / max(hi - lo, 1e-6), 0.0, 1.0) + out[valid & np.isfinite(arr)] = (scaled[valid & np.isfinite(arr)] * 255).astype(np.uint8) + Image.fromarray(out).save(path) + + +def save_mask_png(path: Path, mask: np.ndarray) -> None: + Image.fromarray(np.where(mask, 255, 0).astype(np.uint8)).save(path) + + +def save_preview(path: Path, hh_norm: np.ndarray, hv_norm: np.ndarray, mask: np.ndarray, valid: np.ndarray) -> None: + rgb = np.zeros((*mask.shape, 3), dtype=np.uint8) + rgb[..., 0] = np.nan_to_num(hh_norm * 255.0, nan=0.0).astype(np.uint8) + rgb[..., 1] = np.nan_to_num(hv_norm * 255.0, nan=0.0).astype(np.uint8) + rgb[..., 2] = np.nan_to_num((1.0 - 0.5 * (hh_norm + hv_norm)) * 255.0, nan=0.0).astype(np.uint8) + rgb[mask] = (0.35 * rgb[mask] + np.array([0, 120, 255]) * 0.65).astype(np.uint8) + rgb[~valid] = 0 + Image.fromarray(rgb).save(path) + + +def save_class_preview(path: Path, classified: np.ndarray, hh_norm: np.ndarray, hv_norm: np.ndarray, valid: np.ndarray) -> None: + rgb = np.zeros((*classified.shape, 3), dtype=np.uint8) + base = np.nan_to_num((0.55 * hh_norm + 0.45 * hv_norm) * 180.0, nan=0.0).astype(np.uint8) + rgb[..., 0] = base + rgb[..., 1] = base + rgb[..., 2] = base + colors = { + CLASS_HIGH_CONFIDENCE_WATER: np.array([0, 92, 230], dtype=np.uint8), + CLASS_KNOWN_WATER: np.array([0, 170, 255], dtype=np.uint8), + CLASS_PADDY_WATER_LIKE: np.array([0, 210, 170], dtype=np.uint8), + CLASS_LOW_CONFIDENCE_WATER: np.array([245, 166, 35], dtype=np.uint8), + } + for class_id, color in colors.items(): + idx = classified == class_id + rgb[idx] = (0.30 * rgb[idx] + 0.70 * color).astype(np.uint8) + rgb[~valid] = 0 + Image.fromarray(rgb).save(path) + + +def save_dltb_zone_preview(path: Path, water_mask: np.ndarray, paddy_mask: np.ndarray, strict_mask: np.ndarray) -> None: + rgb = np.zeros((*water_mask.shape, 3), dtype=np.uint8) + rgb[water_mask] = (0, 120, 255) + rgb[paddy_mask] = (0, 210, 170) + rgb[strict_mask] = (220, 80, 40) + Image.fromarray(rgb).save(path) + diff --git a/backend/app/processors/gf3_water/processor.py b/backend/app/processors/gf3_water/processor.py new file mode 100644 index 0000000..3e532cc --- /dev/null +++ b/backend/app/processors/gf3_water/processor.py @@ -0,0 +1,14 @@ +"""Backward-compatible processor facade. + +New code should import from ``gf3_water.pipeline`` or ``gf3_water.cli``. This +module remains so existing scripts and integrations that import +``gf3_water.processor`` continue to work. +""" + +from __future__ import annotations + +from .cli import build_arg_parser, main +from .pipeline import run_from_args + +__all__ = ["build_arg_parser", "main", "run_from_args"] + diff --git a/backend/app/processors/gf3_water/raster_ops.py b/backend/app/processors/gf3_water/raster_ops.py new file mode 100644 index 0000000..ea874bf --- /dev/null +++ b/backend/app/processors/gf3_water/raster_ops.py @@ -0,0 +1,98 @@ +"""Raster transforms, thresholding, and morphology.""" + +from __future__ import annotations + +import math + +import numpy as np +from scipy import ndimage as ndi + + +def slope_degrees(dem: np.ndarray, dx_deg: float, dy_deg: float, center_lat: float) -> np.ndarray: + meters_per_deg_lat = 111_320.0 + meters_per_deg_lon = 111_320.0 * math.cos(math.radians(center_lat)) + dz_dy, dz_dx = np.gradient(dem.astype(np.float32), dy_deg * meters_per_deg_lat, dx_deg * meters_per_deg_lon) + slope = np.degrees(np.arctan(np.sqrt(dz_dx * dz_dx + dz_dy * dz_dy))) + slope[~np.isfinite(slope)] = np.nan + return slope.astype(np.float32) + + +def to_db(arr: np.ndarray, eps: float = 1e-8) -> np.ndarray: + out = np.full(arr.shape, np.nan, dtype=np.float32) + valid = np.isfinite(arr) & (arr > 0) + out[valid] = 10.0 * np.log10(arr[valid] + eps) + return out + + +def robust_normalize(arr: np.ndarray, valid: np.ndarray, q_low: float = 2.0, q_high: float = 98.0) -> tuple[np.ndarray, float, float]: + values = arr[valid] + lo, hi = np.nanpercentile(values, [q_low, q_high]) + if not np.isfinite(lo) or not np.isfinite(hi) or hi <= lo: + lo, hi = float(np.nanmin(values)), float(np.nanmax(values)) + norm = np.clip((arr - lo) / max(hi - lo, 1e-6), 0.0, 1.0) + norm[~valid] = np.nan + return norm.astype(np.float32), float(lo), float(hi) + + +def otsu_threshold(values: np.ndarray, bins: int = 512) -> float: + values = values[np.isfinite(values)] + if values.size == 0: + raise ValueError("No finite values for thresholding") + hist, edges = np.histogram(values, bins=bins) + centers = (edges[:-1] + edges[1:]) * 0.5 + weight1 = np.cumsum(hist).astype(np.float64) + weight2 = np.cumsum(hist[::-1]).astype(np.float64)[::-1] + mean1 = np.cumsum(hist * centers) / np.maximum(weight1, 1.0) + mean2 = (np.cumsum((hist * centers)[::-1]) / np.maximum(weight2[::-1], 1.0))[::-1] + variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2 + return float(centers[:-1][np.argmax(variance12)]) + + +def remove_small_components(mask: np.ndarray, min_pixels: int) -> np.ndarray: + if min_pixels <= 1: + return mask + labels, count = ndi.label(mask) + if count == 0: + return mask + sizes = np.bincount(labels.ravel()) + keep = sizes >= min_pixels + keep[0] = False + return keep[labels] + + +def fill_small_holes(mask: np.ndarray, max_pixels: int) -> np.ndarray: + if max_pixels <= 0: + return mask + inv = ~mask + labels, count = ndi.label(inv) + if count == 0: + return mask + border = np.unique(np.concatenate([labels[0, :], labels[-1, :], labels[:, 0], labels[:, -1]])) + sizes = np.bincount(labels.ravel()) + fill = sizes <= max_pixels + fill[border] = False + out = mask.copy() + out[fill[labels]] = True + return out + + +def disk_structure(radius: int) -> np.ndarray: + if radius <= 0: + return np.ones((1, 1), dtype=bool) + y, x = np.ogrid[-radius : radius + 1, -radius : radius + 1] + return (x * x + y * y) <= radius * radius + + +def close_mask(mask: np.ndarray, radius: int, valid: np.ndarray) -> np.ndarray: + if radius <= 0: + return mask + closed = ndi.binary_closing(mask & valid, structure=disk_structure(radius)) + return closed & valid + + +def open_mask(mask: np.ndarray, radius: int, valid: np.ndarray) -> np.ndarray: + if radius <= 0: + return mask + opened = ndi.binary_opening(mask & valid, structure=disk_structure(radius)) + return opened & valid + diff --git a/backend/app/processors/gf3_water/vector_io.py b/backend/app/processors/gf3_water/vector_io.py new file mode 100644 index 0000000..9f69c11 --- /dev/null +++ b/backend/app/processors/gf3_water/vector_io.py @@ -0,0 +1,88 @@ +"""Vector loading and rasterization helpers.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from shapely.geometry import box, shape +from shapely.ops import transform as shapely_transform + +from .envi import EnviInfo +from .geo import meters_to_degrees + +try: + import fiona +except Exception: # pragma: no cover - optional at runtime + fiona = None + +try: + from rasterio.features import rasterize +except Exception: # pragma: no cover - optional at runtime + rasterize = None + + +def load_vector_geometries(paths: list[Path], bounds: tuple[float, float, float, float], line_buffer_meters: float) -> list: + if not paths: + return [] + if fiona is None: + raise RuntimeError("fiona is required for vector inputs") + + roi = box(*bounds) + center_lat = (bounds[1] + bounds[3]) * 0.5 + buffer_lon, buffer_lat = meters_to_degrees(max(line_buffer_meters, 0.1), center_lat) + search_roi = box(bounds[0] - buffer_lon, bounds[1] - buffer_lat, bounds[2] + buffer_lon, bounds[3] + buffer_lat) + geoms = [] + + for path in paths: + with fiona.open(path) as src: + for feat in src: + if not feat.get("geometry"): + continue + geom = shape(feat["geometry"]) + if geom.is_empty or not geom.intersects(search_roi): + continue + if geom.geom_type in ("LineString", "MultiLineString"): + geom = shapely_transform(lambda x, y, z=None: (np.asarray(x) / buffer_lon, np.asarray(y) / buffer_lat), geom) + geom = geom.buffer(1.0) + geom = shapely_transform(lambda x, y, z=None: (np.asarray(x) * buffer_lon, np.asarray(y) * buffer_lat), geom) + geom = geom.intersection(search_roi) + if not geom.is_empty and geom.intersects(roi): + geoms.append(geom) + return geoms + + +def rasterize_vector_mask(paths: list[Path], info: EnviInfo, line_buffer_meters: float, label: str) -> np.ndarray: + if rasterize is None: + raise RuntimeError(f"rasterio.features.rasterize is required for {label}") + geoms = load_vector_geometries(paths, info.bounds, line_buffer_meters) + if not geoms: + return np.zeros((info.lines, info.samples), dtype=bool) + return rasterize( + [(geom, 1) for geom in geoms], + out_shape=(info.lines, info.samples), + transform=info.transform, + fill=0, + dtype="uint8", + all_touched=True, + ).astype(bool) + + +def rasterize_water_prior(paths: list[Path], info: EnviInfo, river_buffer_meters: float) -> np.ndarray: + return rasterize_vector_mask(paths, info, river_buffer_meters, "--water-vector") + + +def rasterize_geometries(geoms: list, info: EnviInfo) -> np.ndarray: + if rasterize is None: + raise RuntimeError("rasterio.features.rasterize is required to rasterize geometry priors") + if not geoms: + return np.zeros((info.lines, info.samples), dtype=bool) + return rasterize( + [(geom, 1) for geom in geoms if not geom.is_empty], + out_shape=(info.lines, info.samples), + transform=info.transform, + fill=0, + dtype="uint8", + all_touched=True, + ).astype(bool) + diff --git a/backend/app/processors/gf3_water/vector_products.py b/backend/app/processors/gf3_water/vector_products.py new file mode 100644 index 0000000..02f30e7 --- /dev/null +++ b/backend/app/processors/gf3_water/vector_products.py @@ -0,0 +1,269 @@ +"""Cartographic vector product writing.""" + +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +from shapely.geometry import MultiPolygon, Polygon, shape +from shapely.ops import transform as shapely_transform + +from .constants import ( + CLASS_CARTOGRAPHIC_WATER, + CLASS_HIGH_CONFIDENCE_WATER, + CLASS_KNOWN_WATER, + CLASS_LOW_CONFIDENCE_WATER, + CLASS_NAMES, + CLASS_NON_WATER, + CLASS_PADDY_WATER_LIKE, +) +from .envi import EnviInfo +from .geo import meters_to_degrees, pixel_area_m2, polygon_area_m2 + +try: + import fiona +except Exception: # pragma: no cover - optional at runtime + fiona = None + +try: + from rasterio.features import shapes +except Exception: # pragma: no cover - optional at runtime + shapes = None + + +def remove_small_polygon_holes(geom, min_hole_area_m2: float, center_lat: float): + if min_hole_area_m2 <= 0 or geom.is_empty: + return geom + + def clean_polygon(poly: Polygon) -> Polygon: + interiors = [] + for ring in poly.interiors: + hole = Polygon(ring) + if polygon_area_m2(hole, center_lat) >= min_hole_area_m2: + interiors.append(ring) + return Polygon(poly.exterior, interiors) + + if geom.geom_type == "Polygon": + return clean_polygon(geom) + if geom.geom_type == "MultiPolygon": + parts = [clean_polygon(poly) for poly in geom.geoms if not poly.is_empty] + return MultiPolygon(parts) if parts else geom + return geom + + +def smooth_geometry_meters(geom, smooth_meters: float, center_lat: float): + if smooth_meters <= 0 or geom.is_empty: + return geom + smooth_lon, smooth_lat = meters_to_degrees(smooth_meters, center_lat) + scaled = shapely_transform(lambda x, y, z=None: (np.asarray(x) / smooth_lon, np.asarray(y) / smooth_lat), geom) + smoothed = scaled.buffer(1.0).buffer(-1.0) + return shapely_transform(lambda x, y, z=None: (np.asarray(x) * smooth_lon, np.asarray(y) * smooth_lat), smoothed) + + +def write_classified_vectors( + gpkg_path: Path | None, + shp_dir: Path | None, + classified: np.ndarray, + cartographic_water: np.ndarray | None, + info: EnviInfo, + score: np.ndarray, + hh_db: np.ndarray, + hv_db: np.ndarray, + slope: np.ndarray | None, + known_water: np.ndarray, + paddy: np.ndarray, + min_area_m2: float, + simplify_meters: float, + smooth_meters: float, + min_hole_area_m2: float, +) -> dict: + if fiona is None or shapes is None: + raise RuntimeError("fiona and rasterio.features.shapes are required for vector output") + + if gpkg_path is not None: + gpkg_path.parent.mkdir(parents=True, exist_ok=True) + if gpkg_path.exists(): + gpkg_path.unlink() + if shp_dir is not None: + shp_dir.mkdir(parents=True, exist_ok=True) + + px_area = pixel_area_m2(info) + center_lat = info.y0 - info.lines * info.dy * 0.5 + simplify_lon, _ = meters_to_degrees(simplify_meters, center_lat) + schema = { + "geometry": "Polygon", + "properties": { + "class_id": "int", + "class_name": "str:32", + "confidence": "str:16", + "area_m2": "float", + "pixels": "int", + "mean_score": "float", + "mean_hh_db": "float", + "mean_hv_db": "float", + "mean_slope": "float", + "known_water": "int", + "paddy": "int", + "review_flag": "int", + }, + } + crs = "EPSG:4326" + counts = {name: 0 for name in CLASS_NAMES.values()} + layers = { + CLASS_CARTOGRAPHIC_WATER: "cartographic_water", + CLASS_HIGH_CONFIDENCE_WATER: "high_confidence_water", + CLASS_KNOWN_WATER: "known_water", + CLASS_PADDY_WATER_LIKE: "paddy_water_like", + CLASS_LOW_CONFIDENCE_WATER: "review_candidates", + } + handles = {} + shp_handles = {} + try: + if gpkg_path is not None: + for class_id, layer_name in layers.items(): + handles[class_id] = fiona.open(gpkg_path, "w", driver="GPKG", layer=layer_name, crs=crs, schema=schema) + if shp_dir is not None: + for class_id, layer_name in layers.items(): + shp_path = shp_dir / f"{layer_name}.shp" + for suffix in (".shp", ".shx", ".dbf", ".prj", ".cpg"): + sidecar = shp_path.with_suffix(suffix) + if sidecar.exists(): + sidecar.unlink() + shp_handles[class_id] = fiona.open(shp_path, "w", driver="ESRI Shapefile", crs=crs, schema=schema, encoding="UTF-8") + + vector_classes = classified.copy() + class_mask = np.isin(vector_classes, [CLASS_HIGH_CONFIDENCE_WATER, CLASS_KNOWN_WATER, CLASS_PADDY_WATER_LIKE, CLASS_LOW_CONFIDENCE_WATER]) + if cartographic_water is not None: + class_mask |= cartographic_water + for geom_mapping, value in shapes(vector_classes.astype(np.uint8), mask=class_mask, transform=info.transform): + class_id = int(value) + if cartographic_water is not None and class_id == CLASS_NON_WATER: + continue + if class_id not in layers: + continue + geom = shape(geom_mapping) + if geom.is_empty: + continue + geom = smooth_geometry_meters(geom, smooth_meters, center_lat) + geom = remove_small_polygon_holes(geom, min_hole_area_m2, center_lat) + if geom.is_empty: + continue + if simplify_meters > 0: + geom = geom.simplify(simplify_lon, preserve_topology=True) + if geom.is_empty: + continue + geom_area_m2 = polygon_area_m2(geom, center_lat) + + minx, miny, maxx, maxy = geom.bounds + col0 = max(0, int(math.floor((minx - info.x0) / info.dx)) - 1) + col1 = min(info.samples, int(math.ceil((maxx - info.x0) / info.dx)) + 1) + row0 = max(0, int(math.floor((info.y0 - maxy) / info.dy)) - 1) + row1 = min(info.lines, int(math.ceil((info.y0 - miny) / info.dy)) + 1) + if col1 <= col0 or row1 <= row0: + continue + if class_id == CLASS_CARTOGRAPHIC_WATER and cartographic_water is not None: + window = cartographic_water[row0:row1, col0:col1] + else: + window = classified[row0:row1, col0:col1] == class_id + pixels = int(window.sum()) + area_m2 = pixels * px_area + if max(area_m2, geom_area_m2) < min_area_m2: + continue + + score_window = score[row0:row1, col0:col1] + hh_window = hh_db[row0:row1, col0:col1] + hv_window = hv_db[row0:row1, col0:col1] + slope_window = slope[row0:row1, col0:col1] if slope is not None else None + mean_slope = float(np.nanmean(slope_window[window])) if slope_window is not None and np.any(np.isfinite(slope_window[window])) else -9999.0 + review_flag = 1 if class_id in (CLASS_PADDY_WATER_LIKE, CLASS_LOW_CONFIDENCE_WATER, CLASS_CARTOGRAPHIC_WATER) else 0 + confidence = "map" if class_id == CLASS_CARTOGRAPHIC_WATER else ("high" if class_id in (CLASS_HIGH_CONFIDENCE_WATER, CLASS_KNOWN_WATER) else "review") + feature = { + "geometry": geom.__geo_interface__, + "properties": { + "class_id": class_id, + "class_name": CLASS_NAMES[class_id], + "confidence": confidence, + "area_m2": float(geom_area_m2), + "pixels": pixels, + "mean_score": float(np.nanmean(score_window[window])), + "mean_hh_db": float(np.nanmean(hh_window[window])), + "mean_hv_db": float(np.nanmean(hv_window[window])), + "mean_slope": mean_slope, + "known_water": int(np.any(known_water[row0:row1, col0:col1] & window)), + "paddy": int(np.any(paddy[row0:row1, col0:col1] & window)), + "review_flag": review_flag, + }, + } + if class_id in handles: + handles[class_id].write(feature) + if class_id in shp_handles: + shp_handles[class_id].write(feature) + counts[CLASS_NAMES[class_id]] += 1 + + if cartographic_water is not None: + for geom_mapping, value in shapes(np.where(cartographic_water, CLASS_CARTOGRAPHIC_WATER, CLASS_NON_WATER).astype(np.uint8), mask=cartographic_water, transform=info.transform): + class_id = int(value) + geom = shape(geom_mapping) + if geom.is_empty: + continue + geom = smooth_geometry_meters(geom, smooth_meters, center_lat) + geom = remove_small_polygon_holes(geom, min_hole_area_m2, center_lat) + if geom.is_empty: + continue + if simplify_meters > 0: + geom = geom.simplify(simplify_lon, preserve_topology=True) + if geom.is_empty: + continue + geom_area_m2 = polygon_area_m2(geom, center_lat) + if geom_area_m2 < min_area_m2: + continue + minx, miny, maxx, maxy = geom.bounds + col0 = max(0, int(math.floor((minx - info.x0) / info.dx)) - 1) + col1 = min(info.samples, int(math.ceil((maxx - info.x0) / info.dx)) + 1) + row0 = max(0, int(math.floor((info.y0 - maxy) / info.dy)) - 1) + row1 = min(info.lines, int(math.ceil((info.y0 - miny) / info.dy)) + 1) + if col1 <= col0 or row1 <= row0: + continue + window = cartographic_water[row0:row1, col0:col1] + pixels = int(window.sum()) + score_window = score[row0:row1, col0:col1] + hh_window = hh_db[row0:row1, col0:col1] + hv_window = hv_db[row0:row1, col0:col1] + slope_window = slope[row0:row1, col0:col1] if slope is not None else None + mean_slope = float(np.nanmean(slope_window[window])) if slope_window is not None and np.any(np.isfinite(slope_window[window])) else -9999.0 + feature = { + "geometry": geom.__geo_interface__, + "properties": { + "class_id": CLASS_CARTOGRAPHIC_WATER, + "class_name": CLASS_NAMES[CLASS_CARTOGRAPHIC_WATER], + "confidence": "map", + "area_m2": float(geom_area_m2), + "pixels": pixels, + "mean_score": float(np.nanmean(score_window[window])), + "mean_hh_db": float(np.nanmean(hh_window[window])), + "mean_hv_db": float(np.nanmean(hv_window[window])), + "mean_slope": mean_slope, + "known_water": int(np.any(known_water[row0:row1, col0:col1] & window)), + "paddy": int(np.any(paddy[row0:row1, col0:col1] & window)), + "review_flag": 1, + }, + } + if CLASS_CARTOGRAPHIC_WATER in handles: + handles[CLASS_CARTOGRAPHIC_WATER].write(feature) + if CLASS_CARTOGRAPHIC_WATER in shp_handles: + shp_handles[CLASS_CARTOGRAPHIC_WATER].write(feature) + counts[CLASS_NAMES[CLASS_CARTOGRAPHIC_WATER]] += 1 + finally: + for handle in list(handles.values()) + list(shp_handles.values()): + handle.close() + return { + "gpkg_path": str(gpkg_path) if gpkg_path is not None else None, + "shp_dir": str(shp_dir) if shp_dir is not None else None, + "layers": layers, + "feature_counts": counts, + "min_area_m2": float(min_area_m2), + "simplify_meters": float(simplify_meters), + "smooth_meters": float(smooth_meters), + "min_hole_area_m2": float(min_hole_area_m2), + } diff --git a/backend/app/routers/flood.py b/backend/app/routers/flood.py index 97ebdc2..6490f38 100644 --- a/backend/app/routers/flood.py +++ b/backend/app/routers/flood.py @@ -24,6 +24,10 @@ class FloodPreprocessRequest(BaseModel): class FloodWaterExtractionRequest(BaseModel): scene_id: Optional[int] = Field(default=None, description="SARSceneGeoORM primary key") input_path: Optional[str] = Field(default=None, description="Direct GeoTIFF/ENVI input path") + processor: Optional[str] = Field(default=None, description="Water extraction processor code") + hh_path: Optional[str] = Field(default=None, description="GF-3 HH geocoded ENVI path") + hv_path: Optional[str] = Field(default=None, description="GF-3 HV geocoded ENVI path") + processor_params: Optional[dict[str, Any]] = Field(default=None, description="Processor-specific parameters") class FloodPairSearchRequest(BaseModel): diff --git a/backend/app/services/flood_analysis_service.py b/backend/app/services/flood_analysis_service.py index 99f1291..0e0594a 100644 --- a/backend/app/services/flood_analysis_service.py +++ b/backend/app/services/flood_analysis_service.py @@ -134,6 +134,72 @@ def _scene_analysis_path(scene: SARSceneGeoORM | None) -> str | None: return scene.analysis_tif_path +def _normalize_processor(value: Any) -> str: + processor = str(value or "").strip().lower() + if processor in {"gf3_water", "gf3_water_hh_hv", "gf3_hh_hv", "hh_hv"}: + return "gf3_hh_hv" + return processor or "otsu" + + +def _metadata_dict(value: Any) -> dict[str, Any]: + if isinstance(value, dict): + return value + if isinstance(value, str) and value.strip(): + try: + parsed = json.loads(value) + return parsed if isinstance(parsed, dict) else {} + except Exception: + return {} + return {} + + +def _find_standard_asset_by_pol(metadata: dict[str, Any], polarization: str) -> str | None: + target = str(polarization or "").upper() + candidates = [] + candidates.extend(metadata.get("standard_assets") or []) + nested = metadata.get("metadata") + if isinstance(nested, dict): + candidates.extend(nested.get("standard_assets") or []) + for asset in candidates: + if not isinstance(asset, dict): + continue + if str(asset.get("polarization") or "").upper() != target: + continue + for key in ("source_native", "path"): + path = str(asset.get(key) or "").strip() + if path: + return path + return None + + +def _resolve_gf3_hh_hv_inputs( + *, + req: Any, + scene: SARSceneGeoORM | None, + radar: RadarDataORM | None, +) -> tuple[str | None, str | None, dict[str, Any]]: + hh_path = str(getattr(req, "hh_path", "") or "").strip() or None + hv_path = str(getattr(req, "hv_path", "") or "").strip() or None + resolution: dict[str, Any] = {"source": "request" if hh_path or hv_path else "metadata"} + if hh_path and hv_path: + return hh_path, hv_path, resolution + + scene_meta = _metadata_dict(scene.analysis_metadata_json if scene else None) + radar_meta = _metadata_dict(radar.metadata_json if radar else None) + for metadata in (scene_meta, radar_meta): + hh_path = hh_path or _find_standard_asset_by_pol(metadata, "HH") + hv_path = hv_path or _find_standard_asset_by_pol(metadata, "HV") + resolution.update( + { + "scene_metadata_used": bool(scene_meta), + "radar_metadata_used": bool(radar_meta), + "hh_auto_resolved": bool(hh_path), + "hv_auto_resolved": bool(hv_path), + } + ) + return hh_path, hv_path, resolution + + def _resolve_aoi_wkt_from_request(req: Any) -> tuple[str, dict[str, Any], dict[str, Any]]: """Resolve region/GeoJSON AOI using the same parser as the management search page.""" aoi_geojson = getattr(req, "aoi_geojson", None) @@ -393,22 +459,44 @@ async def list_scenes(limit: int, offset: int, db: AsyncSession) -> dict[str, An async def submit_water_extraction(req: Any, db: AsyncSession) -> dict[str, Any]: input_path = req.input_path scene_id = req.scene_id + processor = _normalize_processor(getattr(req, "processor", None)) + processor_params = dict(getattr(req, "processor_params", None) or {}) + scene: SARSceneGeoORM | None = None + radar: RadarDataORM | None = None if scene_id: scene = await db.get(SARSceneGeoORM, scene_id) if not scene: raise HTTPException(status_code=404, detail=f"SARSceneGeoORM id={scene_id} not found") - input_path = _scene_analysis_path(scene) - if not input_path: + if scene.radar_data_id: + radar = await db.get(RadarDataORM, scene.radar_data_id) + if processor == "gf3_hh_hv": + input_path = input_path or _scene_analysis_path(scene) + else: + input_path = _scene_analysis_path(scene) + if processor != "gf3_hh_hv" and not input_path: raise HTTPException(status_code=400, detail="Scene has no analysis-ready GeoTIFF") - if not input_path: + hh_path = None + hv_path = None + input_resolution: dict[str, Any] = {} + if processor == "gf3_hh_hv": + hh_path, hv_path, input_resolution = _resolve_gf3_hh_hv_inputs(req=req, scene=scene, radar=radar) + if not hh_path or not hv_path: + raise HTTPException(status_code=400, detail="GF3 HH/HV water extraction requires hh_path and hv_path") + input_path = input_path or hh_path + elif not input_path: raise HTTPException(status_code=400, detail="scene_id or input_path is required") extraction = WaterExtractionORM( scene_id=scene_id, - processor=getattr(req, "processor", None) or "otsu", + processor=processor, input_path=input_path, + metadata_json={ + "processor_params": processor_params, + "input_assets": {"hh": hh_path, "hv": hv_path} if processor == "gf3_hh_hv" else {}, + "input_resolution": input_resolution, + }, status="PENDING", ) db.add(extraction) @@ -421,7 +509,13 @@ async def submit_water_extraction(req: Any, db: AsyncSession) -> dict[str, Any]: job_type=JOB_TYPE_WATER_DETECT, task_type=f"FLOOD_WATER_EXTRACTION_{extraction_id}", task_name=f"Flood water extraction id={extraction_id}", - payload={"extraction_id": extraction_id, "processor": extraction.processor}, + payload={ + "extraction_id": extraction_id, + "processor": extraction.processor, + "hh_path": hh_path, + "hv_path": hv_path, + "processor_params": processor_params, + }, ) async with db.begin(): queued_extraction = await db.get(WaterExtractionORM, extraction_id) diff --git a/backend/app/services/gf3_water_extraction_service.py b/backend/app/services/gf3_water_extraction_service.py new file mode 100644 index 0000000..c542795 --- /dev/null +++ b/backend/app/services/gf3_water_extraction_service.py @@ -0,0 +1,243 @@ +"""GF-3 HH/HV water extraction adapter for the flood-analysis job chain.""" +from __future__ import annotations + +import json +import math +import os +from pathlib import Path +from typing import Any + +from ..config import settings +from ..processors.gf3_water import WaterExtractionConfig, run_water_extraction + + +GF3_HH_HV_PROCESSOR = "gf3_hh_hv" + + +def _existing_path(value: str | os.PathLike[str] | None, *, label: str) -> Path: + if not value: + raise ValueError(f"{label} is required") + path = Path(str(value)) + if not path.exists(): + raise FileNotFoundError(f"{label} does not exist: {path}") + return path + + +def _optional_existing_path(value: str | os.PathLike[str] | None, *, label: str) -> Path | None: + if not value: + return None + path = Path(str(value)) + if not path.exists(): + raise FileNotFoundError(f"{label} does not exist: {path}") + return path + + +def _as_bool(value: Any, default: bool) -> bool: + if value is None: + return default + if isinstance(value, bool): + return value + text = str(value).strip().lower() + if text in {"1", "true", "yes", "y", "on"}: + return True + if text in {"0", "false", "no", "n", "off"}: + return False + return default + + +def _numeric_param(params: dict[str, Any], name: str, default: Any, cast: type) -> Any: + value = params.get(name, default) + if value is None or value == "": + return default + try: + return cast(value) + except (TypeError, ValueError): + return default + + +def _path_param_list(params: dict[str, Any], name: str) -> list[Path]: + raw = params.get(name) or [] + if isinstance(raw, (str, os.PathLike)): + raw = [raw] + paths: list[Path] = [] + for item in raw: + if not item: + continue + paths.append(_existing_path(item, label=name)) + return paths + + +def _first_existing(*paths: Path) -> str | None: + for path in paths: + if path.exists(): + return str(path) + return None + + +def _read_metadata(path: Path) -> dict[str, Any]: + try: + return json.loads(path.read_text(encoding="utf-8")) + except FileNotFoundError: + return {} + except json.JSONDecodeError as exc: + raise ValueError(f"GF3 water metadata is not valid JSON: {path}") from exc + + +def _pixel_area_km2(transform: Any, crs: Any, bounds: Any) -> float: + px_w = abs(float(transform.a)) + px_h = abs(float(transform.e)) + if crs and getattr(crs, "is_geographic", False): + lat_center = (float(bounds.top) + float(bounds.bottom)) / 2.0 + px_w_m = px_w * math.cos(math.radians(lat_center)) * 111_320.0 + px_h_m = px_h * 111_320.0 + else: + px_w_m, px_h_m = px_w, px_h + return max(0.0, (px_w_m * px_h_m) / 1_000_000.0) + + +def _water_area_km2(mask_path: str | None, *, water_pixel_count: int | None = None) -> float | None: + if not mask_path or not os.path.isfile(mask_path): + return None + try: + import rasterio + + with rasterio.open(mask_path) as src: + data = src.read(1) + pixel_count = int((data > 0).sum()) if water_pixel_count is None else int(water_pixel_count) + return round(pixel_count * _pixel_area_km2(src.transform, src.crs, src.bounds), 4) + except Exception: + return None + + +def _vector_runtime_available() -> bool: + try: + import fiona # noqa: F401 + import rasterio.features # noqa: F401 + + return True + except Exception: + return False + + +def run_gf3_hh_hv_water_extraction( + *, + hh_path: str, + hv_path: str, + output_dir: str, + job_id: str | None = None, + params: dict[str, Any] | None = None, +) -> dict[str, Any]: + """Run the embedded GF-3 HH/HV water extractor and normalize outputs.""" + params = dict(params or {}) + out_dir = Path(output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + hh = _existing_path(hh_path, label="HH input") + hv = _existing_path(hv_path, label="HV input") + dltb_cache_dir = _optional_existing_path( + params.get("dltb_cache_dir") or settings.GF3_WATER_DLTB_CACHE_DIR, + label="GF3_WATER_DLTB_CACHE_DIR", + ) + dem = _optional_existing_path( + params.get("dem") or params.get("dem_path") or settings.GF3_WATER_DEM_PATH, + label="GF3_WATER_DEM_PATH", + ) + + cartographic_water = _as_bool( + params.get("cartographic_water"), + bool(settings.GF3_WATER_DEFAULT_CARTOGRAPHIC), + ) + out_vector = _as_bool(params.get("out_vector"), bool(settings.GF3_WATER_DEFAULT_OUT_VECTOR)) + vector_runtime_available = _vector_runtime_available() + vector_output_enabled = out_vector and vector_runtime_available + + config = WaterExtractionConfig( + hh=hh, + hv=hv, + out_dir=out_dir, + dem=dem, + dltb_cache_dir=dltb_cache_dir, + dltb_mode=str(params.get("dltb_mode") or "soft"), + water_vector=_path_param_list(params, "water_vector"), + paddy_vector=_path_param_list(params, "paddy_vector"), + threshold_method=str(params.get("threshold_method") or "percentile"), + score_percentile=_numeric_param(params, "score_percentile", 95.0, float), + hv_percentile=_numeric_param(params, "hv_percentile", 20.0, float), + candidate_score_percentile=_numeric_param(params, "candidate_score_percentile", 90.0, float), + candidate_hv_percentile=_numeric_param(params, "candidate_hv_percentile", 50.0, float), + close_pixels=_numeric_param(params, "close_pixels", 2, int), + open_pixels=_numeric_param(params, "open_pixels", 1, int), + fill_hole_pixels=_numeric_param(params, "fill_hole_pixels", 2048, int), + min_component_pixels=_numeric_param(params, "min_component_pixels", 4096, int), + candidate_open_pixels=_numeric_param(params, "candidate_open_pixels", 1, int), + cartographic_water=cartographic_water, + cartographic_close_pixels=_numeric_param(params, "cartographic_close_pixels", 4, int), + cartographic_fill_hole_pixels=_numeric_param(params, "cartographic_fill_hole_pixels", 20000, int), + cartographic_min_component_pixels=_numeric_param(params, "cartographic_min_component_pixels", 4096, int), + out_vector_shp_dir=out_dir / "shp" if vector_output_enabled else None, + min_polygon_area_m2=_numeric_param(params, "min_polygon_area_m2", 50000.0, float), + simplify_meters=_numeric_param(params, "simplify_meters", 3.0, float), + smooth_meters=_numeric_param(params, "smooth_meters", 5.0, float), + min_hole_area_m2=_numeric_param(params, "min_hole_area_m2", 5000.0, float), + ) + + exit_code = run_water_extraction(config) + metadata_path = out_dir / "metadata.json" + metadata = _read_metadata(metadata_path) + if exit_code != 0: + return { + "ok": False, + "processor": GF3_HH_HV_PROCESSOR, + "error": f"GF3 HH/HV water extraction failed with exit code {exit_code}", + "metadata_json": metadata, + } + + output_path = _first_existing( + out_dir / "cartographic_water.tif", + out_dir / "water_mask.tif", + out_dir / "classified_water.tif", + ) + preview_path = _first_existing(out_dir / "preview_overlay.png", out_dir / "classified_preview.png") + vector_path = _first_existing(out_dir / "shp" / "cartographic_water.shp", out_dir / "water_products.gpkg") + water_pixels = ( + metadata.get("cartographic_water_pixels") + if metadata.get("cartographic_water_enabled") + else metadata.get("water_pixels") + ) + if water_pixels is None: + water_pixels = metadata.get("water_pixels") + water_pixel_count = int(water_pixels or 0) + area_km2 = _water_area_km2(output_path, water_pixel_count=water_pixel_count) + + normalized_metadata = { + **metadata, + "processor": GF3_HH_HV_PROCESSOR, + "job_id": job_id, + "metadata_path": str(metadata_path), + "output_path": output_path, + "preview_path": preview_path, + "vector_path": vector_path, + "input_assets": { + "hh": str(hh), + "hv": str(hv), + "dem": str(dem) if dem else None, + "dltb_cache_dir": str(dltb_cache_dir) if dltb_cache_dir else None, + }, + "runtime": { + "vector_requested": out_vector, + "vector_runtime_available": vector_runtime_available, + "vector_output_enabled": vector_output_enabled, + }, + } + + return { + "ok": True, + "processor": GF3_HH_HV_PROCESSOR, + "output_path": output_path, + "preview_path": preview_path, + "vector_path": vector_path, + "water_area_km2": area_km2, + "water_pixel_count": water_pixel_count, + "threshold_value": metadata.get("score_threshold"), + "metadata_json": normalized_metadata, + } diff --git a/backend/app/services/job_handlers.py b/backend/app/services/job_handlers.py index 7d7764f..3c4b481 100644 --- a/backend/app/services/job_handlers.py +++ b/backend/app/services/job_handlers.py @@ -3490,6 +3490,7 @@ async def _handle_flood_detection(job: SystemJobORM) -> None: async def _handle_water_detect(job: SystemJobORM) -> None: """水体检测 job handler（Otsu + DEM + 形态学 + 连通分量）。""" + from .gf3_water_extraction_service import GF3_HH_HV_PROCESSOR, run_gf3_hh_hv_water_extraction from .water_extraction_service import run_otsu_water_extraction payload = job.payload or {} @@ -3508,6 +3509,15 @@ async def _handle_water_detect(job: SystemJobORM) -> None: model_name = "WaterExtractionORM" if use_extraction_table else "WaterDetectionORM" raise ValueError(f"{model_name} id={record_id} 不存在") input_path = det.input_path + metadata = det.metadata_json if use_extraction_table and isinstance(det.metadata_json, dict) else {} + processor = str(payload.get("processor") or getattr(det, "processor", None) or "otsu").strip().lower() + if processor in {"gf3_water", "gf3_water_hh_hv", "hh_hv"}: + processor = GF3_HH_HV_PROCESSOR + input_assets = metadata.get("input_assets") if isinstance(metadata.get("input_assets"), dict) else {} + processor_params = dict(metadata.get("processor_params") or {}) + processor_params.update(dict(payload.get("processor_params") or {})) + hh_path = payload.get("hh_path") or input_assets.get("hh") + hv_path = payload.get("hv_path") or input_assets.get("hv") det.status = "RUNNING" if use_extraction_table and hasattr(det, "task_id"): det.task_id = job.task_id @@ -3518,17 +3528,28 @@ async def _handle_water_detect(job: SystemJobORM) -> None: mirror.task_id = job.task_id await db.commit() - if not input_path: + if processor == GF3_HH_HV_PROCESSOR: + if not hh_path or not hv_path: + raise ValueError("GF3 HH/HV water extraction requires hh_path and hv_path") + elif not input_path: raise ValueError("水体检测缺少输入路径 input_path") output_name = f"water_extraction_{record_id}" if use_extraction_table else f"water_detect_{record_id}" output_root = settings.WATER_RESULTS_DIR or os.path.join(settings.BACKEND_DIR, "water_results") - output_dir = os.path.join(output_root, output_name) + output_dir = os.path.join(output_root, processor, output_name) if use_extraction_table else os.path.join(output_root, output_name) os.makedirs(output_dir, exist_ok=True) await task_service.update_task(job.task_id, progress=10, message="启动水体检测算法...") def _run() -> Dict[str, Any]: + if processor == GF3_HH_HV_PROCESSOR: + return run_gf3_hh_hv_water_extraction( + hh_path=hh_path, + hv_path=hv_path, + output_dir=output_dir, + job_id=job.job_id, + params=processor_params, + ) return run_otsu_water_extraction( input_path=input_path, output_dir=output_dir, @@ -3557,16 +3578,24 @@ async def _handle_water_detect(job: SystemJobORM) -> None: if det: if result.get("ok"): det.output_path = result.get("output_path") + if hasattr(det, "preview_path"): + det.preview_path = result.get("preview_path") + if hasattr(det, "vector_path"): + det.vector_path = result.get("vector_path") det.water_area_km2 = result.get("water_area_km2") det.water_pixel_count = result.get("water_pixel_count") if use_extraction_table: det.processor = result.get("processor") or det.processor or "otsu" det.threshold_value = result.get("threshold_value") - det.metadata_json = { - "legacy_otsu_threshold_db": result.get("otsu_threshold_db"), - "value_transform": result.get("value_transform"), - "job_id": job.job_id, - } + result_metadata = result.get("metadata_json") + if isinstance(result_metadata, dict): + det.metadata_json = result_metadata + else: + det.metadata_json = { + "legacy_otsu_threshold_db": result.get("otsu_threshold_db"), + "value_transform": result.get("value_transform"), + "job_id": job.job_id, + } else: det.otsu_threshold_db = result.get("otsu_threshold_db") det.status = "DONE" @@ -3578,6 +3607,8 @@ async def _handle_water_detect(job: SystemJobORM) -> None: mirror = await db.get(WaterExtractionORM, int(record_id)) if mirror: mirror.output_path = det.output_path + mirror.preview_path = result.get("preview_path") + mirror.vector_path = result.get("vector_path") mirror.water_area_km2 = det.water_area_km2 mirror.water_pixel_count = det.water_pixel_count mirror.threshold_value = result.get("threshold_value") or result.get("otsu_threshold_db") diff --git a/docs/GF3_WATER_EXTRACTION_INTEGRATION_20260615.md b/docs/GF3_WATER_EXTRACTION_INTEGRATION_20260615.md new file mode 100644 index 0000000..2987fd7 --- /dev/null +++ b/docs/GF3_WATER_EXTRACTION_INTEGRATION_20260615.md @@ -0,0 +1,157 @@ +# GF3 Water Extraction Integration Design + +Date: 2026-06-15 + +## Scope + +This document describes how to embed the GF-3 HH/HV water extraction work from `D:\Code\Water` into the management system. + +The goal is integration, not bulk import. Code that is reusable should be copied into the backend; code that is workflow-specific should be rewritten around the existing job system; large data and generated products should be transferred or registered as runtime assets outside Git. + +## Source Project Inventory + +`D:\Code\Water` contains: + +- `gf3_water/`: active Python package for GF-3 HH/HV water extraction. +- `scripts/`: thin CLI/data-preparation wrappers. +- `docs/`: algorithm and integration notes. +- `data/`: local scenes and prior data. +- `outputs/`: generated products. +- `gf3-water-ai4g-unet/`: legacy deep-learning experiment/checkpoints. + +The active production-facing implementation is the non-DL `gf3_water` package. It consumes SARscape ENVI HH/HV geocoded assets, for example `_hh_geo` and `_hv_geo`, and writes raster, preview, vector and `metadata.json` outputs. + +## Current System Entry Points + +The existing system already has a suitable production chain: + +1. `/flood/water-extractions` creates a `WaterExtractionORM` record. +2. The record is submitted as a `WATER_DETECT` job. +3. `backend/app/services/job_handlers.py::_handle_water_detect` runs the processor. +4. Results are written back to `WaterExtractionORM`. +5. `FloodAnalysisWorkspace.jsx` lists and previews extraction results. + +The GF-3 HH/HV processor should be embedded into this chain instead of creating a parallel table or router. + +## Copy, Rewrite, Transfer + +### Copy Into The Backend + +Copy the reusable algorithm package: + +- From: `D:\Code\Water\gf3_water` +- To: `backend/app/processors/gf3_water` + +The copied package should remain close to the original algorithm code so it can be compared and upgraded later. System-specific behavior should live in a separate service wrapper. + +Optional later copy: + +- `D:\Code\Water\scripts\water_baseline_hh_hv.py` only if a local CLI smoke-test entry is needed. + +Do not copy: + +- `outputs/` +- `data/scenes/` +- `data/raw/` +- `gf3-water-ai4g-unet/*.pt` + +### Rewrite In The System + +System integration should be written around existing services: + +- Add a backend wrapper service, for example `backend/app/services/gf3_water_extraction_service.py`. +- Extend `/flood/water-extractions` request schema with: + - `processor` + - `hh_path` + - `hv_path` + - `processor_params` +- Add HH/HV asset resolution from `SARSceneGeoORM.analysis_metadata_json.standard_assets` or the related `RadarDataORM.metadata_json.standard_assets`. +- Dispatch `WATER_DETECT` by `processor`: + - `otsu`: keep existing single-raster processor. + - `gf3_hh_hv`: run the embedded GF-3 HH/HV package. +- Map GF-3 outputs into `WaterExtractionORM`: + - `output_path`: prefer `cartographic_water.tif`, fallback `water_mask.tif`. + - `preview_path`: prefer `preview_overlay.png`, fallback `classified_preview.png`. + - `vector_path`: prefer `shp/cartographic_water.shp`, fallback `water_products.gpkg`. + - `water_pixel_count`: prefer `cartographic_water_pixels` when cartographic output is enabled. + - `threshold_value`: `score_threshold`. + - `metadata_json`: original metadata plus resolved input/output paths and processor version. +- Update frontend water extraction UI: + - Processor selector: fast Otsu vs GF-3 HH/HV. + - For GF-3 HH/HV, show whether HH/HV assets are auto-resolved for the selected scene. + - Allow manual HH/HV paths for early operations and debugging. + - Display processor/output type in result rows. + +### Transfer Or Register As Runtime Assets + +Do not store runtime data in Git. + +Required/valuable runtime assets: + +- `D:\Code\Water\data\priors\dltb_cache\heilongjiang` + - Current size observed: 25 files, about 8.5 GB. + - Transfer to a managed runtime asset path, for example `D:\production_assets\gf3_water\priors\dltb_cache\heilongjiang`. + - Configure by `GF3_WATER_DLTB_CACHE_DIR`. + +Optional runtime assets: + +- Hydro prior vectors under `data/priors/hydro`. +- DEM path if slope filtering should be enabled. + +Build-only assets: + +- `data/raw/dltb/DLTB_2025.gdb` +- Heilongjiang boundary Shapefile used by cache-build scripts. + +These are only needed to rebuild the DLTB cache and should be kept in external storage. + +Legacy DL assets: + +- `gf3-water-ai4g-unet/best.pt` +- `gf3-water-ai4g-unet/last.pt` + +These are large model artifacts and should stay outside the application until a separate model registry/runtime is designed. + +## Recommended Runtime Configuration + +Add configuration keys: + +- `GF3_WATER_DLTB_CACHE_DIR` +- `GF3_WATER_DEM_PATH` +- `GF3_WATER_DEFAULT_CARTOGRAPHIC=true` +- `GF3_WATER_DEFAULT_OUT_VECTOR=true` + +`WATER_RESULTS_DIR` remains the output root for generated extraction products. + +## Production Workflow + +```mermaid +flowchart TD + A["GF3 SARscape native production"] --> B["GF3 standardization manifest"] + B --> C["SARSceneGeoORM DONE"] + C --> D["User submits water extraction"] + D --> E{"processor"} + E -->|otsu| F["Existing single-raster Otsu"] + E -->|gf3_hh_hv| G["Resolve HH/HV assets"] + G --> H["Run embedded gf3_water package"] + F --> I["Update WaterExtractionORM"] + H --> I + I --> J["Preview, list, publish through flood UI"] +``` + +## Implementation Order + +1. Copy `gf3_water` package into `backend/app/processors/gf3_water`. +2. Add GF-3 HH/HV wrapper service and output mapping. +3. Add config keys and `.env.example` entries. +4. Extend request schema and submission metadata. +5. Add processor dispatch in `WATER_DETECT`. +6. Add HH/HV auto-resolution from standard GF3 assets. +7. Update frontend extraction controls/result display. +8. Verify with backend compile/import checks and frontend build. + +## Open Decisions + +- Whether runtime execution should be in-process Python API first or always subprocess CLI. Initial integration should use in-process API because it fits the existing worker model and keeps job accounting simple. +- Whether to transfer the 8.5 GB DLTB cache automatically. This should be a deployment operation, not a Git operation. +- Whether legacy AI4G U-Net should be supported later. It should not block the current HH/HV production chain. diff --git a/frontend/src/FloodAnalysisWorkspace.jsx b/frontend/src/FloodAnalysisWorkspace.jsx index c974406..b201198 100644 --- a/frontend/src/FloodAnalysisWorkspace.jsx +++ b/frontend/src/FloodAnalysisWorkspace.jsx @@ -316,6 +316,7 @@ function SceneRow({ scene, readOnly, onShowMap, onExtractWater, onReset }) { } function WaterResultRow({ item, onShowMap }) { + const thresholdValue = item.threshold_value ?? item.otsu_threshold_db; return (
@@ -324,6 +325,7 @@ function WaterResultRow({ item, onShowMap }) { 水体提取 #{item.id} {item.scene_id && 场景 #{item.scene_id}} + {item.processor && {item.processor}} {item.satellite && {item.satellite}} {item.imaging_date && {formatYmd(item.imaging_date, 'zh')}} {item.polarization && {item.polarization}} @@ -331,8 +333,13 @@ function WaterResultRow({ item, onShowMap }) {
- +
+ {(item.preview_path || item.vector_path) && ( +
+ {[item.preview_path ? `preview: ${item.preview_path}` : null, item.vector_path ? `vector: ${item.vector_path}` : null].filter(Boolean).join(' | ')} +
+ )} {item.error_msg &&
{item.error_msg}
}
@@ -451,6 +458,7 @@ export default function FloodAnalysisWorkspace({ const [waterTotal, setWaterTotal] = useState(0); const [waterPage, setWaterPage] = useState(0); const [waterLoading, setWaterLoading] = useState(false); + const [waterProcessor, setWaterProcessor] = useState('otsu'); const [sourceAoiMode, setSourceAoiMode] = useState('none'); const [sourceRegionOptions, setSourceRegionOptions] = useState({ provinces: [], cities: [] }); @@ -774,7 +782,7 @@ export default function FloodAnalysisWorkspace({ showMessage('info', `正在提交场景 #${scene.id} 的水体提取任务...`); onTaskStart?.(null, '正在提交水体提取任务...'); try { - const res = await submitFloodWaterExtraction({ scene_id: scene.id }); + const res = await submitFloodWaterExtraction({ scene_id: scene.id, processor: waterProcessor }); if (res.data?.task_id) onTaskStart?.(res.data.task_id, '水体提取任务已启动'); await loadWaterResults(0); showMessage('success', `场景 #${scene.id} 的水体提取任务已提交。`); @@ -1181,6 +1189,13 @@ export default function FloodAnalysisWorkspace({
loadScenes(scenesPage)}>刷新场景} /> +
+ 水体提取处理器 + +
{scenesLoading && 场景加载中...} {!scenesLoading && scenes.length === 0 && 暂无分析就绪场景。}