feat: engineer SBAS timeseries production workflow

This commit is contained in:
2026-04-29 14:43:31 +08:00
parent dace8b20f6
commit 4c0d1f2c2b
54 changed files with 5843 additions and 201 deletions
+279 -22
View File
@@ -2,6 +2,7 @@
from __future__ import annotations
import argparse
import json
from pathlib import Path
import numpy as np
@@ -11,9 +12,12 @@ gdal.UseExceptions()
DEFAULT_WAVELENGTH = 0.23793052222222222
DEFAULT_NODATA = -9999.0
DEFAULT_REFERENCE_MODE = "none"
DEFAULT_REFERENCE_MODE = "coh_median"
DEFAULT_REFERENCE_COH_THRESHOLD = 0.30
DEFAULT_DERAMP_MODE = "plane"
DEFAULT_DERAMP_COH_THRESHOLD = 0.30
REFERENCE_MODE_CHOICES = ("none", "coh_median")
DERAMP_MODE_CHOICES = ("none", "plane")
def parse_args() -> argparse.Namespace:
@@ -54,7 +58,7 @@ def parse_args() -> argparse.Namespace:
type=str,
choices=REFERENCE_MODE_CHOICES,
default=DEFAULT_REFERENCE_MODE,
help="Optional reference normalization mode for debug exports",
help="Reference normalization mode applied before final displacement export",
)
parser.add_argument(
"--reference-coh-threshold",
@@ -62,10 +66,23 @@ def parse_args() -> argparse.Namespace:
default=DEFAULT_REFERENCE_COH_THRESHOLD,
help="Minimum coherence used to select reference pixels for normalization",
)
parser.add_argument(
"--deramp-mode",
type=str,
choices=DERAMP_MODE_CHOICES,
default=DEFAULT_DERAMP_MODE,
help="Optional long-wavelength ramp removal applied after reference normalization",
)
parser.add_argument(
"--deramp-coh-threshold",
type=float,
default=DEFAULT_DERAMP_COH_THRESHOLD,
help="Minimum coherence used when selecting pixels for deramp fitting",
)
parser.add_argument(
"--include-disp-full",
action="store_true",
help="Also export the coherence-unmasked displacement GeoTIFF for debugging",
help="Also export the coherence-unmasked final displacement GeoTIFF",
)
return parser.parse_args()
@@ -93,6 +110,51 @@ def write_geotiff(array: np.ndarray, ref_ds: gdal.Dataset, out_path: Path, nodat
ds = None
def _resolve_phase_source(work_dir: Path) -> dict[str, str | bool]:
ionosphere_phase = work_dir / "ionosphere" / "nondispersive.bil.unwCor.filt.geo.vrt"
ionosphere_mask = work_dir / "ionosphere" / "mask.bil.geo.vrt"
full_unwrap = work_dir / "interferogram" / "filt_topophase.unw.geo.vrt"
if ionosphere_phase.exists():
return {
"phase_path": str(ionosphere_phase),
"phase_source": "ionosphere_nondispersive",
"mask_path": str(ionosphere_mask) if ionosphere_mask.exists() else "",
"ionosphere_corrected": True,
}
return {
"phase_path": str(full_unwrap),
"phase_source": "interferogram_unwrapped",
"mask_path": "",
"ionosphere_corrected": False,
}
def _select_support_mask(
*,
base_mask: np.ndarray,
amp_valid: np.ndarray,
disp_valid: np.ndarray,
coh: np.ndarray,
selection_threshold: float,
) -> tuple[np.ndarray, dict[str, float | int | str]]:
fallback = ""
support_mask = base_mask & (coh >= selection_threshold)
if not support_mask.any():
support_mask = base_mask & (coh > 0)
fallback = "coh>0"
if not support_mask.any():
support_mask = amp_valid & disp_valid
fallback = "amp_only"
stats: dict[str, float | int | str] = {
"selection_threshold": float(selection_threshold),
"fallback": fallback,
"support_ratio": float(base_mask.mean()),
"support_count": int(support_mask.sum()),
"support_mask_ratio": float(support_mask.mean()),
}
return support_mask, stats
def compute_reference_offset(
disp_m_raw: np.ndarray,
amp: np.ndarray,
@@ -100,7 +162,7 @@ def compute_reference_offset(
coh_threshold: float,
reference_mode: str,
reference_coh_threshold: float,
) -> tuple[float, dict[str, float | int | str]]:
) -> tuple[float, np.ndarray, dict[str, float | int | str]]:
amp_valid = np.isfinite(amp) & (amp != 0)
coh_finite = np.isfinite(coh)
disp_valid = np.isfinite(disp_m_raw)
@@ -121,17 +183,16 @@ def compute_reference_offset(
"fallback": "",
}
if normalized_mode == "none":
return 0.0, stats
return 0.0, base_mask, stats
selection_threshold = min(1.0, max(0.0, max(float(coh_threshold), float(reference_coh_threshold))))
reference_mask = base_mask & (coh >= selection_threshold)
fallback = ""
if not reference_mask.any():
reference_mask = base_mask & (coh > 0)
fallback = "coh>0"
if not reference_mask.any():
reference_mask = amp_valid & disp_valid
fallback = "amp_only"
reference_mask, mask_stats = _select_support_mask(
base_mask=base_mask,
amp_valid=amp_valid,
disp_valid=disp_valid,
coh=coh,
selection_threshold=selection_threshold,
)
reference_count = int(reference_mask.sum())
if reference_count <= 0:
@@ -141,11 +202,101 @@ def compute_reference_offset(
{
"reference_count": reference_count,
"reference_ratio": float(reference_mask.mean()),
"selection_threshold": float(selection_threshold),
"fallback": fallback,
"selection_threshold": float(mask_stats["selection_threshold"]),
"fallback": str(mask_stats["fallback"]),
}
)
return float(np.median(disp_m_raw[reference_mask])), stats
return float(np.median(disp_m_raw[reference_mask])), reference_mask, stats
def compute_deramp_surface(
disp_m: np.ndarray,
amp: np.ndarray,
coh: np.ndarray,
coh_threshold: float,
deramp_mode: str,
deramp_coh_threshold: float,
) -> tuple[np.ndarray, np.ndarray, dict[str, float | int | str | bool]]:
amp_valid = np.isfinite(amp) & (amp != 0)
coh_finite = np.isfinite(coh)
disp_valid = np.isfinite(disp_m)
base_mask = amp_valid & coh_finite & disp_valid
normalized_mode = str(deramp_mode or DEFAULT_DERAMP_MODE).strip().lower()
if normalized_mode not in DERAMP_MODE_CHOICES:
raise ValueError(f"Unsupported deramp mode: {deramp_mode}")
empty_surface = np.zeros_like(disp_m, dtype=np.float32)
stats: dict[str, float | int | str | bool] = {
"mode": normalized_mode,
"applied": False,
"fit_count": 0,
"fit_ratio": 0.0,
"selection_threshold": 0.0,
"fallback": "",
"sample_step": 0,
"sample_count": 0,
}
if normalized_mode == "none":
return empty_surface, base_mask, stats
if not base_mask.any():
return empty_surface, base_mask, stats
selection_threshold = min(1.0, max(0.0, max(float(coh_threshold), float(deramp_coh_threshold))))
fit_mask, mask_stats = _select_support_mask(
base_mask=base_mask,
amp_valid=amp_valid,
disp_valid=disp_valid,
coh=coh,
selection_threshold=selection_threshold,
)
fit_count = int(fit_mask.sum())
stats.update(
{
"fit_count": fit_count,
"fit_ratio": float(fit_mask.mean()),
"selection_threshold": float(mask_stats["selection_threshold"]),
"fallback": str(mask_stats["fallback"]),
}
)
if fit_count < 3:
stats["fallback"] = "insufficient_support"
return empty_surface, fit_mask, stats
yy, xx = np.indices(disp_m.shape, dtype=np.float64)
xs = xx[fit_mask]
ys = yy[fit_mask]
zs = disp_m[fit_mask].astype(np.float64)
sample_step = max(1, fit_count // 250_000)
if sample_step > 1:
xs = xs[::sample_step]
ys = ys[::sample_step]
zs = zs[::sample_step]
sample_count = int(zs.size)
stats["sample_step"] = int(sample_step)
stats["sample_count"] = sample_count
if sample_count < 3:
stats["fallback"] = "insufficient_sample"
return empty_surface, fit_mask, stats
design = np.column_stack([xs, ys, np.ones_like(xs)])
coeffs, _, _, _ = np.linalg.lstsq(design, zs, rcond=None)
plane = (
coeffs[0] * xx
+ coeffs[1] * yy
+ coeffs[2]
).astype(np.float32)
stats.update(
{
"applied": True,
"coef_x_per_pixel": float(coeffs[0]),
"coef_y_per_pixel": float(coeffs[1]),
"intercept_m": float(coeffs[2]),
"left_right_delta_m": float(coeffs[0] * max(disp_m.shape[1] - 1, 0)),
"top_bottom_delta_m": float(coeffs[1] * max(disp_m.shape[0] - 1, 0)),
}
)
return plane, fit_mask, stats
def export_products(
@@ -156,32 +307,48 @@ def export_products(
coh_threshold: float,
reference_mode: str = DEFAULT_REFERENCE_MODE,
reference_coh_threshold: float = DEFAULT_REFERENCE_COH_THRESHOLD,
deramp_mode: str = DEFAULT_DERAMP_MODE,
deramp_coh_threshold: float = DEFAULT_DERAMP_COH_THRESHOLD,
include_disp_full: bool = False,
nodata: float = DEFAULT_NODATA,
) -> dict[str, Path]:
unw_path = work_dir / "interferogram" / "filt_topophase.unw.geo.vrt"
cor_path = work_dir / "interferogram" / "topophase.cor.geo.vrt"
phase_source = _resolve_phase_source(work_dir)
phase_path = Path(str(phase_source["phase_path"]))
mask_path = Path(str(phase_source["mask_path"])) if str(phase_source["mask_path"]) else None
if not unw_path.exists():
raise FileNotFoundError(f"Missing unwrapped product: {unw_path}")
if not cor_path.exists():
raise FileNotFoundError(f"Missing coherence product: {cor_path}")
if not phase_path.exists():
raise FileNotFoundError(f"Missing phase source product: {phase_path}")
unw_ds = gdal.Open(str(unw_path))
cor_ds = gdal.Open(str(cor_path))
if unw_ds is None or cor_ds is None:
phase_ds = gdal.Open(str(phase_path))
mask_ds = gdal.Open(str(mask_path)) if mask_path is not None else None
if unw_ds is None or cor_ds is None or phase_ds is None:
raise RuntimeError("Failed to open ISCE2 geo products with GDAL.")
amp = unw_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
phase = unw_ds.GetRasterBand(2).ReadAsArray().astype(np.float32)
if bool(phase_source["ionosphere_corrected"]):
phase = phase_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
else:
phase = unw_ds.GetRasterBand(2).ReadAsArray().astype(np.float32)
coh_band = 2 if cor_ds.RasterCount >= 2 else 1
coh = cor_ds.GetRasterBand(coh_band).ReadAsArray().astype(np.float32)
coh_valid = np.isfinite(coh) & (coh > 0)
amp_valid = np.isfinite(amp) & (amp != 0)
ionosphere_mask_valid = None
if mask_ds is not None:
ionosphere_mask = mask_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
ionosphere_mask_valid = np.isfinite(ionosphere_mask) & (ionosphere_mask > 0)
disp_m_raw = phase * wavelength / (4.0 * np.pi)
reference_offset_m, reference_stats = compute_reference_offset(
reference_offset_m, reference_mask, reference_stats = compute_reference_offset(
disp_m_raw=disp_m_raw,
amp=amp,
coh=coh,
@@ -189,10 +356,25 @@ def export_products(
reference_mode=reference_mode,
reference_coh_threshold=reference_coh_threshold,
)
disp_m = disp_m_raw - reference_offset_m
disp_m_ref = disp_m_raw - reference_offset_m
deramp_surface_m, deramp_mask, deramp_stats = compute_deramp_surface(
disp_m=disp_m_ref,
amp=amp,
coh=coh,
coh_threshold=coh_threshold,
deramp_mode=deramp_mode,
deramp_coh_threshold=deramp_coh_threshold,
)
disp_m = disp_m_ref - deramp_surface_m
disp_m_full = disp_m.copy()
mask = (~amp_valid) | (~np.isfinite(disp_m)) | (~np.isfinite(coh)) | (coh < coh_threshold)
if ionosphere_mask_valid is not None:
mask |= ~ionosphere_mask_valid
disp_m_raw_masked = disp_m_raw.copy()
disp_m_raw_masked[mask] = nodata
disp_m_ref_masked = disp_m_ref.copy()
disp_m_ref_masked[mask] = nodata
disp_m_masked = disp_m.copy()
disp_m_masked[mask] = nodata
disp_m_full[(~amp_valid) | (~np.isfinite(disp_m_full))] = nodata
@@ -202,8 +384,13 @@ def export_products(
output_dir.mkdir(parents=True, exist_ok=True)
out_disp = output_dir / f"{prefix}_disp.tif"
out_disp_raw = output_dir / f"{prefix}_disp_raw.tif"
out_disp_ref = output_dir / f"{prefix}_disp_ref.tif"
out_coh = output_dir / f"{prefix}_coh.tif"
out_meta = output_dir / f"{prefix}_disp_meta.json"
write_geotiff(disp_m_raw_masked, unw_ds, out_disp_raw, nodata)
write_geotiff(disp_m_ref_masked, unw_ds, out_disp_ref, nodata)
write_geotiff(disp_m_masked, unw_ds, out_disp, nodata)
write_geotiff(coh_out, cor_ds, out_coh, nodata)
out_disp_full = None
@@ -211,14 +398,50 @@ def export_products(
out_disp_full = output_dir / f"{prefix}_disp_full.tif"
write_geotiff(disp_m_full, unw_ds, out_disp_full, nodata)
valid_raw_masked = disp_m_raw_masked[disp_m_raw_masked != nodata]
valid_ref_masked = disp_m_ref_masked[disp_m_ref_masked != nodata]
valid_disp = disp_m_masked[disp_m_masked != nodata]
valid_coh = coh_out[coh_out != nodata]
valid_full = disp_m_full[disp_m_full != nodata] if include_disp_full else np.array([], dtype=np.float32)
valid_raw = disp_m_raw[amp_valid & np.isfinite(disp_m_raw)]
using_reference = str(reference_stats["mode"]) != "none"
using_deramp = bool(deramp_stats["applied"])
meta_payload = {
"work_dir": str(work_dir),
"output_dir": str(output_dir),
"prefix": prefix,
"coh_threshold": float(coh_threshold),
"phase_source": {
"kind": str(phase_source["phase_source"]),
"path": str(phase_path),
"ionosphere_corrected": bool(phase_source["ionosphere_corrected"]),
"mask_path": str(mask_path) if mask_path is not None else "",
"mask_applied": bool(ionosphere_mask_valid is not None),
"mask_valid_ratio": float(ionosphere_mask_valid.mean()) if ionosphere_mask_valid is not None else None,
},
"reference": {
**reference_stats,
"offset_m": float(reference_offset_m),
"support_count": int(reference_mask.sum()),
},
"deramp": {
**deramp_stats,
"support_count": int(deramp_mask.sum()),
},
"ranges_m": {
"raw_valid": [float(valid_raw.min()), float(valid_raw.max())] if valid_raw.size else [],
"raw_masked": [float(valid_raw_masked.min()), float(valid_raw_masked.max())] if valid_raw_masked.size else [],
"ref_masked": [float(valid_ref_masked.min()), float(valid_ref_masked.max())] if valid_ref_masked.size else [],
"final_masked": [float(valid_disp.min()), float(valid_disp.max())] if valid_disp.size else [],
"final_full": [float(valid_full.min()), float(valid_full.max())] if valid_full.size else [],
},
}
out_meta.write_text(json.dumps(meta_payload, indent=2, ensure_ascii=False), encoding="utf-8")
print(f"Work dir: {work_dir}")
print(f"Output prefix: {prefix}")
print(f"Phase source: {phase_source['phase_source']}")
print(f"Coherence threshold: {coh_threshold}")
print(f"Reference mode: {reference_stats['mode']}")
if using_reference:
@@ -233,29 +456,61 @@ def export_products(
print(f"Reference offset: {reference_offset_m:.4f} m")
if reference_stats["fallback"]:
print(f"Reference fallback: {reference_stats['fallback']}")
print(f"Deramp mode: {deramp_stats['mode']}")
if using_deramp:
print(
"Deramp coh floor: "
f"{float(deramp_stats['selection_threshold']):.2f}"
)
print(
"Deramp pixel ratio: "
f"{float(deramp_stats['fit_ratio'])*100:.2f}%"
)
print(
"Deramp plane delta: "
f"dx={float(deramp_stats['left_right_delta_m']):.4f} m, "
f"dy={float(deramp_stats['top_bottom_delta_m']):.4f} m"
)
if deramp_stats["fallback"]:
print(f"Deramp fallback: {deramp_stats['fallback']}")
elif deramp_stats["fallback"]:
print(f"Deramp fallback: {deramp_stats['fallback']}")
print(f"Unwrap support ratio: {amp_valid.mean()*100:.2f}%")
print(f"Coherence support ratio: {coh_valid.mean()*100:.2f}%")
print(f"Masked disp ratio: {(disp_m_masked != nodata).mean()*100:.2f}%")
if valid_raw.size:
print(f"Raw disp range: [{valid_raw.min():.4f}, {valid_raw.max():.4f}] m")
if valid_raw_masked.size:
print(f"Raw masked range: [{valid_raw_masked.min():.4f}, {valid_raw_masked.max():.4f}] m")
if valid_ref_masked.size:
label = "Ref disp range" if using_reference else "Ref disp range"
print(f"{label + ':':24}[{valid_ref_masked.min():.4f}, {valid_ref_masked.max():.4f}] m")
if valid_disp.size:
label = "Norm disp range" if using_reference else "Disp range"
label = "Final disp range" if using_reference or using_deramp else "Disp range"
print(f"{label + ':':24}[{valid_disp.min():.4f}, {valid_disp.max():.4f}] m")
if include_disp_full and valid_full.size:
label = "Norm full disp range" if using_reference else "Full disp range"
label = "Final full disp range" if using_reference or using_deramp else "Full disp range"
print(f"{label + ':':24}[{valid_full.min():.4f}, {valid_full.max():.4f}] m")
if valid_coh.size:
print(f"Coherence range: [{valid_coh.min():.4f}, {valid_coh.max():.4f}]")
print(f"Wrote: {out_disp_raw}")
print(f"Wrote: {out_disp_ref}")
print(f"Wrote: {out_disp}")
if out_disp_full is not None:
print(f"Wrote: {out_disp_full}")
print(f"Wrote: {out_coh}")
print(f"Wrote: {out_meta}")
unw_ds = None
cor_ds = None
phase_ds = None
mask_ds = None
outputs: dict[str, Path] = {
"disp_raw": out_disp_raw,
"disp_ref": out_disp_ref,
"disp": out_disp,
"coh": out_coh,
"meta": out_meta,
}
if out_disp_full is not None:
outputs["disp_full"] = out_disp_full
@@ -276,6 +531,8 @@ def main() -> int:
coh_threshold=args.coh_threshold,
reference_mode=args.reference_mode,
reference_coh_threshold=args.reference_coh_threshold,
deramp_mode=args.deramp_mode,
deramp_coh_threshold=args.deramp_coh_threshold,
include_disp_full=args.include_disp_full,
)
return 0
@@ -2,9 +2,11 @@
from __future__ import annotations
import sys
import re
import xml.etree.ElementTree as ET
from dataclasses import dataclass
from pathlib import Path
from typing import Callable, Iterable, Mapping, Optional
from typing import Any, Callable, Iterable, Mapping, Optional
try:
from .convert_lt1_orbit_to_isce_xml import (
@@ -37,6 +39,7 @@ DEFAULT_WSL_DEM_CANDIDATES = (
"/mnt/d/SRTM30m/SRTMDEM_RSP_SARscape",
)
DEFAULT_WINDOWS_ORBIT_POOL_CANDIDATES = (r"D:\orbit_pools\isce2",)
DEM_SIDECAR_PROPERTY_NAMES = ("file_name", "metadata_location", "extra_file_name")
@dataclass(frozen=True)
@@ -126,6 +129,24 @@ def resolve_prepared_dem_path(
return resolve_existing_prepared_file(candidates, path_transform=path_transform)
def repair_related_dem_sidecars(
dem_path: Path,
*,
write_changes: bool = True,
) -> list[dict[str, Any]]:
reports: list[dict[str, Any]] = []
seen: set[str] = set()
for candidate in _related_dem_sidecar_candidates(dem_path):
key = str(candidate)
if key in seen:
continue
seen.add(key)
report = repair_dem_sidecar_paths(candidate, write_changes=write_changes)
if report.get("exists"):
reports.append(report)
return reports
def resolve_existing_directory(
candidates: Iterable[str | Path],
path_transform: PathTransform = identity_path_transform,
@@ -162,6 +183,63 @@ def resolve_existing_prepared_file(
return None
def repair_dem_sidecar_paths(
dem_path: Path,
*,
write_changes: bool = True,
) -> dict[str, Any]:
normalized_path = Path(str(dem_path))
xml_path = Path(str(normalized_path) + ".xml")
vrt_path = Path(str(normalized_path) + ".vrt")
report: dict[str, Any] = {
"dem_path": str(normalized_path),
"xml_path": str(xml_path),
"vrt_path": str(vrt_path),
"exists": xml_path.exists(),
"changed": False,
"updated_fields": [],
"expected": {},
"current": {},
}
if not xml_path.exists():
return report
expected_values = {
"file_name": _to_isce_sidecar_path(normalized_path),
"metadata_location": _to_isce_sidecar_path(xml_path),
"extra_file_name": _to_isce_sidecar_path(vrt_path) if vrt_path.exists() else "",
}
tree = ET.parse(xml_path)
root = tree.getroot()
updates: list[str] = []
for prop in root.findall("property"):
name = str(prop.get("name") or "").strip()
if name not in DEM_SIDECAR_PROPERTY_NAMES:
continue
value_node = prop.find("value")
if value_node is None:
value_node = ET.SubElement(prop, "value")
current_value = str(value_node.text or "").strip()
expected_value = expected_values.get(name, "")
report["current"][name] = current_value
report["expected"][name] = expected_value
if not expected_value:
continue
if current_value == expected_value:
continue
value_node.text = expected_value
updates.append(name)
if updates and write_changes:
ET.indent(tree, space=" ")
tree.write(xml_path, encoding="utf-8")
report["changed"] = bool(updates)
report["updated_fields"] = updates
return report
def ensure_lt1_orbit_xml(
date_yyyymmdd: str,
satellite: str,
@@ -251,3 +329,24 @@ def _prepared_dem_variants(value: str | Path) -> tuple[str | Path, ...]:
if text.lower().endswith(".wgs84"):
return (value,)
return (f"{text}.wgs84", value)
def _related_dem_sidecar_candidates(dem_path: Path) -> tuple[Path, ...]:
text = str(dem_path).strip()
if not text:
return ()
if text.lower().endswith(".wgs84"):
raw_path = Path(text[:-6])
return (dem_path, raw_path)
prepared_path = Path(text + ".wgs84")
return (dem_path, prepared_path)
def _to_isce_sidecar_path(path: Path) -> str:
text = str(path).strip()
match = re.match(r"^([A-Za-z]):[\\/](.*)$", text)
if match:
drive = match.group(1).lower()
rest = match.group(2).replace("\\", "/")
return f"/mnt/{drive}/{rest}"
return Path(text).as_posix()
@@ -0,0 +1,91 @@
#!/usr/bin/env python3
from __future__ import annotations
import argparse
import sys
from pathlib import Path
try:
from .lt1_input_resolver import repair_dem_sidecar_paths
except ImportError:
SCRIPT_DIR = Path(__file__).resolve().parent
if str(SCRIPT_DIR) not in sys.path:
sys.path.insert(0, str(SCRIPT_DIR))
from lt1_input_resolver import repair_dem_sidecar_paths # type: ignore
def parse_args() -> argparse.Namespace:
parser = argparse.ArgumentParser(
description="Audit and optionally repair moved ISCE DEM XML sidecars."
)
parser.add_argument(
"--root",
type=Path,
required=True,
help="Directory containing DEM files and sidecars",
)
parser.add_argument(
"--repair",
action="store_true",
help="Write repaired file_name / metadata_location / extra_file_name values back to XML",
)
return parser.parse_args()
def iter_dem_sidecars(root: Path) -> list[Path]:
sidecars: list[Path] = []
for xml_path in sorted(root.rglob("*.xml")):
if xml_path.name.lower().endswith(".aux.xml"):
continue
dem_path = Path(str(xml_path)[:-4])
if dem_path.exists():
sidecars.append(dem_path)
return sidecars
def main() -> int:
args = parse_args()
root = args.root.resolve()
if not root.exists() or not root.is_dir():
raise FileNotFoundError(f"DEM root directory not found: {root}")
sidecars = iter_dem_sidecars(root)
changed_count = 0
mismatch_count = 0
print(f"DEM root: {root}")
print(f"Sidecars: {len(sidecars)}")
for dem_path in sidecars:
report = repair_dem_sidecar_paths(dem_path, write_changes=bool(args.repair))
updated_fields = list(report.get("updated_fields") or [])
if updated_fields:
changed_count += 1
mismatch_count += 1
print(
f"[fixed] {report['xml_path']} -> {', '.join(updated_fields)}"
if args.repair
else f"[mismatch] {report['xml_path']} -> {', '.join(updated_fields)}"
)
continue
current = report.get("current") or {}
expected = report.get("expected") or {}
mismatched = [
key
for key, expected_value in expected.items()
if expected_value and str(current.get(key) or "").strip() != str(expected_value).strip()
]
if mismatched:
mismatch_count += 1
print(f"[mismatch] {report['xml_path']} -> {', '.join(mismatched)}")
else:
print(f"[ok] {report['xml_path']}")
print(f"Mismatched: {mismatch_count}")
if args.repair:
print(f"Repaired: {changed_count}")
return 0
if __name__ == "__main__":
raise SystemExit(main())
@@ -3,6 +3,7 @@ from __future__ import annotations
import argparse
import ast
import importlib.util
import os
import re
import shutil
@@ -14,15 +15,19 @@ from dataclasses import dataclass
from pathlib import Path
from export_isce_geotiff import (
DEFAULT_DERAMP_COH_THRESHOLD,
DEFAULT_DERAMP_MODE,
DEFAULT_REFERENCE_COH_THRESHOLD,
DEFAULT_REFERENCE_MODE,
DEFAULT_WAVELENGTH,
DERAMP_MODE_CHOICES,
REFERENCE_MODE_CHOICES,
export_products,
)
from lt1_input_resolver import (
DEFAULT_WSL_DEM_CANDIDATES,
ensure_lt1_orbit_xml,
repair_related_dem_sidecars,
resolve_prepared_dem_path,
)
@@ -35,7 +40,22 @@ RESUME_STAGE_CHOICES = PIPELINE_STAGE_ORDER[1:]
DEFAULT_EXPORT_GEOCODE_PRODUCTS = [
"interferogram/filt_topophase.unw",
"interferogram/topophase.cor",
"ionosphere/dispersive.bil.unwCor.filt",
"ionosphere/nondispersive.bil.unwCor.filt",
"ionosphere/mask.bil",
]
DEFAULT_EXPORT_GEOCODE_PRODUCTS_NO_IONO = [
"interferogram/filt_topophase.unw",
"interferogram/topophase.cor",
]
DEFAULT_RUBBER_SHEET_SNR_THRESHOLD = 5.0
DEFAULT_RUBBER_SHEET_FILTER_SIZE = 9
DEFAULT_DENSE_WINDOW_WIDTH = 64
DEFAULT_DENSE_WINDOW_HEIGHT = 64
DEFAULT_DENSE_SEARCH_WIDTH = 20
DEFAULT_DENSE_SEARCH_HEIGHT = 20
DEFAULT_DENSE_SKIP_WIDTH = 32
DEFAULT_DENSE_SKIP_HEIGHT = 32
@dataclass
@@ -59,6 +79,18 @@ class PipelineConfig:
target_grid_size_m: int
geo_posting_deg: float
geocode_products: list[str] | None
ionosphere_correction: bool
dense_offsets: bool
rubbersheet_range: bool
rubbersheet_azimuth: bool
rubber_sheet_snr_threshold: float
rubber_sheet_filter_size: int
dense_window_width: int
dense_window_height: int
dense_search_width: int
dense_search_height: int
dense_skip_width: int
dense_skip_height: int
def parse_args() -> argparse.Namespace:
@@ -66,7 +98,7 @@ def parse_args() -> argparse.Namespace:
repo_root = script_dir.parent
parser = argparse.ArgumentParser(
description="Run an LT-1 ISCE2 DInSAR production pipeline with SNAPHU."
description="Run an LT-1 ISCE2 DInSAR production pipeline with the standard stripmap workflow."
)
parser.add_argument("task_dir", help="Task directory, for example Task_20250112_20250309")
parser.add_argument(
@@ -156,7 +188,7 @@ def parse_args() -> argparse.Namespace:
"--reference-mode",
choices=REFERENCE_MODE_CHOICES,
default=DEFAULT_REFERENCE_MODE,
help="Optional reference normalization mode used only for debug exports",
help="Reference normalization mode applied during final displacement export",
)
parser.add_argument(
"--reference-coh-threshold",
@@ -164,6 +196,18 @@ def parse_args() -> argparse.Namespace:
default=DEFAULT_REFERENCE_COH_THRESHOLD,
help="Minimum coherence used when selecting reference pixels for export normalization",
)
parser.add_argument(
"--deramp-mode",
choices=DERAMP_MODE_CHOICES,
default=DEFAULT_DERAMP_MODE,
help="Optional ramp-removal mode applied after reference normalization",
)
parser.add_argument(
"--deramp-coh-threshold",
type=float,
default=DEFAULT_DERAMP_COH_THRESHOLD,
help="Minimum coherence used when selecting pixels for deramp fitting",
)
parser.add_argument(
"--target-grid-size-m",
type=int,
@@ -180,6 +224,76 @@ def parse_args() -> argparse.Namespace:
action="store_true",
help="Let ISCE2 geocode its full default product list instead of the reduced export-only list.",
)
parser.add_argument(
"--no-ionosphere-correction",
action="store_false",
dest="ionosphere_correction",
help="Disable split-spectrum dispersive correction and export the standard unwrapped interferogram.",
)
parser.set_defaults(ionosphere_correction=True)
parser.add_argument(
"--dense-offsets",
action="store_true",
help="Enable ISCE2 dense offset estimation before fine resampling.",
)
parser.add_argument(
"--rubbersheet-range",
action="store_true",
help="Enable ISCE2 range rubbersheeting using dense offsets.",
)
parser.add_argument(
"--rubbersheet-azimuth",
action="store_true",
help="Enable ISCE2 azimuth rubbersheeting using dense offsets.",
)
parser.add_argument(
"--rubber-sheet-snr-threshold",
type=float,
default=DEFAULT_RUBBER_SHEET_SNR_THRESHOLD,
help="SNR threshold used by ISCE2 rubbersheet offset masking.",
)
parser.add_argument(
"--rubber-sheet-filter-size",
type=int,
default=DEFAULT_RUBBER_SHEET_FILTER_SIZE,
help="Median filter size used by ISCE2 rubbersheet offset masking.",
)
parser.add_argument(
"--dense-window-width",
type=int,
default=DEFAULT_DENSE_WINDOW_WIDTH,
help="Dense offset correlation window width.",
)
parser.add_argument(
"--dense-window-height",
type=int,
default=DEFAULT_DENSE_WINDOW_HEIGHT,
help="Dense offset correlation window height.",
)
parser.add_argument(
"--dense-search-width",
type=int,
default=DEFAULT_DENSE_SEARCH_WIDTH,
help="Dense offset search window width.",
)
parser.add_argument(
"--dense-search-height",
type=int,
default=DEFAULT_DENSE_SEARCH_HEIGHT,
help="Dense offset search window height.",
)
parser.add_argument(
"--dense-skip-width",
type=int,
default=DEFAULT_DENSE_SKIP_WIDTH,
help="Dense offset sampling stride in range direction.",
)
parser.add_argument(
"--dense-skip-height",
type=int,
default=DEFAULT_DENSE_SKIP_HEIGHT,
help="Dense offset sampling stride in azimuth direction.",
)
parser.add_argument(
"--resume-from",
choices=RESUME_STAGE_CHOICES,
@@ -219,11 +333,75 @@ def parse_args() -> argparse.Namespace:
raise ValueError("--target-grid-size-m must be greater than 0")
if args.reference_coh_threshold < 0 or args.reference_coh_threshold > 1:
raise ValueError("--reference-coh-threshold must be between 0 and 1")
if args.deramp_coh_threshold < 0 or args.deramp_coh_threshold > 1:
raise ValueError("--deramp-coh-threshold must be between 0 and 1")
if args.rubber_sheet_snr_threshold < 0:
raise ValueError("--rubber-sheet-snr-threshold must be non-negative")
if args.rubber_sheet_filter_size <= 0:
raise ValueError("--rubber-sheet-filter-size must be greater than 0")
for field_name in (
"dense_window_width",
"dense_window_height",
"dense_search_width",
"dense_search_height",
"dense_skip_width",
"dense_skip_height",
):
if int(getattr(args, field_name)) <= 0:
raise ValueError(f"--{field_name.replace('_', '-')} must be greater than 0")
if args.force and args.resume_from:
raise ValueError("--force cannot be used together with --resume-from")
return args
def _find_python_module(module_name: str) -> bool:
try:
return importlib.util.find_spec(module_name) is not None
except ModuleNotFoundError:
return False
def validate_runtime_dependencies(args: argparse.Namespace) -> str:
errors: list[str] = []
env_for_cli = build_process_env()
ionosphere_correction = bool(getattr(args, "ionosphere_correction", True))
if ionosphere_correction:
missing_ionosphere_modules: list[str] = []
if not _find_python_module("cv2"):
missing_ionosphere_modules.append("cv2")
if not _find_python_module("scipy"):
missing_ionosphere_modules.append("scipy")
if missing_ionosphere_modules:
errors.append(
"Missing Python dependencies for the ISCE2 stripmap ionosphere step: "
+ ", ".join(missing_ionosphere_modules)
+ ". The managed LT-1 workflow enables split-spectrum dispersive correction "
"before geocode. Install the missing packages in the WSL runtime, for example: "
"conda install -n insar_wsl_v1 -c conda-forge opencv scipy."
)
if (args.rubbersheet_range or args.rubbersheet_azimuth) and not _find_python_module(
"astropy.convolution"
):
errors.append(
"Missing Python dependency 'astropy.convolution'. "
"ISCE2 stripmap rubbersheeting imports astropy.convolution in "
"runRubbersheetRange.py. Install astropy in the WSL runtime, for example: "
"conda install -n insar_wsl_v1 -c conda-forge astropy."
)
if ionosphere_correction and not shutil.which("imageMath.py", path=str(env_for_cli.get("PATH") or "")):
errors.append(
"Missing CLI dependency 'imageMath.py' on PATH. "
"ISCE2 stripmap shells out to imageMath.py in the ionosphere step, so a missing PATH entry "
"will only surface late in the run. Export the active conda env bin directory into PATH "
"before launching production."
)
return "\n".join(errors)
def locate_stripmap_app() -> Path:
import isce
@@ -233,6 +411,29 @@ def locate_stripmap_app() -> Path:
return app_path
def locate_isce_applications_dir() -> Path | None:
spec = importlib.util.find_spec("isce")
if not spec or not spec.origin:
return None
app_dir = Path(spec.origin).resolve().parent / "applications"
if app_dir.exists():
return app_dir
return None
def build_process_env(base_env: dict[str, str] | None = None) -> dict[str, str]:
env = dict(base_env or os.environ.copy())
path_prefixes = [Path(sys.executable).resolve().parent.as_posix()]
app_dir = locate_isce_applications_dir()
if app_dir:
path_prefixes.append(app_dir.as_posix())
current_path = str(env.get("PATH") or "")
env["PATH"] = ":".join(path_prefixes + ([current_path] if current_path else []))
return env
def normalize_linux_path(value: str | Path) -> Path:
text = str(value).strip()
if text.startswith("\\\\"):
@@ -374,6 +575,14 @@ def resolve_dem(dem_value: str | None) -> Path:
path_transform=normalize_linux_path,
)
if dem_path is not None:
repair_reports = repair_related_dem_sidecars(dem_path)
for report in repair_reports:
if not report.get("changed"):
continue
print(
"Repaired DEM sidecar paths: "
f"{report['xml_path']} -> {', '.join(report['updated_fields'])}"
)
return dem_path
searched = ", ".join(str(path) for path in DEFAULT_WSL_DEM_CANDIDATES)
@@ -502,9 +711,30 @@ def meters_to_geoposting_degrees(target_grid_size_m: int) -> float:
return float(target_grid_size_m) / METERS_PER_DEGREE
def build_default_geocode_products(*, ionosphere_correction: bool) -> list[str]:
return list(
DEFAULT_EXPORT_GEOCODE_PRODUCTS
if ionosphere_correction
else DEFAULT_EXPORT_GEOCODE_PRODUCTS_NO_IONO
)
def write_stripmap_xml(xml_path: Path, config: PipelineConfig) -> None:
bbox_xml = render_bbox(config.bbox)
geocode_list_xml = render_string_list("geocode list", config.geocode_products)
enhancement_props = (
f" <property name=\"do denseoffsets\">{str(config.dense_offsets)}</property>\n"
f" <property name=\"do rubbersheetingRange\">{str(config.rubbersheet_range)}</property>\n"
f" <property name=\"do rubbersheetingAzimuth\">{str(config.rubbersheet_azimuth)}</property>\n"
f" <property name=\"rubber sheet SNR Threshold\">{config.rubber_sheet_snr_threshold}</property>\n"
f" <property name=\"rubber sheet filter size\">{config.rubber_sheet_filter_size}</property>\n"
f" <property name=\"dense window width\">{config.dense_window_width}</property>\n"
f" <property name=\"dense window height\">{config.dense_window_height}</property>\n"
f" <property name=\"dense search width\">{config.dense_search_width}</property>\n"
f" <property name=\"dense search height\">{config.dense_search_height}</property>\n"
f" <property name=\"dense skip width\">{config.dense_skip_width}</property>\n"
f" <property name=\"dense skip height\">{config.dense_skip_height}</property>\n"
)
text = (
"<stripmapApp>\n"
" <component name=\"stripmapApp\">\n"
@@ -514,10 +744,13 @@ def write_stripmap_xml(xml_path: Path, config: PipelineConfig) -> None:
" <property name=\"renderer\">xml</property>\n"
" <property name=\"do unwrap\">True</property>\n"
" <property name=\"unwrapper name\">snaphu</property>\n"
f" <property name=\"do split spectrum\">{str(config.ionosphere_correction)}</property>\n"
f" <property name=\"do dispersive\">{str(config.ionosphere_correction)}</property>\n"
f" <property name=\"posting\">{config.target_grid_size_m}</property>\n"
f" <property name=\"geoPosting\">{config.geo_posting_deg:.12f}</property>\n"
f"{bbox_xml}"
f"{geocode_list_xml}"
f"{enhancement_props}"
f" <property name=\"demFilename\">{config.dem_path.as_posix()}</property>\n"
"\n"
" <component name=\"Reference\">\n"
@@ -553,7 +786,7 @@ def run_logged(stage_name: str, cmd: list[str], cwd: Path, log_path: Path) -> No
handle.write(f"Log: {log_path}\n")
handle.flush()
child_env = os.environ.copy()
child_env = build_process_env()
child_env["PYTHONUNBUFFERED"] = "1"
proc = subprocess.Popen(
cmd,
@@ -730,59 +963,57 @@ def should_run_stage(start_stage: str, stage_name: str) -> bool:
return stage_index >= start_index
def prepare_snaphu_resume(work_dir: Path, bbox: list[float] | None) -> None:
def has_pickle_state(work_dir: Path, state_name: str) -> bool:
pickle_dir = work_dir / "PICKLE"
src = pickle_dir / "filter"
src_xml = pickle_dir / "filter.xml"
dst = pickle_dir / "filter_high_band"
dst_xml = pickle_dir / "filter_high_band.xml"
if not src.exists() or not src_xml.exists():
raise FileNotFoundError("filter step output is missing; cannot prepare SNAPHU resume state.")
shutil.copy2(src, dst)
shutil.copy2(src_xml, dst_xml)
root = ET.fromstring(dst_xml.read_text(encoding="utf-8"))
props = {prop.attrib.get("name"): prop for prop in root.findall("property")}
required = {
"referenceslccroppedproduct": "reference_slc.xml",
"secondaryslccroppedproduct": "secondary_slc.xml",
"referenceslcproduct": "reference_slc.xml",
"secondaryslcproduct": "secondary_slc.xml",
"referencegeometrysystem": "Zero Doppler",
"secondarygeometrysystem": "Zero Doppler",
}
if bbox is not None:
required["estimatedboundingbox"] = str(bbox)
for name, value in required.items():
if name in props:
node = props[name].find("value")
if node is None:
node = ET.SubElement(props[name], "value")
node.text = value
continue
prop = ET.SubElement(root, "property", {"name": name})
ET.SubElement(prop, "value").text = value
dst_xml.write_text(ET.tostring(root, encoding="unicode"), encoding="utf-8")
return (pickle_dir / state_name).exists() and (pickle_dir / f"{state_name}.xml").exists()
def prepare_geocode_resume(work_dir: Path) -> None:
pickle_dir = work_dir / "PICKLE"
unwrap = pickle_dir / "unwrap"
unwrap_xml = pickle_dir / "unwrap.xml"
ionosphere = pickle_dir / "ionosphere"
ionosphere_xml = pickle_dir / "ionosphere.xml"
def resolve_unwrap_start_step(work_dir: Path, *, ionosphere_correction: bool) -> str:
if ionosphere_correction:
if has_pickle_state(work_dir, "ionosphere"):
return "ionosphere"
if has_pickle_state(work_dir, "unwrap_low_band") and has_pickle_state(
work_dir, "unwrap_high_band"
):
return "ionosphere"
if has_pickle_state(work_dir, "filter_low_band") and has_pickle_state(
work_dir, "filter_high_band"
):
return "unwrap"
if has_pickle_state(work_dir, "filter"):
return "filter_low_band"
raise FileNotFoundError(
"Unable to resume the ISCE2 unwrap/ionosphere stage. Missing PICKLE state for "
"filter, filter_low_band/filter_high_band, unwrap_low_band/unwrap_high_band, or ionosphere."
)
if not unwrap.exists() or not unwrap_xml.exists():
raise FileNotFoundError("unwrap step output is missing; cannot prepare geocode resume state.")
if has_pickle_state(work_dir, "unwrap"):
return "unwrap"
if has_pickle_state(work_dir, "filter"):
return "unwrap"
raise FileNotFoundError(
"Unable to resume the ISCE2 unwrap stage. Missing PICKLE state for filter or unwrap."
)
shutil.copy2(unwrap, ionosphere)
shutil.copy2(unwrap_xml, ionosphere_xml)
def resolve_geocode_start_step(work_dir: Path, *, ionosphere_correction: bool) -> str:
if ionosphere_correction:
if has_pickle_state(work_dir, "ionosphere"):
return "geocode"
if has_pickle_state(work_dir, "unwrap_low_band") and has_pickle_state(
work_dir, "unwrap_high_band"
):
return "ionosphere"
raise FileNotFoundError(
"Unable to resume the ISCE2 geocode stage. Missing PICKLE state for ionosphere or "
"unwrap_low_band/unwrap_high_band."
)
if has_pickle_state(work_dir, "unwrap"):
return "geocode"
raise FileNotFoundError(
"Unable to resume the ISCE2 geocode stage. Missing PICKLE state for unwrap."
)
def print_summary(
@@ -804,6 +1035,25 @@ def print_summary(
print(f"BBox: {config.bbox if config.bbox is not None else 'auto'}")
print(f"Target grid: {config.target_grid_size_m} m")
print(f"Geo posting: {config.geo_posting_deg:.12f} deg")
print(
"Enhancement: "
f"split_spectrum={config.ionosphere_correction}, "
f"ionosphere={config.ionosphere_correction}, "
f"dense_offsets={config.dense_offsets}, "
f"rubbersheet_range={config.rubbersheet_range}, "
f"rubbersheet_azimuth={config.rubbersheet_azimuth}"
)
print(
"Dense params: "
f"window={config.dense_window_width}x{config.dense_window_height}, "
f"search={config.dense_search_width}x{config.dense_search_height}, "
f"skip={config.dense_skip_width}x{config.dense_skip_height}"
)
print(
"Rubber mask: "
f"snr_threshold={config.rubber_sheet_snr_threshold}, "
f"filter_size={config.rubber_sheet_filter_size}"
)
print(
"Geocode list: "
+ (
@@ -816,6 +1066,11 @@ def print_summary(
def main() -> int:
args = parse_args()
if not args.dry_run:
dependency_error = validate_runtime_dependencies(args)
if dependency_error:
print(dependency_error, file=sys.stderr)
return 2
resume_from = str(args.resume_from or "").strip().lower()
start_stage = resume_from or PIPELINE_STAGE_ORDER[0]
task_dir = normalize_linux_path(args.task_dir).resolve()
@@ -873,12 +1128,30 @@ def main() -> int:
bbox=bbox,
target_grid_size_m=args.target_grid_size_m,
geo_posting_deg=geo_posting_deg,
geocode_products=None if args.full_geocode else list(DEFAULT_EXPORT_GEOCODE_PRODUCTS),
geocode_products=(
None
if args.full_geocode
else build_default_geocode_products(
ionosphere_correction=bool(args.ionosphere_correction)
)
),
ionosphere_correction=bool(args.ionosphere_correction),
dense_offsets=bool(args.dense_offsets),
rubbersheet_range=bool(args.rubbersheet_range),
rubbersheet_azimuth=bool(args.rubbersheet_azimuth),
rubber_sheet_snr_threshold=float(args.rubber_sheet_snr_threshold),
rubber_sheet_filter_size=int(args.rubber_sheet_filter_size),
dense_window_width=int(args.dense_window_width),
dense_window_height=int(args.dense_window_height),
dense_search_width=int(args.dense_search_width),
dense_search_height=int(args.dense_search_height),
dense_skip_width=int(args.dense_skip_width),
dense_skip_height=int(args.dense_skip_height),
)
if start_stage == PIPELINE_STAGE_ORDER[0]:
guard_large_unprepared_base_dem(config.dem_path)
if resume_from in {"unwrap", "geocode", "export"}:
if resume_from in {"unwrap", "geocode"}:
ensure_geocode_bbox(work_dir, config, args.bbox_margin)
if should_run_stage(start_stage, "geocode"):
prepare_geocode_dem(work_dir, config)
@@ -908,20 +1181,42 @@ def main() -> int:
write_stripmap_xml(xml_path, config)
if should_run_stage(start_stage, "unwrap"):
prepare_snaphu_resume(work_dir, config.bbox)
unwrap_start_step = resolve_unwrap_start_step(
work_dir,
ionosphere_correction=config.ionosphere_correction,
)
unwrap_end_step = "ionosphere" if config.ionosphere_correction else "unwrap"
unwrap_stage_name = "02_to_ionosphere" if config.ionosphere_correction else "02_to_unwrap"
run_logged(
"02_unwrap_snaphu",
[sys.executable, app_py.as_posix(), xml_path.as_posix(), "--steps", "--start=unwrap", "--end=unwrap"],
unwrap_stage_name,
[
sys.executable,
app_py.as_posix(),
xml_path.as_posix(),
"--steps",
f"--start={unwrap_start_step}",
f"--end={unwrap_end_step}",
],
cwd=work_dir,
log_path=work_dir / "02_unwrap_snaphu.log",
log_path=work_dir / f"{unwrap_stage_name}.log",
)
if should_run_stage(start_stage, "geocode"):
prepare_geocode_resume(work_dir)
geocode_start_step = resolve_geocode_start_step(
work_dir,
ionosphere_correction=config.ionosphere_correction,
)
cleanup_geocode_outputs(work_dir, config.geocode_products)
run_logged(
"03_geocode",
[sys.executable, app_py.as_posix(), xml_path.as_posix(), "--steps", "--start=geocode", "--end=geocode"],
[
sys.executable,
app_py.as_posix(),
xml_path.as_posix(),
"--steps",
f"--start={geocode_start_step}",
"--end=geocode",
],
cwd=work_dir,
log_path=work_dir / "03_geocode.log",
)
@@ -936,6 +1231,8 @@ def main() -> int:
coh_threshold=args.coh_threshold,
reference_mode=args.reference_mode,
reference_coh_threshold=args.reference_coh_threshold,
deramp_mode=args.deramp_mode,
deramp_coh_threshold=args.deramp_coh_threshold,
include_disp_full=args.include_disp_full,
)