Compare commits

...
5 Commits
32 changed files with 3121 additions and 43 deletions
+10 -4
View File
@@ -59,15 +59,18 @@ ALLOWED_EXPORT_DIRS=
# 源数据目录 # 源数据目录
# ----------------------------------------------------------------------------- # -----------------------------------------------------------------------------
UNPACK_SOURCE_DIRS=D:\Archives UNPACK_SOURCE_DIRS=D:\Archives
SOURCE_PRODUCT_DIRS=D:\LuTan1_Image_Pool;D:\Sentinel1_Image_Pool_ZIP TASK_POOL_ROOT=D:\Task_Pool
DINSAR_TASK_POOL_ROOT=D:\Task_Pool\DInSAR
SBAS_TASK_POOL_ROOT=D:\Task_Pool\SBAS
SOURCE_PRODUCT_DIRS=D:\LuTan1_Image_Pool;D:\Sentinel1_Image_Pool_ZIP;\\DESKTOP-N16HJ84\InSAR_Storage_2\LuTan-1\Archive;\\DESKTOP-N16HJ84\InSAR_Storage_2\Sentinel-1\Archive
SENTINEL1_STORAGE_DIRS=D:\Sentinel1_Image_Pool SENTINEL1_STORAGE_DIRS=D:\Sentinel1_Image_Pool
INSAR_STORAGE_DIRS=D:\LuTan1_Image_Pool INSAR_STORAGE_DIRS=D:\LuTan1_Image_Pool
MONITOR_RADAR_DIRS=D:\LuTan1_Image_Pool MONITOR_RADAR_DIRS=D:\LuTan1_Image_Pool
MONITOR_DINSAR_DIRS=D:\DInSARResult MONITOR_DINSAR_DIRS=D:\DInSARResult
ORBIT_SOURCE_DIRS=D:\LT1_data_lsarorbit;D:\Sentinel1_EOF_Pool ORBIT_SOURCE_DIRS=D:\LT1_data_lsarorbit;D:\Sentinel1_EOF_Pool;\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\LuTan-1;\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\Sentinel-1
MONITOR_ORBIT_DIR=D:\LT1_data_lsarorbit MONITOR_ORBIT_DIR=D:\LT1_data_lsarorbit
GF3_ARCHIVE_SOURCE_DIRS=D:\production_inputs\gf3\archives GF3_ARCHIVE_SOURCE_DIRS=\\DESKTOP-N16HJ84\InSAR_Storage_1\GaoFen-3
GF3_ARCHIVE_EXTS=.zip,.tar,.tar.gz,.tgz GF3_ARCHIVE_EXTS=.zip,.tar,.tar.gz,.tgz
GF3_UNPACK_DELETE_ARCHIVE=true GF3_UNPACK_DELETE_ARCHIVE=true
GF3_LEGACY_GDAL_ENABLED=false GF3_LEGACY_GDAL_ENABLED=false
@@ -119,6 +122,9 @@ IDL_DINSAR_DEM_BASE_FILE=D:\SRTM30m\SRTMDEM_RSP_SARscape
SRTM_DEM_DIR=D:\SRTM30m SRTM_DEM_DIR=D:\SRTM30m
GF3_GEO_DEM_PATH=D:\DEM\gf3_dem.jp2 GF3_GEO_DEM_PATH=D:\DEM\gf3_dem.jp2
WATER_RESULTS_DIR=D:\WaterResult WATER_RESULTS_DIR=D:\WaterResult
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_READY_ROOT=D:\production_results\sar_analysis_ready
SAR_ANALYSIS_WORK_ROOT=D:\production_runtime\sar_analysis_work SAR_ANALYSIS_WORK_ROOT=D:\production_runtime\sar_analysis_work
SAR_ANALYSIS_NODATA_VALUE=-9999 SAR_ANALYSIS_NODATA_VALUE=-9999
@@ -272,7 +278,7 @@ GAMMA_SBAS_RUNTIME_ID=gamma_sbas_runtime_v1
GAMMA_SBAS_WSL_DISTRO=Ubuntu-24.04 GAMMA_SBAS_WSL_DISTRO=Ubuntu-24.04
GAMMA_SBAS_PYTHON=/home/administrator/miniconda3/envs/insar_wsl_v1/bin/python GAMMA_SBAS_PYTHON=/home/administrator/miniconda3/envs/insar_wsl_v1/bin/python
GAMMA_SBAS_ENV_SCRIPT=D:\Code\Insar_management_system_v2\deploy\wsl\profiles\gamma_env.sh GAMMA_SBAS_ENV_SCRIPT=D:\Code\Insar_management_system_v2\deploy\wsl\profiles\gamma_env.sh
GAMMA_SBAS_WORK_ROOT=D:\production_runtime\sbas_insar_work GAMMA_SBAS_WORK_ROOT=D:\Task_Pool\SBAS
GAMMA_SBAS_PRODUCT_ROOT=D:\production_results\timeseries\sbas GAMMA_SBAS_PRODUCT_ROOT=D:\production_results\timeseries\sbas
GAMMA_SBAS_TRIAL_ROOT=D:\production_runtime\gamma_ipta_trials GAMMA_SBAS_TRIAL_ROOT=D:\production_runtime\gamma_ipta_trials
GAMMA_SBAS_SCRIPT_TEMPLATE_ROOT=D:\Code\Insar_management_system_v2\backend\templates\gamma_sbas GAMMA_SBAS_SCRIPT_TEMPLATE_ROOT=D:\Code\Insar_management_system_v2\backend\templates\gamma_sbas
+24 -1
View File
@@ -93,6 +93,14 @@ def _default_runtime_root(project_root: str) -> str:
return os.path.join(normalized_root, "runtime") return os.path.join(normalized_root, "runtime")
def _default_task_pool_root(project_root: str) -> str:
normalized_root = os.path.normpath(project_root)
drive, _tail = os.path.splitdrive(normalized_root)
if drive:
return os.path.join(drive + os.sep, "Task_Pool")
return os.path.join(normalized_root, "Task_Pool")
def _default_runtime_dir(project_root: str, *parts: str) -> str: def _default_runtime_dir(project_root: str, *parts: str) -> str:
return os.path.join(_default_runtime_root(project_root), *parts) return os.path.join(_default_runtime_root(project_root), *parts)
@@ -187,6 +195,9 @@ class Settings(BaseSettings):
DB_SCHEMA_RESET_CONFIRM: bool = False DB_SCHEMA_RESET_CONFIRM: bool = False
UNPACK_SOURCE_DIRS: str = "" UNPACK_SOURCE_DIRS: str = ""
TASK_POOL_ROOT: str = ""
DINSAR_TASK_POOL_ROOT: str = ""
SBAS_TASK_POOL_ROOT: str = ""
SOURCE_PRODUCT_DIRS: str = "" SOURCE_PRODUCT_DIRS: str = ""
SENTINEL1_STORAGE_DIRS: str = "" SENTINEL1_STORAGE_DIRS: str = ""
ORBIT_SOURCE_DIRS: str = "" ORBIT_SOURCE_DIRS: str = ""
@@ -215,6 +226,9 @@ class Settings(BaseSettings):
RADAR_PREVIEW_BUILD_ON_DEMAND: bool = True RADAR_PREVIEW_BUILD_ON_DEMAND: bool = True
WATER_RESULTS_DIR: str = "" WATER_RESULTS_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_READY_ROOT: str = ""
SAR_ANALYSIS_WORK_ROOT: str = "" SAR_ANALYSIS_WORK_ROOT: str = ""
SAR_ANALYSIS_NODATA_VALUE: float = -9999.0 SAR_ANALYSIS_NODATA_VALUE: float = -9999.0
@@ -438,6 +452,12 @@ class Settings(BaseSettings):
"WATER_RESULTS_DIR", "WATER_RESULTS_DIR",
os.path.join(backend_dir, "water_results"), os.path.join(backend_dir, "water_results"),
) )
if not self.TASK_POOL_ROOT:
object.__setattr__(self, "TASK_POOL_ROOT", _default_task_pool_root(project_root))
if not self.DINSAR_TASK_POOL_ROOT:
object.__setattr__(self, "DINSAR_TASK_POOL_ROOT", os.path.join(self.TASK_POOL_ROOT, "DInSAR"))
if not self.SBAS_TASK_POOL_ROOT:
object.__setattr__(self, "SBAS_TASK_POOL_ROOT", os.path.join(self.TASK_POOL_ROOT, "SBAS"))
if not self.SAR_ANALYSIS_READY_ROOT: if not self.SAR_ANALYSIS_READY_ROOT:
object.__setattr__( object.__setattr__(
self, self,
@@ -777,7 +797,7 @@ class Settings(BaseSettings):
object.__setattr__( object.__setattr__(
self, self,
"GAMMA_SBAS_WORK_ROOT", "GAMMA_SBAS_WORK_ROOT",
os.path.join(_default_runtime_root(project_root), "sbas_insar_work"), self.SBAS_TASK_POOL_ROOT,
) )
if not self.GAMMA_SBAS_PRODUCT_ROOT: if not self.GAMMA_SBAS_PRODUCT_ROOT:
object.__setattr__( object.__setattr__(
@@ -987,6 +1007,9 @@ class Settings(BaseSettings):
os.makedirs(settings.RESULT_QUARANTINE_ROOT, exist_ok=True) os.makedirs(settings.RESULT_QUARANTINE_ROOT, exist_ok=True)
os.makedirs(settings.SAR_ANALYSIS_READY_ROOT, exist_ok=True) os.makedirs(settings.SAR_ANALYSIS_READY_ROOT, exist_ok=True)
os.makedirs(settings.SAR_ANALYSIS_WORK_ROOT, exist_ok=True) os.makedirs(settings.SAR_ANALYSIS_WORK_ROOT, exist_ok=True)
os.makedirs(settings.TASK_POOL_ROOT, exist_ok=True)
os.makedirs(settings.DINSAR_TASK_POOL_ROOT, exist_ok=True)
os.makedirs(settings.SBAS_TASK_POOL_ROOT, exist_ok=True)
for path in split_env_paths(settings.GF3_ARCHIVE_SOURCE_DIRS): for path in split_env_paths(settings.GF3_ARCHIVE_SOURCE_DIRS):
os.makedirs(path, exist_ok=True) os.makedirs(path, exist_ok=True)
for path in split_env_paths(settings.GF3_SARSCAPE_NATIVE_DIRS): for path in split_env_paths(settings.GF3_SARSCAPE_NATIVE_DIRS):
+118 -2
View File
@@ -2,6 +2,7 @@
import shutil import shutil
import asyncio import asyncio
import tempfile import tempfile
import tarfile
import zipfile import zipfile
import json import json
import hashlib import hashlib
@@ -43,6 +44,100 @@ def find_dinsar_source_to_copy(path: str) -> str:
return path return path
_DINSAR_ARCHIVE_SUFFIXES = (".tar.gz", ".tgz", ".zip", ".tar")
def _is_supported_archive(path: str) -> bool:
lower = str(path or "").lower()
return any(lower.endswith(suffix) for suffix in _DINSAR_ARCHIVE_SUFFIXES)
def _safe_archive_member_name(member_name: str, archive_path: str) -> str:
name = str(member_name or "").replace("\\", "/").strip("/")
if not name or name.startswith("../") or "/../" in f"/{name}/":
raise ValueError(f"Unsafe archive member path in {archive_path}: {member_name}")
if os.path.isabs(name) or os.path.splitdrive(name)[0]:
raise ValueError(f"Unsafe archive member path in {archive_path}: {member_name}")
return name
def _extract_archive_to_dir(archive_path: str, dest_dir: str) -> int:
if zipfile.is_zipfile(archive_path):
extracted = 0
with zipfile.ZipFile(archive_path) as zip_obj:
for info in zip_obj.infolist():
rel_name = _safe_archive_member_name(info.filename, archive_path)
dest_path = os.path.abspath(os.path.join(dest_dir, rel_name))
if not dest_path.startswith(os.path.abspath(dest_dir) + os.sep):
raise ValueError(f"Unsafe ZIP member path: {info.filename}")
if info.is_dir():
os.makedirs(dest_path, exist_ok=True)
continue
os.makedirs(os.path.dirname(dest_path), exist_ok=True)
with zip_obj.open(info, "r") as source, open(dest_path, "wb") as target:
shutil.copyfileobj(source, target, length=1024 * 1024)
extracted += 1
return extracted
if tarfile.is_tarfile(archive_path):
extracted = 0
with tarfile.open(archive_path, "r:*") as tar_obj:
for member in tar_obj:
rel_name = _safe_archive_member_name(member.name, archive_path)
dest_path = os.path.abspath(os.path.join(dest_dir, rel_name))
if not dest_path.startswith(os.path.abspath(dest_dir) + os.sep):
raise ValueError(f"Unsafe TAR member path: {member.name}")
if member.isdir():
os.makedirs(dest_path, exist_ok=True)
continue
if not member.isfile():
continue
source = tar_obj.extractfile(member)
if source is None:
continue
os.makedirs(os.path.dirname(dest_path), exist_ok=True)
with source, open(dest_path, "wb") as target:
shutil.copyfileobj(source, target, length=1024 * 1024)
extracted += 1
return extracted
raise ValueError(f"Unsupported archive format: {archive_path}")
def _materialize_dinsar_source(source_path: str, dest_dir: str) -> Dict[str, Any]:
normalized = os.path.normpath(os.path.abspath(str(source_path or "")))
if not os.path.exists(normalized):
raise FileNotFoundError(normalized)
if os.path.isdir(normalized):
shutil.copytree(normalized, dest_dir, dirs_exist_ok=True)
return {"mode": "copy_directory", "source_path": normalized}
if os.path.isfile(normalized) and _is_supported_archive(normalized):
os.makedirs(dest_dir, exist_ok=True)
extracted = _extract_archive_to_dir(normalized, dest_dir)
if extracted <= 0:
raise OSError(f"Archive extraction produced no files: {normalized}")
return {
"mode": "extract_archive",
"source_path": normalized,
"archive_path": normalized,
"extracted_files": extracted,
}
if os.path.isfile(normalized):
os.makedirs(dest_dir, exist_ok=True)
target = os.path.join(dest_dir, os.path.basename(normalized))
shutil.copy2(normalized, target)
return {
"mode": "copy_file",
"source_path": normalized,
"relative_path": os.path.basename(target),
}
raise FileNotFoundError(normalized)
def _resolve_orbit_dest_path( def _resolve_orbit_dest_path(
orbit_dir: str, orbit_dir: str,
role: str, role: str,
@@ -191,6 +286,7 @@ def _build_dinsar_pair_metadata(
package_format: str, package_format: str,
include_orbit_files: bool, include_orbit_files: bool,
orbit_entries: List[Dict[str, Any]], orbit_entries: List[Dict[str, Any]],
source_materialization: Optional[Dict[str, Any]] = None,
) -> Dict[str, Any]: ) -> Dict[str, Any]:
return { return {
"pair_key": item.get("pair_key"), "pair_key": item.get("pair_key"),
@@ -214,6 +310,7 @@ def _build_dinsar_pair_metadata(
"master_orbit_file_path": item.get("master_orbit_file_path"), "master_orbit_file_path": item.get("master_orbit_file_path"),
"slave_orbit_file_path": item.get("slave_orbit_file_path"), "slave_orbit_file_path": item.get("slave_orbit_file_path"),
"orbit_files": orbit_entries, "orbit_files": orbit_entries,
"source_materialization": source_materialization or {},
"scene_pair_uid": item.get("scene_pair_uid") or item.get("pair_uid"), "scene_pair_uid": item.get("scene_pair_uid") or item.get("pair_uid"),
"pair_uid": item.get("pair_uid") or item.get("scene_pair_uid"), "pair_uid": item.get("pair_uid") or item.get("scene_pair_uid"),
"network_run_id": item.get("network_run_id"), "network_run_id": item.get("network_run_id"),
@@ -1094,8 +1191,26 @@ async def run_dinsar_copy_items(
failed_count += 1 failed_count += 1
continue continue
await asyncio.to_thread(shutil.copytree, master_src_path, master_dir, dirs_exist_ok=True) master_materialization = await asyncio.to_thread(
await asyncio.to_thread(shutil.copytree, slave_src_path, slave_dir, dirs_exist_ok=True) _materialize_dinsar_source,
master_src_path,
master_dir,
)
slave_materialization = await asyncio.to_thread(
_materialize_dinsar_source,
slave_src_path,
slave_dir,
)
source_materialization = {
"master": {
**master_materialization,
"target_relative_path": "master",
},
"slave": {
**slave_materialization,
"target_relative_path": "slave",
},
}
orbit_entries = await _copy_dinsar_orbit_files( orbit_entries = await _copy_dinsar_orbit_files(
task_id, task_id,
item, item,
@@ -1112,6 +1227,7 @@ async def run_dinsar_copy_items(
"zip" if export_zip else "folder", "zip" if export_zip else "folder",
include_orbit_files, include_orbit_files,
orbit_entries, orbit_entries,
source_materialization,
), ),
) )
if export_zip and zip_path: if export_zip and zip_path:
View File
@@ -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"]
+25
View File
@@ -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))
+64
View File
@@ -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))
@@ -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
@@ -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",
}
+146
View File
@@ -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,
)
+188
View File
@@ -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)
+25
View File
@@ -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)
@@ -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
@@ -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)
@@ -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"]
@@ -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
@@ -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)
@@ -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),
}
+35
View File
@@ -1,5 +1,6 @@
from __future__ import annotations from __future__ import annotations
from datetime import datetime
from typing import Any, Dict, List, Optional from typing import Any, Dict, List, Optional
from fastapi import APIRouter, Depends, HTTPException, Query from fastapi import APIRouter, Depends, HTTPException, Query
@@ -39,6 +40,11 @@ class S1BatchUnpackRequest(BaseModel):
scan_before_unpack: bool = True scan_before_unpack: bool = True
class SourceMaterializeRequest(BaseModel):
target_root: Optional[str] = None
overwrite: bool = False
@router.get("/inventory/status") @router.get("/inventory/status")
async def get_asset_inventory_status( async def get_asset_inventory_status(
current_user: AuthUserORM = Depends(_get_current_user), current_user: AuthUserORM = Depends(_get_current_user),
@@ -199,6 +205,35 @@ async def unpack_sentinel1_source_asset(
return {"message": "Sentinel-1 unpack task queued", "task_id": task_id, "job_id": job_id} return {"message": "Sentinel-1 unpack task queued", "task_id": task_id, "job_id": job_id}
@router.post("/sources/{asset_id}/materialize")
async def materialize_source_asset(
asset_id: int,
request: Optional[SourceMaterializeRequest] = None,
admin_user: AuthUserORM = Depends(_require_admin),
db: AsyncSession = Depends(get_db),
):
_ = admin_user
asset = await db.get(SourceProductAssetORM, asset_id)
if asset is None:
raise HTTPException(status_code=404, detail="Source product asset not found.")
request_data = request or SourceMaterializeRequest()
try:
result = asset_inventory_service.materialize_source_asset(
asset,
target_root=request_data.target_root,
overwrite=bool(request_data.overwrite),
)
except Exception as exc:
raise HTTPException(status_code=400, detail=str(exc)) from exc
metadata = dict(asset.metadata_json or {})
metadata["last_materialized_dir"] = result.get("safe_dir") or result.get("target_dir")
metadata["last_materialized_at"] = datetime.utcnow().isoformat()
metadata["last_materialized_status"] = result.get("status")
asset.metadata_json = metadata
await db.commit()
return {"message": "Source asset materialized", "result": result}
@router.post("/inventory/unpack-sentinel1", status_code=202) @router.post("/inventory/unpack-sentinel1", status_code=202)
async def run_sentinel1_unpack_batch( async def run_sentinel1_unpack_batch(
request: Optional[S1BatchUnpackRequest] = None, request: Optional[S1BatchUnpackRequest] = None,
+4
View File
@@ -24,6 +24,10 @@ class FloodPreprocessRequest(BaseModel):
class FloodWaterExtractionRequest(BaseModel): class FloodWaterExtractionRequest(BaseModel):
scene_id: Optional[int] = Field(default=None, description="SARSceneGeoORM primary key") scene_id: Optional[int] = Field(default=None, description="SARSceneGeoORM primary key")
input_path: Optional[str] = Field(default=None, description="Direct GeoTIFF/ENVI input path") 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): class FloodPairSearchRequest(BaseModel):
+511 -22
View File
@@ -5,8 +5,10 @@ import hashlib
import os import os
import re import re
import shutil import shutil
import tarfile
import zipfile import zipfile
from datetime import datetime, timedelta from datetime import datetime, timedelta
from pathlib import PurePosixPath
from typing import Any, Callable, Dict, Iterable, List, Optional, Sequence, Tuple from typing import Any, Callable, Dict, Iterable, List, Optional, Sequence, Tuple
from geoalchemy2.shape import from_shape from geoalchemy2.shape import from_shape
@@ -30,6 +32,7 @@ from ..models import (
from ..utils import ( from ..utils import (
find_xml_file, find_xml_file,
normalize_satellite_family, normalize_satellite_family,
parse_gf3_l2_dirname,
parse_lt1_radar_filename, parse_lt1_radar_filename,
parse_xml_metadata, parse_xml_metadata,
) )
@@ -68,6 +71,8 @@ _LT1_ORBIT_RE = re.compile(
r"^(?P<satellite>LT1[A-Z]?)_GpsData_GAS_C_(?P<date>\d{8})\.txt$", r"^(?P<satellite>LT1[A-Z]?)_GpsData_GAS_C_(?P<date>\d{8})\.txt$",
re.IGNORECASE, re.IGNORECASE,
) )
_LT1_ARCHIVE_EXTS = (".tar.gz", ".tgz", ".zip", ".tar")
_GF3_ARCHIVE_EXTS = (".tar.gz", ".tgz", ".zip", ".tar")
def _parse_bool(value: Any, default: bool = False) -> bool: def _parse_bool(value: Any, default: bool = False) -> bool:
@@ -190,13 +195,167 @@ def _asset_uid(prefix: str, path: str) -> str:
def _strip_known_suffix(name: str) -> str: def _strip_known_suffix(name: str) -> str:
lower = name.lower() lower = name.lower()
for suffix in (".tar.gz", ".tgz"):
if lower.endswith(suffix):
return name[: -len(suffix)]
if lower.endswith(".zip"): if lower.endswith(".zip"):
return name[:-4] return name[:-4]
if lower.endswith(".tar"):
return name[:-4]
if lower.endswith(".safe"): if lower.endswith(".safe"):
return name[:-5] return name[:-5]
return name return name
def _has_archive_suffix(name: str, suffixes: Sequence[str]) -> bool:
lower = str(name or "").lower()
return any(lower.endswith(suffix) for suffix in suffixes)
def _archive_member_base_name(member_name: str) -> str:
text = str(member_name or "").replace("\\", "/").strip("/")
return PurePosixPath(text).name
def _archive_member_scene_name(member_name: str, fallback: str) -> str:
parts = [part for part in str(member_name or "").replace("\\", "/").split("/") if part]
for part in parts:
stem = _strip_known_suffix(part)
if parse_lt1_radar_filename(stem) or parse_gf3_l2_dirname(stem):
return stem
return fallback
def _archive_read_first_matching(path: str, predicate: Callable[[str], bool]) -> Tuple[Optional[str], Optional[bytes], List[str]]:
members: List[str] = []
if zipfile.is_zipfile(path):
with zipfile.ZipFile(path) as archive:
for info in archive.infolist():
name = info.filename
if info.is_dir():
continue
members.append(name)
if predicate(name):
return name, archive.read(info), members
return None, None, members
if tarfile.is_tarfile(path):
with tarfile.open(path, "r:*") as archive:
for member in archive:
name = member.name
if not member.isfile():
continue
members.append(name)
if predicate(name):
source = archive.extractfile(member)
if source is None:
continue
with source:
return name, source.read(), members
return None, None, members
return None, None, members
def _archive_list_matching(path: str, predicate: Callable[[str], bool], *, limit: int = 20) -> List[str]:
matches: List[str] = []
def _visit(name: str) -> None:
if predicate(name):
matches.append(name)
if zipfile.is_zipfile(path):
with zipfile.ZipFile(path) as archive:
for info in archive.infolist():
if info.is_dir():
continue
_visit(info.filename)
if len(matches) >= limit:
break
return matches
if tarfile.is_tarfile(path):
with tarfile.open(path, "r:*") as archive:
for member in archive:
if not member.isfile():
continue
_visit(member.name)
if len(matches) >= limit:
break
return matches
def _safe_archive_member_name(member_name: str, archive_path: str) -> str:
name = str(member_name or "").replace("\\", "/").strip("/")
if not name or name.startswith("../") or "/../" in f"/{name}/":
raise ValueError(f"Unsafe archive member path in {archive_path}: {member_name}")
if os.path.isabs(name) or os.path.splitdrive(name)[0]:
raise ValueError(f"Unsafe archive member path in {archive_path}: {member_name}")
return name
def _extract_archive_to_dir(archive_path: str, target_dir: str, *, overwrite: bool = False) -> Dict[str, Any]:
archive = _normalize_path(archive_path)
target = _normalize_path(target_dir)
if not os.path.isfile(archive):
raise FileNotFoundError(f"Archive does not exist: {archive}")
if os.path.exists(target):
if not overwrite:
return {"status": "EXISTS", "archive_path": archive, "target_dir": target, "extracted": False, "member_count": None}
shutil.rmtree(target)
tmp_dir = target + ".materialize_tmp"
if os.path.exists(tmp_dir):
shutil.rmtree(tmp_dir)
os.makedirs(tmp_dir, exist_ok=True)
extracted = 0
try:
if zipfile.is_zipfile(archive):
with zipfile.ZipFile(archive) as zip_obj:
for info in zip_obj.infolist():
rel_name = _safe_archive_member_name(info.filename, archive)
destination = os.path.abspath(os.path.join(tmp_dir, rel_name))
if not destination.startswith(os.path.abspath(tmp_dir) + os.sep):
raise ValueError(f"Unsafe ZIP member path: {info.filename}")
if info.is_dir():
os.makedirs(destination, exist_ok=True)
continue
os.makedirs(os.path.dirname(destination), exist_ok=True)
with zip_obj.open(info, "r") as source, open(destination, "wb") as target_stream:
shutil.copyfileobj(source, target_stream, length=1024 * 1024)
extracted += 1
elif tarfile.is_tarfile(archive):
with tarfile.open(archive, "r:*") as tar_obj:
for member in tar_obj:
rel_name = _safe_archive_member_name(member.name, archive)
destination = os.path.abspath(os.path.join(tmp_dir, rel_name))
if not destination.startswith(os.path.abspath(tmp_dir) + os.sep):
raise ValueError(f"Unsafe TAR member path: {member.name}")
if member.isdir():
os.makedirs(destination, exist_ok=True)
continue
if not member.isfile():
continue
source = tar_obj.extractfile(member)
if source is None:
continue
os.makedirs(os.path.dirname(destination), exist_ok=True)
with source, open(destination, "wb") as target_stream:
shutil.copyfileobj(source, target_stream, length=1024 * 1024)
extracted += 1
else:
raise ValueError(f"Unsupported archive format: {archive}")
if extracted <= 0:
raise OSError(f"Archive extraction produced no files: {archive}")
os.makedirs(os.path.dirname(target), exist_ok=True)
os.replace(tmp_dir, target)
return {"status": "EXTRACTED", "archive_path": archive, "target_dir": target, "extracted": True, "member_count": extracted}
finally:
if os.path.exists(tmp_dir):
shutil.rmtree(tmp_dir, ignore_errors=True)
def _parse_datetime_token(value: Optional[str]) -> Optional[datetime]: def _parse_datetime_token(value: Optional[str]) -> Optional[datetime]:
text = str(value or "").strip() text = str(value or "").strip()
if not text: if not text:
@@ -216,6 +375,105 @@ def _parse_datetime_token(value: Optional[str]) -> Optional[datetime]:
return None return None
def _xml_text_by_local_names(root: etree._Element, names: Sequence[str]) -> Optional[str]:
name_set = {str(item).lower() for item in names}
for element in root.iter():
local_name = etree.QName(element).localname.lower()
if local_name not in name_set:
continue
text = str(element.text or "").strip()
if text:
return text
return None
def _xml_text_under_local_path(root: etree._Element, parent_name: str, child_name: str) -> Optional[str]:
parent_key = parent_name.lower()
child_key = child_name.lower()
for parent in root.iter():
if etree.QName(parent).localname.lower() != parent_key:
continue
for child in parent.iter():
if child is parent:
continue
if etree.QName(child).localname.lower() == child_key:
text = str(child.text or "").strip()
if text:
return text
return None
def _xml_float(value: Optional[str]) -> Optional[float]:
try:
if value is None or str(value).strip() == "":
return None
return float(str(value).strip())
except (TypeError, ValueError):
return None
def _parse_radar_xml_metadata_bytes(data: bytes) -> Tuple[Optional[List[Tuple[float, float]]], Dict[str, Any]]:
parser = _xml_parser()
root = etree.fromstring(data, parser=parser)
corners: List[Tuple[float, float]] = []
for element in root.iter():
if etree.QName(element).localname.lower() != "scenecornercoord":
continue
lon = _xml_float(_xml_text_under_local_path(element, "sceneCornerCoord", "lon") or _xml_text_by_local_names(element, ["lon"]))
lat = _xml_float(_xml_text_under_local_path(element, "sceneCornerCoord", "lat") or _xml_text_by_local_names(element, ["lat"]))
if lon is not None and lat is not None:
corners.append((lon, lat))
coverage_polygon: Optional[List[Tuple[float, float]]] = None
if len(corners) >= 4:
coverage_polygon = corners[:4]
if coverage_polygon[0] != coverage_polygon[-1]:
coverage_polygon.append(coverage_polygon[0])
start_time = (
_xml_text_under_local_path(root, "start", "timeUTC")
or _xml_text_by_local_names(root, ["startTime", "start_time", "beginPosition"])
)
stop_time = (
_xml_text_under_local_path(root, "stop", "timeUTC")
or _xml_text_by_local_names(root, ["stopTime", "stop_time", "endPosition"])
)
center_lon = _xml_float(_xml_text_under_local_path(root, "sceneCenterCoord", "lon"))
center_lat = _xml_float(_xml_text_under_local_path(root, "sceneCenterCoord", "lat"))
metadata = {
"orbit_direction": (_xml_text_by_local_names(root, ["pass", "orbitDirection"]) or "").upper() or None,
"imaging_mode": _xml_text_under_local_path(root, "acquisitionInfo", "imagingMode")
or _xml_text_under_local_path(root, "orderInfo", "imagingMode")
or _xml_text_by_local_names(root, ["imagingMode"]),
"polarization": _xml_text_under_local_path(root, "acquisitionInfo", "polarisationMode")
or _xml_text_under_local_path(root, "polarisationList", "polLayer")
or _xml_text_under_local_path(root, "polList", "polLayer")
or _xml_text_by_local_names(root, ["polarisationMode", "polarization", "polarisation", "polLayer"]),
"receiving_station": _xml_text_under_local_path(root, "generationInfo", "receivingStation")
or _xml_text_by_local_names(root, ["receivingStation"]),
"satellite_mode": _xml_text_by_local_names(root, ["satelliteMode"]),
"orbit_circle": _xml_text_under_local_path(root, "missionInfo", "absOrbit")
or _xml_text_by_local_names(root, ["absOrbit", "absoluteOrbit"]),
"relative_orbit": _xml_text_under_local_path(root, "missionInfo", "relOrbit")
or _xml_text_by_local_names(root, ["relOrbit", "relativeOrbit"]),
"scene_center_lon": center_lon,
"scene_center_lat": center_lat,
"acquisition_time_utc": start_time,
"acquisition_stop_time_utc": stop_time,
"product_type": _xml_text_under_local_path(root, "imageDataInfo", "imageDataType")
or _xml_text_under_local_path(root, "orderInfo", "productVariant")
or _xml_text_by_local_names(root, ["productType", "imageDataType", "productVariant"]),
"image_data_format": _xml_text_under_local_path(root, "imageDataInfo", "imageDataFormat")
or _xml_text_by_local_names(root, ["imageDataFormat"]),
"product_level": _xml_text_by_local_names(root, ["productLevel", "itemName"]),
"product_unique_id": _xml_text_by_local_names(root, ["logicalProductID", "sceneID", "productID"]),
"look_direction": (_xml_text_under_local_path(root, "acquisitionInfo", "lookDirection") or "").upper() or None,
"coverage_polygon": coverage_polygon,
}
return coverage_polygon, {key: value for key, value in metadata.items() if value not in (None, "", [])}
def _date_start_stop(date_yyyymmdd: str) -> Tuple[Optional[datetime], Optional[datetime]]: def _date_start_stop(date_yyyymmdd: str) -> Tuple[Optional[datetime], Optional[datetime]]:
try: try:
start = datetime.strptime(date_yyyymmdd, "%Y%m%d") start = datetime.strptime(date_yyyymmdd, "%Y%m%d")
@@ -416,6 +674,62 @@ def _parse_s1_safe_manifest(path: str) -> Dict[str, Any]:
} }
def _parse_lt1_archive_metadata(path: str) -> Dict[str, Any]:
archive_stem = _strip_known_suffix(os.path.basename(path))
xml_member, xml_data, members = _archive_read_first_matching(
path,
lambda name: _archive_member_base_name(name).lower().endswith(".meta.xml"),
)
tiff_members = _archive_list_matching(
path,
lambda name: _archive_member_base_name(name).lower().endswith((".tiff", ".tif")),
limit=8,
)
if not xml_member or not xml_data:
return {
"archive_parse_status": "MISSING_XML",
"archive_member_count_scanned": len(members),
"contained_tiff_members": tiff_members,
}
coverage_polygon, xml_meta = _parse_radar_xml_metadata_bytes(xml_data)
return {
"archive_parse_status": "OK",
"archive_xml_member": xml_member,
"archive_scene_name": _archive_member_scene_name(xml_member, archive_stem),
"contained_tiff_members": tiff_members,
"coverage_polygon": coverage_polygon,
**xml_meta,
}
def _parse_gf3_archive_metadata(path: str) -> Dict[str, Any]:
archive_stem = _strip_known_suffix(os.path.basename(path))
xml_member, xml_data, members = _archive_read_first_matching(
path,
lambda name: _archive_member_base_name(name).lower().endswith(".xml"),
)
quicklooks = _archive_list_matching(
path,
lambda name: _archive_member_base_name(name).lower().endswith((".jpg", ".jpeg", ".png", ".bmp", "_ql.tif", "_ql.tiff")),
limit=8,
)
if not xml_member or not xml_data:
return {
"archive_parse_status": "MISSING_XML",
"archive_member_count_scanned": len(members),
"quicklook_members": quicklooks,
}
coverage_polygon, xml_meta = _parse_radar_xml_metadata_bytes(xml_data)
return {
"archive_parse_status": "OK",
"archive_xml_member": xml_member,
"archive_scene_name": _archive_member_scene_name(xml_member, archive_stem),
"quicklook_members": quicklooks,
"coverage_polygon": coverage_polygon,
**xml_meta,
}
def _parse_s1_eof_header(path: str) -> Dict[str, Any]: def _parse_s1_eof_header(path: str) -> Dict[str, Any]:
try: try:
root = etree.parse(path, parser=_xml_parser()).getroot() root = etree.parse(path, parser=_xml_parser()).getroot()
@@ -443,6 +757,7 @@ def _parse_source_entry(path: str, root: ManagedRootORM) -> Optional[Dict[str, A
lower_name = name.lower() lower_name = name.lower()
stat = _stat_path(path) stat = _stat_path(path)
now = _utcnow() now = _utcnow()
name_stem = _strip_known_suffix(name)
if lower_name.endswith(".zip") and name.upper().startswith("S1"): if lower_name.endswith(".zip") and name.upper().startswith("S1"):
name_meta = _parse_s1_source_name(name) name_meta = _parse_s1_source_name(name)
@@ -459,6 +774,55 @@ def _parse_source_entry(path: str, root: ManagedRootORM) -> Optional[Dict[str, A
manifest_meta = {"manifest_parse_status": "FAILED", "manifest_parse_error": str(exc)} manifest_meta = {"manifest_parse_status": "FAILED", "manifest_parse_error": str(exc)}
return _build_s1_source_asset(path, root, name_meta, manifest_meta, stat, parse_status, parse_error, now) return _build_s1_source_asset(path, root, name_meta, manifest_meta, stat, parse_status, parse_error, now)
if _has_archive_suffix(name, _LT1_ARCHIVE_EXTS) and name_stem.upper().startswith("LT1"):
parsed = parse_lt1_radar_filename(name_stem)
if not parsed:
return None
parse_status = "OK"
parse_error = None
archive_meta: Dict[str, Any] = {}
try:
archive_meta = _parse_lt1_archive_metadata(path)
if archive_meta.get("archive_parse_status") != "OK":
parse_status = "PARTIAL"
parse_error = str(archive_meta.get("archive_parse_status") or "archive metadata incomplete")
except Exception as exc:
parse_status = "PARTIAL"
parse_error = str(exc)
archive_meta = {"archive_parse_status": "FAILED", "archive_parse_error": str(exc)}
return _build_lt1_source_asset(
path,
root,
parsed,
archive_meta,
archive_meta.get("coverage_polygon"),
stat,
now,
source_format="LT1_ARCHIVE",
archive_path=path,
parser_name="lt1_archive_metadata",
parse_status=parse_status,
parse_error=parse_error,
)
if _has_archive_suffix(name, _GF3_ARCHIVE_EXTS) and name_stem.upper().startswith("GF3"):
parsed = parse_gf3_l2_dirname(name_stem)
if not parsed:
return None
parse_status = "OK"
parse_error = None
archive_meta = {}
try:
archive_meta = _parse_gf3_archive_metadata(path)
if archive_meta.get("archive_parse_status") != "OK":
parse_status = "PARTIAL"
parse_error = str(archive_meta.get("archive_parse_status") or "archive metadata incomplete")
except Exception as exc:
parse_status = "PARTIAL"
parse_error = str(exc)
archive_meta = {"archive_parse_status": "FAILED", "archive_parse_error": str(exc)}
return _build_gf3_archive_asset(path, root, parsed, archive_meta, stat, parse_status, parse_error, now)
if lower_name.endswith(".safe") and os.path.isdir(path) and name.upper().startswith("S1"): if lower_name.endswith(".safe") and os.path.isdir(path) and name.upper().startswith("S1"):
name_meta = _parse_s1_source_name(name) name_meta = _parse_s1_source_name(name)
if not name_meta: if not name_meta:
@@ -569,6 +933,12 @@ def _build_lt1_source_asset(
coverage_polygon: Optional[List[Tuple[float, float]]], coverage_polygon: Optional[List[Tuple[float, float]]],
stat: Dict[str, Optional[float]], stat: Dict[str, Optional[float]],
now: datetime, now: datetime,
*,
source_format: str = "LT1_DIR",
archive_path: Optional[str] = None,
parser_name: str = "lt1_source_directory",
parse_status: str = "OK",
parse_error: Optional[str] = None,
) -> Dict[str, Any]: ) -> Dict[str, Any]:
metadata = dict(parsed) metadata = dict(parsed)
metadata.update({key: value for key, value in xml_meta.items() if value not in (None, "")}) metadata.update({key: value for key, value in xml_meta.items() if value not in (None, "")})
@@ -582,35 +952,99 @@ def _build_lt1_source_asset(
return { return {
"asset_uid": _asset_uid("source", path), "asset_uid": _asset_uid("source", path),
"logical_product_uid": os.path.basename(path), "logical_product_uid": _strip_known_suffix(os.path.basename(path)),
"satellite_family": normalize_satellite_family(satellite), "satellite_family": normalize_satellite_family(satellite),
"satellite": satellite, "satellite": satellite,
"source_format": "LT1_DIR", "source_format": source_format,
"product_type": parsed.get("product_type"), "product_type": xml_meta.get("product_type") or parsed.get("product_type"),
"product_level": parsed.get("product_level"), "product_level": xml_meta.get("product_level") or parsed.get("product_level"),
"imaging_mode": parsed.get("imaging_mode"), "imaging_mode": xml_meta.get("imaging_mode") or parsed.get("imaging_mode"),
"polarization": parsed.get("polarization"), "polarization": xml_meta.get("polarization") or parsed.get("polarization"),
"absolute_orbit": parsed.get("orbit_circle"), "absolute_orbit": xml_meta.get("orbit_circle") or parsed.get("orbit_circle"),
"relative_orbit": None, "relative_orbit": xml_meta.get("relative_orbit"),
"orbit_direction": xml_meta.get("orbit_direction") or parsed.get("orbit_direction"), "orbit_direction": xml_meta.get("orbit_direction") or parsed.get("orbit_direction"),
"acquisition_start_time_utc": None, "acquisition_start_time_utc": _parse_datetime_token(xml_meta.get("acquisition_time_utc")),
"acquisition_stop_time_utc": None, "acquisition_stop_time_utc": _parse_datetime_token(xml_meta.get("acquisition_stop_time_utc")),
"imaging_date": imaging_date, "imaging_date": imaging_date,
"root_ref_id": root.id, "root_ref_id": root.id,
"root_path": root.path, "root_path": root.path,
"file_path": path, "file_path": path,
"archive_path": None, "archive_path": archive_path,
"path_kind": _path_kind(path), "path_kind": _path_kind(path),
"file_name": os.path.basename(path), "file_name": os.path.basename(path),
"file_stem": os.path.basename(path), "file_stem": _strip_known_suffix(os.path.basename(path)),
"file_ext": "", "file_ext": os.path.splitext(path)[1].lower(),
"size_bytes": stat.get("size_bytes"), "size_bytes": stat.get("size_bytes"),
"mtime_epoch": stat.get("mtime_epoch"), "mtime_epoch": stat.get("mtime_epoch"),
"checksum_status": "NOT_COMPUTED", "checksum_status": "NOT_COMPUTED",
"parser_name": "lt1_source_directory", "parser_name": parser_name,
"parser_version": PARSER_VERSION, "parser_version": PARSER_VERSION,
"parse_status": "OK", "parse_status": parse_status,
"parse_error": None, "parse_error": parse_error,
"parsed_at": now,
"metadata_json": _json_safe(metadata),
"is_active": True,
"missing_since": None,
"updated_at": now,
}
def _build_gf3_archive_asset(
path: str,
root: ManagedRootORM,
parsed: Dict[str, Any],
archive_meta: Dict[str, Any],
stat: Dict[str, Optional[float]],
parse_status: str,
parse_error: Optional[str],
now: datetime,
) -> Dict[str, Any]:
coverage_polygon = archive_meta.get("coverage_polygon")
metadata = dict(parsed)
metadata.update({key: value for key, value in archive_meta.items() if value not in (None, "")})
metadata["coverage_polygon"] = coverage_polygon
metadata["coverage_bbox"] = _bbox_from_polygon(coverage_polygon)
centroid_lon, centroid_lat = _centroid_from_polygon(coverage_polygon)
metadata["scene_center_lon"] = parsed.get("scene_center_lon") if parsed.get("scene_center_lon") is not None else centroid_lon
metadata["scene_center_lat"] = parsed.get("scene_center_lat") if parsed.get("scene_center_lat") is not None else centroid_lat
start_time = _parse_datetime_token(archive_meta.get("acquisition_time_utc"))
stop_time = _parse_datetime_token(archive_meta.get("acquisition_stop_time_utc"))
imaging_date = parsed.get("imaging_date")
if not imaging_date and start_time:
imaging_date = start_time.strftime("%Y%m%d")
stem = _strip_known_suffix(os.path.basename(path))
return {
"asset_uid": _asset_uid("source", path),
"logical_product_uid": archive_meta.get("product_unique_id") or stem,
"satellite_family": "GF3",
"satellite": "GF3",
"source_format": "GF3_ARCHIVE",
"product_type": archive_meta.get("product_type") or parsed.get("product_type") or "L1A",
"product_level": archive_meta.get("product_level") or parsed.get("product_level") or "L1A",
"imaging_mode": archive_meta.get("imaging_mode") or parsed.get("imaging_mode"),
"polarization": archive_meta.get("polarization") or parsed.get("polarization"),
"absolute_orbit": archive_meta.get("orbit_circle") or parsed.get("orbit_circle"),
"relative_orbit": archive_meta.get("relative_orbit"),
"orbit_direction": archive_meta.get("orbit_direction") or parsed.get("orbit_direction"),
"acquisition_start_time_utc": start_time,
"acquisition_stop_time_utc": stop_time,
"imaging_date": imaging_date,
"root_ref_id": root.id,
"root_path": root.path,
"file_path": path,
"archive_path": path,
"path_kind": _path_kind(path),
"file_name": os.path.basename(path),
"file_stem": stem,
"file_ext": os.path.splitext(path)[1].lower(),
"size_bytes": stat.get("size_bytes"),
"mtime_epoch": stat.get("mtime_epoch"),
"checksum_status": "NOT_COMPUTED",
"parser_name": "gf3_archive_metadata",
"parser_version": PARSER_VERSION,
"parse_status": parse_status,
"parse_error": parse_error,
"parsed_at": now, "parsed_at": now,
"metadata_json": _json_safe(metadata), "metadata_json": _json_safe(metadata),
"is_active": True, "is_active": True,
@@ -732,6 +1166,9 @@ def _iter_source_candidates(root_path: str) -> Iterable[str]:
yield _normalize_path(entry.path) yield _normalize_path(entry.path)
continue continue
stack.append(entry.path) stack.append(entry.path)
elif entry.is_file(follow_symlinks=False):
if entry.name.upper().startswith("S1") and entry.name.lower().endswith(".zip"):
yield _normalize_path(entry.path)
except OSError: except OSError:
continue continue
except OSError: except OSError:
@@ -749,8 +1186,17 @@ def _iter_s1_zip_candidates(root_path: str) -> Iterable[str]:
if entry.is_dir(follow_symlinks=False): if entry.is_dir(follow_symlinks=False):
stack.append(entry.path) stack.append(entry.path)
elif entry.is_file(follow_symlinks=False): elif entry.is_file(follow_symlinks=False):
if entry.name.upper().startswith("S1") and entry.name.lower().endswith(".zip"): name_upper = entry.name.upper()
stem_upper = _strip_known_suffix(entry.name).upper()
if name_upper.startswith("S1") and entry.name.lower().endswith(".zip"):
yield _normalize_path(entry.path) yield _normalize_path(entry.path)
continue
if stem_upper.startswith("LT1") and _has_archive_suffix(entry.name, _LT1_ARCHIVE_EXTS):
yield _normalize_path(entry.path)
continue
if stem_upper.startswith("GF3") and _has_archive_suffix(entry.name, _GF3_ARCHIVE_EXTS):
yield _normalize_path(entry.path)
continue
except OSError: except OSError:
continue continue
except OSError: except OSError:
@@ -860,6 +1306,13 @@ def _insar_source_ready(row: Dict[str, Any], coverage_polygon: Optional[List[Tup
return True, None return True, None
def _image_data_format_for_source(row: Dict[str, Any]) -> str:
source_format = str(row.get("source_format") or "").upper()
if source_format in {"S1_ZIP", "LT1_ARCHIVE", "GF3_ARCHIVE"}:
return "ARCHIVE"
return "DIRECTORY"
class AssetInventoryService: class AssetInventoryService:
async def _progress(self, task_id: Optional[str], message: str, progress: int) -> None: async def _progress(self, task_id: Optional[str], message: str, progress: int) -> None:
if not task_id: if not task_id:
@@ -876,7 +1329,12 @@ class AssetInventoryService:
type_set = {str(item or "").strip().lower() for item in (inventory_types or []) if str(item or "").strip()} type_set = {str(item or "").strip().lower() for item in (inventory_types or []) if str(item or "").strip()}
roles: List[str] = [] roles: List[str] = []
if not type_set or "source_product" in type_set or "source" in type_set: if not type_set or "source_product" in type_set or "source" in type_set:
roles.append("source_product_pool") roles.extend(
[
"source_product_pool",
"source_pool_gf3_archive",
]
)
if not type_set or "orbit_asset" in type_set or "orbit" in type_set: if not type_set or "orbit_asset" in type_set or "orbit" in type_set:
roles.append("orbit_asset_pool") roles.append("orbit_asset_pool")
stmt = ( stmt = (
@@ -921,7 +1379,7 @@ class AssetInventoryService:
for index, root in enumerate(roots, start=1): for index, root in enumerate(roots, start=1):
progress = 5 + int((index - 1) / max(1, total_roots) * 75) progress = 5 + int((index - 1) / max(1, total_roots) * 75)
await self._progress(task_id, f"Scanning {root.display_name}: {root.path}", progress) await self._progress(task_id, f"Scanning {root.display_name}: {root.path}", progress)
if root.root_role == "source_product_pool": if root.root_role in {"source_product_pool", "source_pool_gf3_archive"}:
result = await self.scan_source_root(db, root) result = await self.scan_source_root(db, root)
totals["source_roots"] += 1 totals["source_roots"] += 1
totals["source_assets"] += int(result.get("asset_count") or 0) totals["source_assets"] += int(result.get("asset_count") or 0)
@@ -1453,6 +1911,39 @@ class AssetInventoryService:
"member_count": len(names), "member_count": len(names),
} }
def materialize_source_asset(
self,
asset: SourceProductAssetORM,
*,
target_root: Optional[str] = None,
overwrite: bool = False,
) -> Dict[str, Any]:
source_format = str(asset.source_format or "").upper()
source_path = _normalize_path(str(asset.archive_path or asset.file_path or ""))
if not source_path:
raise ValueError("Source asset path is empty.")
if source_format == "S1_ZIP":
return self.unpack_sentinel1_archive(source_path, target_root=target_root, overwrite=overwrite)
if source_format not in {"LT1_ARCHIVE", "GF3_ARCHIVE"}:
if os.path.isdir(source_path):
return {
"status": "DIRECTORY_READY",
"source_path": source_path,
"target_dir": source_path,
"extracted": False,
"source_format": source_format,
}
raise ValueError(f"Source format is not materializable from archive: {source_format}")
requested_root = _normalize_path(target_root or "")
if not requested_root:
requested_root = _normalize_path(os.path.join(settings.PYINT_WORK_ROOT, "source_materialized", source_format.lower()))
scene_name = _strip_known_suffix(os.path.basename(source_path))
target_dir = os.path.join(requested_root, scene_name)
result = _extract_archive_to_dir(source_path, target_dir, overwrite=overwrite)
result["source_format"] = source_format
return result
async def run_sentinel1_unpack_task(self, task_id: str, payload: Optional[Dict[str, Any]] = None) -> None: async def run_sentinel1_unpack_task(self, task_id: str, payload: Optional[Dict[str, Any]] = None) -> None:
payload = payload if isinstance(payload, dict) else {} payload = payload if isinstance(payload, dict) else {}
asset_id = payload.get("asset_id") asset_id = payload.get("asset_id")
@@ -1831,7 +2322,7 @@ class AssetInventoryService:
"product_type": row.get("product_type"), "product_type": row.get("product_type"),
"source_product_token": metadata.get("filename_class_token") or metadata.get("source_product_token"), "source_product_token": metadata.get("filename_class_token") or metadata.get("source_product_token"),
"image_data_type": "COMPLEX", "image_data_type": "COMPLEX",
"image_data_format": "ZIP" if row.get("source_format") == "S1_ZIP" else "DIRECTORY", "image_data_format": _image_data_format_for_source(row),
"product_variant": metadata.get("product_variant"), "product_variant": metadata.get("product_variant"),
"product_level": row.get("product_level"), "product_level": row.get("product_level"),
"product_unique_id": row.get("logical_product_uid"), "product_unique_id": row.get("logical_product_uid"),
@@ -2149,7 +2640,6 @@ class AssetInventoryService:
await db.execute( await db.execute(
select(func.count(SourceProductAssetORM.id)).where( select(func.count(SourceProductAssetORM.id)).where(
SourceProductAssetORM.is_active == True, # noqa: E712 SourceProductAssetORM.is_active == True, # noqa: E712
SourceProductAssetORM.source_format != "S1_ZIP",
) )
) )
).scalar_one() ).scalar_one()
@@ -2202,7 +2692,6 @@ class AssetInventoryService:
filters = [] filters = []
if not include_inactive: if not include_inactive:
filters.append(SourceProductAssetORM.is_active == True) # noqa: E712 filters.append(SourceProductAssetORM.is_active == True) # noqa: E712
filters.append(SourceProductAssetORM.source_format != "S1_ZIP")
if satellite_family: if satellite_family:
filters.append(SourceProductAssetORM.satellite_family == satellite_family.upper()) filters.append(SourceProductAssetORM.satellite_family == satellite_family.upper())
if satellite: if satellite:
+98 -4
View File
@@ -134,6 +134,72 @@ def _scene_analysis_path(scene: SARSceneGeoORM | None) -> str | None:
return scene.analysis_tif_path 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]]: 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.""" """Resolve region/GeoJSON AOI using the same parser as the management search page."""
aoi_geojson = getattr(req, "aoi_geojson", None) 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]: async def submit_water_extraction(req: Any, db: AsyncSession) -> dict[str, Any]:
input_path = req.input_path input_path = req.input_path
scene_id = req.scene_id 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: if scene_id:
scene = await db.get(SARSceneGeoORM, scene_id) scene = await db.get(SARSceneGeoORM, scene_id)
if not scene: if not scene:
raise HTTPException(status_code=404, detail=f"SARSceneGeoORM id={scene_id} not found") raise HTTPException(status_code=404, detail=f"SARSceneGeoORM id={scene_id} not found")
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) input_path = _scene_analysis_path(scene)
if not input_path: if processor != "gf3_hh_hv" and not input_path:
raise HTTPException(status_code=400, detail="Scene has no analysis-ready GeoTIFF") 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") raise HTTPException(status_code=400, detail="scene_id or input_path is required")
extraction = WaterExtractionORM( extraction = WaterExtractionORM(
scene_id=scene_id, scene_id=scene_id,
processor=getattr(req, "processor", None) or "otsu", processor=processor,
input_path=input_path, 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", status="PENDING",
) )
db.add(extraction) 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, job_type=JOB_TYPE_WATER_DETECT,
task_type=f"FLOOD_WATER_EXTRACTION_{extraction_id}", task_type=f"FLOOD_WATER_EXTRACTION_{extraction_id}",
task_name=f"Flood water extraction id={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(): async with db.begin():
queued_extraction = await db.get(WaterExtractionORM, extraction_id) queued_extraction = await db.get(WaterExtractionORM, extraction_id)
@@ -0,0 +1,246 @@
"""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")
use_dltb = False
dltb_cache_dir = None
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="off",
water_vector=[],
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,
"water_vector": [],
"paddy_vector": [],
},
"runtime": {
"prior_inputs_enabled": False,
"dltb_enabled": False,
"deep_learning_enabled": False,
"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,
}
+33 -2
View File
@@ -3490,6 +3490,7 @@ async def _handle_flood_detection(job: SystemJobORM) -> None:
async def _handle_water_detect(job: SystemJobORM) -> None: async def _handle_water_detect(job: SystemJobORM) -> None:
"""水体检测 job handlerOtsu + DEM + 形态学 + 连通分量)。""" """水体检测 job handlerOtsu + 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 from .water_extraction_service import run_otsu_water_extraction
payload = job.payload or {} 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" model_name = "WaterExtractionORM" if use_extraction_table else "WaterDetectionORM"
raise ValueError(f"{model_name} id={record_id} 不存在") raise ValueError(f"{model_name} id={record_id} 不存在")
input_path = det.input_path 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" det.status = "RUNNING"
if use_extraction_table and hasattr(det, "task_id"): if use_extraction_table and hasattr(det, "task_id"):
det.task_id = job.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 mirror.task_id = job.task_id
await db.commit() 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") raise ValueError("水体检测缺少输入路径 input_path")
output_name = f"water_extraction_{record_id}" if use_extraction_table else f"water_detect_{record_id}" 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_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) os.makedirs(output_dir, exist_ok=True)
await task_service.update_task(job.task_id, progress=10, message="启动水体检测算法...") await task_service.update_task(job.task_id, progress=10, message="启动水体检测算法...")
def _run() -> Dict[str, Any]: 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( return run_otsu_water_extraction(
input_path=input_path, input_path=input_path,
output_dir=output_dir, output_dir=output_dir,
@@ -3557,11 +3578,19 @@ async def _handle_water_detect(job: SystemJobORM) -> None:
if det: if det:
if result.get("ok"): if result.get("ok"):
det.output_path = result.get("output_path") 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_area_km2 = result.get("water_area_km2")
det.water_pixel_count = result.get("water_pixel_count") det.water_pixel_count = result.get("water_pixel_count")
if use_extraction_table: if use_extraction_table:
det.processor = result.get("processor") or det.processor or "otsu" det.processor = result.get("processor") or det.processor or "otsu"
det.threshold_value = result.get("threshold_value") det.threshold_value = result.get("threshold_value")
result_metadata = result.get("metadata_json")
if isinstance(result_metadata, dict):
det.metadata_json = result_metadata
else:
det.metadata_json = { det.metadata_json = {
"legacy_otsu_threshold_db": result.get("otsu_threshold_db"), "legacy_otsu_threshold_db": result.get("otsu_threshold_db"),
"value_transform": result.get("value_transform"), "value_transform": result.get("value_transform"),
@@ -3578,6 +3607,8 @@ async def _handle_water_detect(job: SystemJobORM) -> None:
mirror = await db.get(WaterExtractionORM, int(record_id)) mirror = await db.get(WaterExtractionORM, int(record_id))
if mirror: if mirror:
mirror.output_path = det.output_path 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_area_km2 = det.water_area_km2
mirror.water_pixel_count = det.water_pixel_count mirror.water_pixel_count = det.water_pixel_count
mirror.threshold_value = result.get("threshold_value") or result.get("otsu_threshold_db") mirror.threshold_value = result.get("threshold_value") or result.get("otsu_threshold_db")
@@ -458,6 +458,35 @@ def _build_root_specs_from_settings() -> List[RootSpec]:
owner_engine="pyint", owner_engine="pyint",
) )
) )
specs.extend(
_iter_single_root_specs(
env_var="TASK_POOL_ROOT",
path=settings.TASK_POOL_ROOT,
root_role="task_pool_root",
display_name="Task Pool Root",
scan_mode="workspace",
)
)
specs.extend(
_iter_single_root_specs(
env_var="DINSAR_TASK_POOL_ROOT",
path=settings.DINSAR_TASK_POOL_ROOT,
root_role="task_pool_dinsar",
display_name="D-InSAR Task Pool",
scan_mode="workspace",
owner_engine="dinsar",
)
)
specs.extend(
_iter_single_root_specs(
env_var="SBAS_TASK_POOL_ROOT",
path=settings.SBAS_TASK_POOL_ROOT,
root_role="task_pool_sbas",
display_name="SBAS Task Pool",
scan_mode="workspace",
owner_engine="sbas",
)
)
specs.extend( specs.extend(
_iter_single_root_specs( _iter_single_root_specs(
env_var="GAMMA_SBAS_WORK_ROOT", env_var="GAMMA_SBAS_WORK_ROOT",
@@ -210,3 +210,22 @@ Task_20250101_20250201
- Task_Pool 根路径需要明确配置项,建议新增 `DINSAR_TASK_POOL_ROOT`,默认 `D:\Task_Pool\DInSAR` - Task_Pool 根路径需要明确配置项,建议新增 `DINSAR_TASK_POOL_ROOT`,默认 `D:\Task_Pool\DInSAR`
- 中间文件清理需要先确认每个引擎的“可删目录”和“必须保留资产”清单,不能用一套规则覆盖全部。 - 中间文件清理需要先确认每个引擎的“可删目录”和“必须保留资产”清单,不能用一套规则覆盖全部。
- 前端 grouped 结果视图需要兼容旧 flat API 一段时间,避免已有页面一次性断裂。 - 前端 grouped 结果视图需要兼容旧 flat API 一段时间,避免已有页面一次性断裂。
## 2026-06-15 Task_Pool Materialize Update
Current direction:
- `TASK_POOL_ROOT` defaults to `D:\Task_Pool`.
- `DINSAR_TASK_POOL_ROOT` defaults to `D:\Task_Pool\DInSAR`.
- `SBAS_TASK_POOL_ROOT` defaults to `D:\Task_Pool\SBAS`.
- D-InSAR distribution materializes source inputs inside the task folder:
- directory sources are copied into `master/` and `slave/`;
- `S1_ZIP`, `LT1_ARCHIVE`, and other supported archives are extracted into `master/` and `slave/`;
- staged orbit files go into `orbit/`.
- Engines must consume local Task_Pool paths, not UNC source archive paths.
- `.dinsar_pair.json` records `source_materialization` so cleanup can distinguish copied directories, extracted archives, and staged files.
Cleanup implication:
- `master/`, `slave/`, `orbit/`, and engine `work/` folders are local materialized inputs/workspace and may be cleaned after all required results are registered.
- `publish/`, manifests, result assets, previews, and catalog metadata are preserved.
@@ -0,0 +1,163 @@
# 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 from the experiment workspace.
- `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.
Production policy:
- Do not use DLTB, hydro, water-vector or paddy-vector priors.
- Do not use deep-learning checkpoints or U-Net inference.
- Use the current HH/HV machine-learning-style threshold, candidate and morphology workflow.
## 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.
Current integration does not use DLTB priors. The GF-3 HH/HV processor runs from SAR backscatter, morphology and optional vector/DEM inputs only.
Not transferred:
- `D:\Code\Water\data\priors\dltb_cache\heilongjiang`
- Current size observed: 25 files, about 8.5 GB.
- Hydro prior vectors under `data/priors/hydro`.
- Water/paddy vector priors.
Optional runtime asset:
- 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_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.
- DLTB, hydro, water-vector and paddy-vector priors are not part of the current workflow.
- Legacy AI4G U-Net/deep-learning checkpoints are not part of the current workflow.
+2
View File
@@ -30,6 +30,8 @@
- [DINSAR_TASK_POOL_THREE_ENGINE_REFACTOR_20260614.md](DINSAR_TASK_POOL_THREE_ENGINE_REFACTOR_20260614.md) - [DINSAR_TASK_POOL_THREE_ENGINE_REFACTOR_20260614.md](DINSAR_TASK_POOL_THREE_ENGINE_REFACTOR_20260614.md)
D-InSAR 保留 ENVI/SARscape、LandSAR、Gamma/PyINT 三引擎,退出 ISCE2,统一 Task_Pool、结果聚合和中间文件清理的当前设计。 D-InSAR 保留 ENVI/SARscape、LandSAR、Gamma/PyINT 三引擎,退出 ISCE2,统一 Task_Pool、结果聚合和中间文件清理的当前设计。
- [UNC_SOURCE_ARCHIVE_AND_MATERIALIZE_DESIGN_20260615.md](UNC_SOURCE_ARCHIVE_AND_MATERIALIZE_DESIGN_20260615.md)
UNC/SMB 源压缩包管理、包内 XML/manifest 资产化、本地 materialize 和 D-InSAR/SBAS 生产边界。
- [DINSAR_PRODUCTION_CORES_OVERVIEW.md](DINSAR_PRODUCTION_CORES_OVERVIEW.md) - [DINSAR_PRODUCTION_CORES_OVERVIEW.md](DINSAR_PRODUCTION_CORES_OVERVIEW.md)
旧版 ENVI/SARscape、ISCE2、Gamma/PyINT D-InSAR 生产核心说明。ISCE2 相关内容仅作历史背景。 旧版 ENVI/SARscape、ISCE2、Gamma/PyINT D-InSAR 生产核心说明。ISCE2 相关内容仅作历史背景。
+25
View File
@@ -252,3 +252,28 @@ GET /api/sbas-insar-products/{product_id}/assets/{asset_id}
- 选中候选可生成 Stack Manifest - 选中候选可生成 Stack Manifest
- 结果发布支持 GeoTIFF、预览图、监测点曲线和点矢量下载; - 结果发布支持 GeoTIFF、预览图、监测点曲线和点矢量下载;
- 前端生产页和结果页构建通过。 - 前端生产页和结果页构建通过。
## 2026-06-15 SBAS Task_Pool Update
Recommended SBAS work root:
```text
SBAS_TASK_POOL_ROOT=D:\Task_Pool\SBAS
GAMMA_SBAS_WORK_ROOT=D:\Task_Pool\SBAS
```
SBAS has its own Task_Pool sub-root and must not reuse `D:\Task_Pool\DInSAR`.
Recommended run layout:
```text
D:\Task_Pool\SBAS\<stack_task>
├─ task_manifest.json
├─ sbas_stack_manifest.json
├─ sources
├─ orbits
├─ work
└─ publish
```
Source archives remain on UNC. Selected scenes and orbit files are materialized under the SBAS task directory before Gamma/LandSAR execution. Cleanup may remove `sources`, `orbits`, and `work` after result registration, but must preserve manifests, `publish`, previews, and catalog assets.
@@ -0,0 +1,307 @@
# UNC Source Archive and Local Materialize Design
## Decision
UNC/SMB storage is treated as the source archive pool. Production engines should not use UNC paths as their working input. D-InSAR, SBAS, Gamma/PyINT, LandSAR, and SARscape should consume local materialized task inputs.
This keeps the 20 TB storage useful for long-term source management while protecting production from SMB disconnects, credential scope, WSL path conversion, and external engine UNC compatibility.
## Current Implementation
- Source asset inventory now recognizes archive assets:
- `S1_ZIP`
- `LT1_ARCHIVE`
- `GF3_ARCHIVE`
- Sentinel-1 ZIP manifest parsing already reads `manifest.safe` directly from the ZIP.
- LT-1 archive parsing reads `*.meta.xml` directly from `.zip`, `.tar`, `.tar.gz`, or `.tgz` and records contained TIFF members.
- GF3 archive parsing reads the first XML member directly from `.zip`, `.tar`, `.tar.gz`, or `.tgz` and records quicklook-like members when present.
- `GF3_ARCHIVE_SOURCE_DIRS` roots are included in asset inventory scans as source pools.
- Source asset listing and inventory counts now include archive assets instead of hiding `S1_ZIP`.
- A generic source materialize endpoint exists:
- `POST /api/assets/sources/{asset_id}/materialize`
- `S1_ZIP` uses the existing Sentinel-1 SAFE unpacker.
- `LT1_ARCHIVE` and `GF3_ARCHIVE` extract to a local materialized directory.
- Directory assets return `DIRECTORY_READY`.
Default source materialization is task-scoped. D-InSAR and SBAS callers should pass a Task_Pool target directory:
```text
D:\Task_Pool\DInSAR\<task>\master
D:\Task_Pool\DInSAR\<task>\slave
D:\Task_Pool\SBAS\<stack>\sources\<YYYYMMDD>
```
The generic materialize endpoint still accepts `target_root` for ad hoc checks. Production callers must provide a Task_Pool destination.
## Production Boundary
D-InSAR and SBAS should store source asset references in task/run manifests, then materialize selected inputs into the run directory before engine execution.
Required next integration points:
- D-InSAR Task_Pool publishing:
- store `source_product_asset_id`, `archive_path`, `source_format`;
- materialize master/slave archive assets into the task directory before engine dispatch.
- Gamma/PyINT:
- always consume local materialized paths because WSL conversion rejects or cannot reliably map UNC paths.
- LandSAR and ENVI/SARscape:
- prefer local materialized paths even when Windows can see UNC, to avoid external engine path and credential issues.
- SBAS:
- stack discovery can use archive metadata;
- selected scenes must be materialized into the SBAS `RAW`/input structure before Gamma commands such as `par_LT1_SLC`.
## GF3 Management
GF3 has two asset layers:
- `GF3_ARCHIVE`: original source archive, suitable for UNC source management and migration tracking.
- GF3 SARscape standardized L2: production result/analysis-ready layer, used for map footprint, preview, radar data management, and water extraction.
Do not replace standardized L2 management with raw archive management. Archive assets should link migration and production status; previews and water extraction should continue to consume standardized L2/analysis-ready products.
## Migration Guidance
1. Register UNC roots first and scan inventory.
2. Verify archive asset counts and parse status.
3. Keep existing local standardized results and D-InSAR/SBAS products in place.
4. Move source archives to UNC and update root configuration.
5. Only after inventory and materialize tests pass, switch D-InSAR/SBAS publishing to archive asset references.
Production safety rule: if a run cannot materialize every selected source asset locally, the run must fail before invoking the engine.
## Recommended UNC Layout
The current deployment uses two SMB shares:
```text
\\DESKTOP-N16HJ84\InSAR_Storage_1
\\DESKTOP-N16HJ84\InSAR_Storage_2
```
Recommended source archive layout:
```text
\\DESKTOP-N16HJ84\InSAR_Storage_1
└─ GaoFen-3
├─ 20260513
│ └─ GF3_*.tar.gz
└─ 20260514
\\DESKTOP-N16HJ84\InSAR_Storage_2
├─ LuTan-1
│ └─ Archive
│ ├─ 20260513
│ │ └─ LT1*.tar.gz / LT1*.tgz / LT1*.zip / LT1*.tar
│ └─ 20260514
├─ Sentinel-1
│ └─ Archive
│ ├─ 20260513
│ │ └─ S1*.zip
│ └─ 20260514
└─ Orbit
├─ LuTan-1
│ ├─ LT1A_GpsData_GAS_C_YYYYMMDD.txt
│ └─ LT1B_GpsData_GAS_C_YYYYMMDD.txt
└─ Sentinel-1
└─ S1*.EOF
```
Recommended local Task_Pool layout:
```text
D:\Task_Pool
├─ DInSAR
│ └─ <pair_task>
│ ├─ task_manifest.json
│ ├─ .dinsar_pair.json
│ ├─ master
│ ├─ slave
│ ├─ orbit
│ ├─ work
│ └─ publish
└─ SBAS
└─ <stack_task>
├─ task_manifest.json
├─ sbas_stack_manifest.json
├─ sources
├─ orbits
├─ work
└─ publish
```
Date folders are optional for the scanner because source and orbit inventory recurse through configured roots. They are recommended for operator readability and migration checks.
## Current Local Configuration Example
The local `.env` should keep legacy local roots and UNC roots side by side during migration:
```text
SOURCE_PRODUCT_DIRS=D:\LuTan1_Image_Pool;D:\Sentinel1_Image_Pool_ZIP;\\DESKTOP-N16HJ84\InSAR_Storage_2\LuTan-1\Archive;\\DESKTOP-N16HJ84\InSAR_Storage_2\Sentinel-1\Archive
ORBIT_SOURCE_DIRS=D:\LT1_data_lsarorbit;D:\Sentinel1_EOF_Pool;\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\LuTan-1;\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\Sentinel-1
GF3_ARCHIVE_SOURCE_DIRS=\\DESKTOP-N16HJ84\InSAR_Storage_1\GaoFen-3
TASK_POOL_ROOT=D:\Task_Pool
DINSAR_TASK_POOL_ROOT=D:\Task_Pool\DInSAR
SBAS_TASK_POOL_ROOT=D:\Task_Pool\SBAS
GAMMA_SBAS_WORK_ROOT=D:\Task_Pool\SBAS
```
Do not store SMB credentials in `.env`. Credentials should be stored in Windows Credential Manager for the account that runs the backend/worker service.
## Orbit Pool Contract
There are two different orbit concepts:
- `ORBIT_SOURCE_DIRS`: source inventory roots. These can be UNC and may be date-organized or flat.
- `ORBIT_POOL_ENVI` / `PYINT_ORBIT_POOL_TXT`: local production orbit pools. These should remain local disk paths.
LT-1 local production orbit pool should support both flat and satellite-split layouts:
```text
D:\orbit_pools\envi
├─ LT1A
│ └─ LT1A_GpsData_GAS_C_YYYYMMDD.txt
├─ LT1B
│ └─ LT1B_GpsData_GAS_C_YYYYMMDD.txt
└─ converted
└─ envi
```
The `LT1A` and `LT1B` names are satellite names, not product levels. ENVI/Gamma/PyINT/SBAS should use local orbit files copied or synchronized from `ORBIT_SOURCE_DIRS`; they should not be required to read UNC directly.
Sentinel-1 EOF files can be indexed from UNC. Gamma/PyINT/SBAS execution should stage required EOF files locally with the selected scenes.
## Migration Phases
### Phase 1: Source archive migration
Move or copy source archives only:
- LT-1 compressed scenes to `\\DESKTOP-N16HJ84\InSAR_Storage_2\LuTan-1\Archive\<YYYYMMDD>\`.
- Sentinel-1 ZIP scenes to `\\DESKTOP-N16HJ84\InSAR_Storage_2\Sentinel-1\Archive\<YYYYMMDD>\`.
- GF3 raw archives to `\\DESKTOP-N16HJ84\InSAR_Storage_1\GaoFen-3\<YYYYMMDD>\`.
Keep current local unpacked scene directories in place until D-InSAR and SBAS archive materialization have been tested.
### Phase 2: Orbit source migration
Copy orbit source files to UNC:
- LT-1 TXT files to `\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\LuTan-1\`.
- Sentinel-1 EOF files to `\\DESKTOP-N16HJ84\InSAR_Storage_2\Orbit\Sentinel-1\`.
Keep `ORBIT_POOL_ENVI` and `PYINT_ORBIT_POOL_TXT` local. Add a later sync/materialize step to populate local orbit pools from the indexed UNC source assets.
### Phase 3: Production cutover
After inventory scan verifies UNC assets:
1. D-InSAR Task_Pool stores source asset IDs and archive paths.
2. Task preparation materializes master/slave scenes and orbit files under `D:\Task_Pool\DInSAR\<task>`.
3. Engines run only against local Task_Pool paths.
4. Results register normally.
5. Local materialized inputs and intermediate products are eligible for cleanup after result registration.
### Phase 4: Retire old local source pools
Only after repeated D-InSAR/SBAS runs succeed from archive materialization:
- remove old local source roots from `SOURCE_PRODUCT_DIRS`;
- keep local work/result roots;
- keep standardized GF3 L2 products unless explicitly migrated and revalidated.
## Local Cleanup Design
After source archives are managed on UNC and production results are registered as assets, local disk can be treated as a cache/work area. Cleanup should be explicit and asset-aware.
### Keep Classes
Cleanup must never delete:
- configured UNC source archive roots;
- local or UNC orbit source roots;
- registered D-InSAR result assets;
- registered SBAS result assets;
- registered GF3 standardized L2 assets;
- `SAR_ANALYSIS_READY_ROOT` products and water extraction result assets;
- current pointers, manifests, previews, and catalog metadata needed to open results.
### Cleanup Classes
Cleanup may delete only these local classes after verification:
- materialized source inputs under `source_materialized`;
- D-InSAR Task_Pool copied inputs after every required engine run is registered;
- D-InSAR engine intermediate folders not listed in the result manifest;
- Gamma/PyINT temporary project work directories after result registration;
- SBAS `RAW`, `SLC`, `RSLC`, `MLI`, `DIFF`, `DIFF1`, script logs, and temporary staging after SBAS product registration;
- GF3 SARscape native intermediates only after standardized L2 registration and optional native-retention policy allows cleanup.
### Safety Contract
Every cleanup operation should run in two phases:
1. `preview`: enumerate candidate paths, classify each path, show size, last modified time, owning task/run/product, and keep/delete reason.
2. `execute`: delete only candidates from a persisted preview token or exact candidate list.
Deletion must require:
- path is inside an approved local work root;
- path is not inside any configured source archive root;
- path is not inside a result publish root unless the exact file is classified as intermediate;
- associated result or standardized asset is registered;
- candidate is older than a configurable minimum age;
- no active task references the path.
### Proposed API
```text
POST /api/maintenance/cleanup/preview
POST /api/maintenance/cleanup/execute
```
Preview request fields:
```json
{
"scope": "dinsar|sbas|gf3|materialized|all",
"root_ids": [],
"older_than_hours": 24,
"require_registered_result": true,
"include_task_pool_inputs": false
}
```
Preview response should include:
```json
{
"preview_id": "...",
"total_bytes": 0,
"candidates": [
{
"path": "D:\\production_runtime\\...",
"class": "materialized_source",
"owner": "task/run/product id",
"size_bytes": 0,
"eligible": true,
"reason": "registered_result_exists"
}
],
"blocked": []
}
```
### Recommended Defaults
- `materialized`: delete after 24 hours if no active task references it.
- `dinsar`: delete engine intermediates after result registration; keep Task_Pool inputs until all selected engines are complete or user opts in.
- `sbas`: delete heavy Gamma working directories after SBAS catalog registration and product assets exist.
- `gf3`: keep standardized L2; clean SARscape native only when `GF3_SARSCAPE_CLEAN_AFTER_SUCCESS=true` and standardized registration is confirmed.
### Implementation Order
1. Add read-only cleanup preview service.
2. Add path classification and approved-root checks.
3. Add execute endpoint with preview token.
4. Add frontend maintenance panel.
5. Wire D-InSAR/SBAS/GF3 run pages to show cleanup eligibility after successful registration.
+17 -2
View File
@@ -316,6 +316,7 @@ function SceneRow({ scene, readOnly, onShowMap, onExtractWater, onReset }) {
} }
function WaterResultRow({ item, onShowMap }) { function WaterResultRow({ item, onShowMap }) {
const thresholdValue = item.threshold_value ?? item.otsu_threshold_db;
return ( return (
<div style={rowStyle}> <div style={rowStyle}>
<div style={{ display: 'grid', gridTemplateColumns: '1fr auto', gap: 8, alignItems: 'center' }}> <div style={{ display: 'grid', gridTemplateColumns: '1fr auto', gap: 8, alignItems: 'center' }}>
@@ -324,6 +325,7 @@ function WaterResultRow({ item, onShowMap }) {
<strong>水体提取 #{item.id}</strong> <strong>水体提取 #{item.id}</strong>
<StatusBadge status={item.status} /> <StatusBadge status={item.status} />
{item.scene_id && <span style={{ color: palette.muted }}>场景 #{item.scene_id}</span>} {item.scene_id && <span style={{ color: palette.muted }}>场景 #{item.scene_id}</span>}
{item.processor && <StatusBadge tone="info">{item.processor}</StatusBadge>}
{item.satellite && <span style={{ color: palette.muted }}>{item.satellite}</span>} {item.satellite && <span style={{ color: palette.muted }}>{item.satellite}</span>}
{item.imaging_date && <span style={{ color: palette.muted }}>{formatYmd(item.imaging_date, 'zh')}</span>} {item.imaging_date && <span style={{ color: palette.muted }}>{formatYmd(item.imaging_date, 'zh')}</span>}
{item.polarization && <span style={{ color: palette.subtle }}>{item.polarization}</span>} {item.polarization && <span style={{ color: palette.subtle }}>{item.polarization}</span>}
@@ -331,8 +333,13 @@ function WaterResultRow({ item, onShowMap }) {
<div style={{ display: 'grid', gridTemplateColumns: 'repeat(3, minmax(0, 1fr))', gap: 8, marginTop: 8 }}> <div style={{ display: 'grid', gridTemplateColumns: 'repeat(3, minmax(0, 1fr))', gap: 8, marginTop: 8 }}>
<KeyValue label="面积" value={formatArea(item.water_area_km2)} strong /> <KeyValue label="面积" value={formatArea(item.water_area_km2)} strong />
<KeyValue label="像素" value={item.water_pixel_count?.toLocaleString?.() || '-'} /> <KeyValue label="像素" value={item.water_pixel_count?.toLocaleString?.() || '-'} />
<KeyValue label="阈值" value={item.otsu_threshold_db != null ? Number(item.otsu_threshold_db).toFixed(2) : '-'} /> <KeyValue label="阈值" value={thresholdValue != null ? Number(thresholdValue).toFixed(2) : '-'} />
</div> </div>
{(item.preview_path || item.vector_path) && (
<div style={{ color: palette.subtle, fontSize: 11, marginTop: 6, overflow: 'hidden', textOverflow: 'ellipsis', whiteSpace: 'nowrap' }}>
{[item.preview_path ? `preview: ${item.preview_path}` : null, item.vector_path ? `vector: ${item.vector_path}` : null].filter(Boolean).join(' | ')}
</div>
)}
{item.error_msg && <div style={{ color: palette.red, fontSize: 11, marginTop: 6 }}>{item.error_msg}</div>} {item.error_msg && <div style={{ color: palette.red, fontSize: 11, marginTop: 6 }}>{item.error_msg}</div>}
</div> </div>
<button type="button" style={buttonStyle('secondary', asStatus(item.status) !== 'DONE')} disabled={asStatus(item.status) !== 'DONE'} onClick={() => onShowMap(item)}>显示图层</button> <button type="button" style={buttonStyle('secondary', asStatus(item.status) !== 'DONE')} disabled={asStatus(item.status) !== 'DONE'} onClick={() => onShowMap(item)}>显示图层</button>
@@ -451,6 +458,7 @@ export default function FloodAnalysisWorkspace({
const [waterTotal, setWaterTotal] = useState(0); const [waterTotal, setWaterTotal] = useState(0);
const [waterPage, setWaterPage] = useState(0); const [waterPage, setWaterPage] = useState(0);
const [waterLoading, setWaterLoading] = useState(false); const [waterLoading, setWaterLoading] = useState(false);
const [waterProcessor, setWaterProcessor] = useState('otsu');
const [sourceAoiMode, setSourceAoiMode] = useState('none'); const [sourceAoiMode, setSourceAoiMode] = useState('none');
const [sourceRegionOptions, setSourceRegionOptions] = useState({ provinces: [], cities: [] }); const [sourceRegionOptions, setSourceRegionOptions] = useState({ provinces: [], cities: [] });
@@ -774,7 +782,7 @@ export default function FloodAnalysisWorkspace({
showMessage('info', `正在提交场景 #${scene.id} 的水体提取任务...`); showMessage('info', `正在提交场景 #${scene.id} 的水体提取任务...`);
onTaskStart?.(null, '正在提交水体提取任务...'); onTaskStart?.(null, '正在提交水体提取任务...');
try { 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, '水体提取任务已启动'); if (res.data?.task_id) onTaskStart?.(res.data.task_id, '水体提取任务已启动');
await loadWaterResults(0); await loadWaterResults(0);
showMessage('success', `场景 #${scene.id} 的水体提取任务已提交。`); showMessage('success', `场景 #${scene.id} 的水体提取任务已提交。`);
@@ -1181,6 +1189,13 @@ export default function FloodAnalysisWorkspace({
<section style={sectionStyle}> <section style={sectionStyle}>
<SectionHeader title={`分析就绪场景 (${scenesTotal})`} actions={<button type="button" style={buttonStyle('quiet', scenesLoading)} disabled={scenesLoading} onClick={() => loadScenes(scenesPage)}>刷新场景</button>} /> <SectionHeader title={`分析就绪场景 (${scenesTotal})`} actions={<button type="button" style={buttonStyle('quiet', scenesLoading)} disabled={scenesLoading} onClick={() => loadScenes(scenesPage)}>刷新场景</button>} />
<div style={{ display: 'flex', gap: 8, flexWrap: 'wrap', alignItems: 'center', justifyContent: 'space-between', marginBottom: 8 }}>
<span style={{ color: palette.muted, fontSize: 12 }}>水体提取处理器</span>
<select value={waterProcessor} onChange={event => setWaterProcessor(event.target.value)} style={{ ...inputStyle, width: 180 }}>
<option value="otsu">Otsu GeoTIFF</option>
<option value="gf3_hh_hv">GF3 HH/HV</option>
</select>
</div>
<div style={{ display: 'grid', gap: 8 }}> <div style={{ display: 'grid', gap: 8 }}>
{scenesLoading && <EmptyState>场景加载中...</EmptyState>} {scenesLoading && <EmptyState>场景加载中...</EmptyState>}
{!scenesLoading && scenes.length === 0 && <EmptyState>暂无分析就绪场景</EmptyState>} {!scenesLoading && scenes.length === 0 && <EmptyState>暂无分析就绪场景</EmptyState>}
+3
View File
@@ -18,5 +18,8 @@ export const listAssetIssues = (params = {}) =>
export const unpackSentinel1Source = (assetId, payload = {}) => export const unpackSentinel1Source = (assetId, payload = {}) =>
apiClient.post(`/assets/sources/${assetId}/unpack-sentinel1`, payload).then(r => r.data); apiClient.post(`/assets/sources/${assetId}/unpack-sentinel1`, payload).then(r => r.data);
export const materializeSourceAsset = (assetId, payload = {}) =>
apiClient.post(`/assets/sources/${assetId}/materialize`, payload).then(r => r.data);
export const unpackSentinel1Batch = (payload = {}) => export const unpackSentinel1Batch = (payload = {}) =>
apiClient.post('/assets/inventory/unpack-sentinel1', payload).then(r => r.data); apiClient.post('/assets/inventory/unpack-sentinel1', payload).then(r => r.data);