#!/usr/bin/env python3 from __future__ import annotations import argparse import ast import importlib.util import os import re import shutil import subprocess import sys import time import xml.etree.ElementTree as ET 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, ) DEFAULT_TARGET_GRID_SIZE_M = 10 METERS_PER_DEGREE = 111320.0 LARGE_BASE_DEM_PIXEL_THRESHOLD = 200_000_000 PIPELINE_STAGE_ORDER = ("filter", "unwrap", "geocode", "export") 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 class Scene: role: str tiff_path: Path meta_path: Path date_yyyymmdd: str satellite: str orbit_xml_path: Path @dataclass class PipelineConfig: task_name: str output_prefix: str dem_path: Path reference: Scene secondary: Scene bbox: list[float] | None 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: script_dir = Path(__file__).resolve().parent repo_root = script_dir.parent parser = argparse.ArgumentParser( 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( "--task-name", default=None, help="Override the task name used for work directory and default outputs", ) parser.add_argument( "--work-root", default=str(script_dir / "jobs"), help="Root directory for ISCE2 work folders", ) parser.add_argument( "--work-dir", default=None, help="Explicit work directory. Overrides --work-root/", ) parser.add_argument( "--output-dir", default=None, help="Directory for final GeoTIFFs. Default: work_dir", ) parser.add_argument( "--output-prefix", default=None, help="Prefix for final output filenames. Default: task name", ) parser.add_argument( "--orbit-root", default=str(repo_root / "orbit"), help="Directory containing LT1A_GpsData_GAS_C_YYYYMMDD.txt", ) parser.add_argument( "--orbit-output-dir", default=None, help="Directory to place generated orbit XML files. Default: work_dir/orbits", ) parser.add_argument( "--dem", default=None, help="DEM base path. Default: auto-detect the prepared WGS84 DEM", ) parser.add_argument( "--bbox", default=None, help="Optional geocode bounding box: south,north,west,east", ) parser.add_argument( "--bbox-margin", type=float, default=0.05, help="Auto-expand topo estimated bbox by this many degrees on each side", ) parser.add_argument( "--orbit-margin-sec", type=float, default=60.0, help="Seconds to expand around scene time when clipping precise orbit, must be between 60 and 120", ) parser.add_argument( "--master-dir-name", default="master", help="Subdirectory name for the reference scene inside the task directory", ) parser.add_argument( "--slave-dir-name", default="slave", help="Subdirectory name for the secondary scene inside the task directory", ) parser.add_argument( "--scene-glob", default="*.tiff", help="Glob pattern used to find scene files inside master/slave directories", ) parser.add_argument( "--prefer-scene-keyword", default="_SLC_", help="Prefer matching files containing this keyword when multiple scene files are present", ) parser.add_argument( "--coh-threshold", type=float, default=0.05, help="Coherence threshold for *_disp.tif export", ) parser.add_argument( "--reference-mode", choices=REFERENCE_MODE_CHOICES, default=DEFAULT_REFERENCE_MODE, help="Reference normalization mode applied during final displacement export", ) parser.add_argument( "--reference-coh-threshold", type=float, 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, default=DEFAULT_TARGET_GRID_SIZE_M, help="Target grid size in meters used to control multilook scale and geocoding spacing", ) parser.add_argument( "--include-disp-full", action="store_true", help="Also export the unmasked displacement GeoTIFF for debugging", ) parser.add_argument( "--full-geocode", 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, default=None, help="Resume from an existing work directory starting at the given stage.", ) parser.add_argument( "--wavelength", type=float, default=DEFAULT_WAVELENGTH, help="Radar wavelength in meters", ) parser.add_argument( "--force", action="store_true", help="Delete an existing work directory before rerunning", ) parser.add_argument( "--reference-satellite", default=None, help="Optional LT-1 satellite for the reference/master scene (LT1A or LT1B)", ) parser.add_argument( "--secondary-satellite", default=None, help="Optional LT-1 satellite for the secondary/slave scene (LT1A or LT1B)", ) parser.add_argument( "--dry-run", action="store_true", help="Resolve inputs and print the planned configuration without running ISCE2", ) args = parser.parse_args() if args.orbit_margin_sec < 60 or args.orbit_margin_sec > 120: raise ValueError("--orbit-margin-sec must be between 60 and 120 seconds") if args.target_grid_size_m <= 0: 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 app_path = Path(isce.__file__).resolve().parent / "applications" / "stripmapApp.py" if not app_path.exists(): raise FileNotFoundError(f"stripmapApp.py not found: {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("\\\\"): raise ValueError("UNC paths are not supported directly. Mount them in WSL first.") match = re.match(r"^([A-Za-z]):[\\/](.*)$", text) if match: drive = match.group(1).lower() rest = match.group(2).replace("\\", "/") return Path(f"/mnt/{drive}/{rest}") return Path(text) def choose_scene_tiff(scene_dir: Path, scene_glob: str, prefer_scene_keyword: str) -> Path: candidates = sorted(scene_dir.glob(scene_glob)) if not candidates: raise FileNotFoundError(f"No scene file matching {scene_glob} found in {scene_dir}") slc_candidates = [path for path in candidates if prefer_scene_keyword in path.name] if len(slc_candidates) == 1: return slc_candidates[0] if len(candidates) == 1: return candidates[0] raise RuntimeError( f"Expected one scene file in {scene_dir}; found {len(candidates)} matches for {scene_glob}" ) def scene_meta_from_tiff(tiff_path: Path) -> Path: candidates = [tiff_path.with_suffix(".meta.xml")] legacy_path = Path(str(tiff_path).replace(".tiff", ".meta.xml")) if legacy_path not in candidates: candidates.append(legacy_path) for meta_path in candidates: if meta_path.exists(): return meta_path searched = ", ".join(str(path) for path in candidates) raise FileNotFoundError(f"Missing meta XML for {tiff_path}. Searched: {searched}") def extract_scene_date(name: str) -> str: match = re.search(r"_(\d{8})_SLC_", name) if not match: match = re.search(r"(\d{8})", name) if not match: raise ValueError(f"Unable to extract scene date from filename: {name}") return match.group(1) def normalize_lt1_satellite(value: str | None) -> str: text = str(value or "").strip().upper().replace("-", "").replace("_", "") if "LT1A" in text or text in {"A", "LTA"}: return "LT1A" if "LT1B" in text or text in {"B", "LTB"}: return "LT1B" return "" def extract_scene_satellite_from_name(name: str) -> str: match = re.search(r"(LT1[AB])", str(name or ""), re.IGNORECASE) return normalize_lt1_satellite(match.group(1) if match else "") def extract_scene_satellite_from_meta(meta_path: Path) -> str: try: root = ET.parse(meta_path).getroot() except Exception: return "" for element in root.iter(): tag = str(element.tag or "").rsplit("}", 1)[-1].strip().lower() if tag not in {"mission", "satellite", "platform", "platformid", "missionid"}: continue satellite = normalize_lt1_satellite(element.text) if satellite: return satellite return "" def resolve_scene_satellite( tiff_path: Path, meta_path: Path, explicit_satellite: str | None = None, ) -> str: explicit = normalize_lt1_satellite(explicit_satellite) name_satellite = extract_scene_satellite_from_name(tiff_path.name) meta_satellite = extract_scene_satellite_from_meta(meta_path) if explicit: if name_satellite and name_satellite != explicit: raise ValueError( f"Explicit satellite {explicit} does not match filename for {tiff_path.name}: {name_satellite}" ) if meta_satellite and meta_satellite != explicit: raise ValueError( f"Explicit satellite {explicit} does not match metadata for {meta_path.name}: {meta_satellite}" ) return explicit if name_satellite and meta_satellite and name_satellite != meta_satellite: raise ValueError( f"Satellite mismatch between filename and metadata for {tiff_path.name}: " f"{name_satellite} vs {meta_satellite}" ) if name_satellite: return name_satellite if meta_satellite: return meta_satellite raise ValueError(f"Unable to resolve LT-1 satellite from {tiff_path} / {meta_path}") def ensure_orbit_xml( date_yyyymmdd: str, satellite: str, annotation_xml: Path, orbit_root: Path, orbit_out_dir: Path, margin_sec: float, ) -> Path: resolution = ensure_lt1_orbit_xml( date_yyyymmdd=date_yyyymmdd, satellite=satellite, annotation_xml=annotation_xml, orbit_root=orbit_root, orbit_output_dir=orbit_out_dir, margin_sec=margin_sec, ) return resolution.path def resolve_dem(dem_value: str | None) -> Path: dem_path = resolve_prepared_dem_path( explicit_path=dem_value, env_values=None, default_candidates=DEFAULT_WSL_DEM_CANDIDATES, 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) raise FileNotFoundError( "Unable to resolve a prepared DEM with ISCE wrappers. " f"Searched: {searched}" ) def _read_xml_property_value(root: ET.Element, name: str) -> str: for prop in root.findall("property"): if str(prop.get("name") or "").strip() != name: continue return str(prop.findtext("value") or "").strip() return "" def read_dem_dimensions(dem_path: Path) -> tuple[int, int] | None: xml_path = Path(str(dem_path) + ".xml") if not xml_path.exists(): return None root = ET.fromstring(xml_path.read_text(encoding="utf-8", errors="ignore")) width_text = _read_xml_property_value(root, "width") length_text = _read_xml_property_value(root, "length") if not width_text or not length_text: return None try: return int(float(width_text)), int(float(length_text)) except ValueError: return None def has_prepared_dem_sibling(dem_path: Path) -> bool: if dem_path.as_posix().lower().endswith(".wgs84"): return True sibling = Path(str(dem_path) + ".wgs84") return sibling.exists() and Path(str(sibling) + ".xml").exists() def guard_large_unprepared_base_dem(dem_path: Path) -> None: if has_prepared_dem_sibling(dem_path): return dimensions = read_dem_dimensions(dem_path) if dimensions is None: return width, length = dimensions pixel_count = width * length if pixel_count < LARGE_BASE_DEM_PIXEL_THRESHOLD: return raise RuntimeError( "Configured DEM resolves to a large base raster without a prepared '.wgs84' sibling. " f"Selected DEM: {dem_path} ({width}x{length}, {pixel_count} pixels). " "A fresh ISCE2 run would spend a very long time rebuilding the geoid-corrected DEM during " "verifyDEM/topo. Prepare '.wgs84' once, or point ISCE2_DEM_PATH directly to the " "prepared file before starting a fresh run." ) def resolve_task( task_dir: Path, orbit_root: Path, orbit_out_dir: Path, margin_sec: float, master_dir_name: str, slave_dir_name: str, scene_glob: str, prefer_scene_keyword: str, reference_satellite: str | None = None, secondary_satellite: str | None = None, ) -> tuple[Scene, Scene]: scenes: list[Scene] = [] satellite_hints = { "master": reference_satellite, "slave": secondary_satellite, } for role, subdir in (("master", master_dir_name), ("slave", slave_dir_name)): scene_dir = task_dir / subdir if not scene_dir.exists(): raise FileNotFoundError(f"Missing task subdirectory: {scene_dir}") tiff_path = choose_scene_tiff(scene_dir, scene_glob, prefer_scene_keyword) meta_path = scene_meta_from_tiff(tiff_path) date_yyyymmdd = extract_scene_date(tiff_path.name) satellite = resolve_scene_satellite( tiff_path=tiff_path, meta_path=meta_path, explicit_satellite=satellite_hints.get(role), ) orbit_xml_path = ensure_orbit_xml( date_yyyymmdd=date_yyyymmdd, satellite=satellite, annotation_xml=meta_path, orbit_root=orbit_root, orbit_out_dir=orbit_out_dir, margin_sec=margin_sec, ) scenes.append( Scene( role=role, tiff_path=tiff_path, meta_path=meta_path, date_yyyymmdd=date_yyyymmdd, satellite=satellite, orbit_xml_path=orbit_xml_path, ) ) return scenes[0], scenes[1] def render_bbox(bbox: list[float] | None) -> str: if bbox is None: return "" values = ", ".join(f"{value:.10f}".rstrip("0").rstrip(".") for value in bbox) return f' [{values}]\n' def render_string_list(name: str, values: list[str] | None) -> str: if not values: return "" rendered = ", ".join(repr(str(value)) for value in values if str(value).strip()) return f' [{rendered}]\n' if rendered else "" 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" {str(config.dense_offsets)}\n" f" {str(config.rubbersheet_range)}\n" f" {str(config.rubbersheet_azimuth)}\n" f" {config.rubber_sheet_snr_threshold}\n" f" {config.rubber_sheet_filter_size}\n" f" {config.dense_window_width}\n" f" {config.dense_window_height}\n" f" {config.dense_search_width}\n" f" {config.dense_search_height}\n" f" {config.dense_skip_width}\n" f" {config.dense_skip_height}\n" ) text = ( "\n" " \n" " LUTAN1\n" " LUTAN1\n" " LUTAN1\n" " xml\n" " True\n" " snaphu\n" f" {str(config.ionosphere_correction)}\n" f" {str(config.ionosphere_correction)}\n" f" {config.target_grid_size_m}\n" f" {config.geo_posting_deg:.12f}\n" f"{bbox_xml}" f"{geocode_list_xml}" f"{enhancement_props}" f" {config.dem_path.as_posix()}\n" "\n" " \n" f" {config.reference.tiff_path.as_posix()}\n" f" {config.reference.orbit_xml_path.as_posix()}\n" " reference\n" " \n" "\n" " \n" f" {config.secondary.tiff_path.as_posix()}\n" f" {config.secondary.orbit_xml_path.as_posix()}\n" " secondary\n" " \n" " \n" "\n" ) xml_path.write_text(text, encoding="utf-8") def run_logged(stage_name: str, cmd: list[str], cwd: Path, log_path: Path) -> None: started_monotonic = time.monotonic() started_text = time.strftime("%Y-%m-%d %H:%M:%S") log_path.parent.mkdir(parents=True, exist_ok=True) print(f"[{stage_name}] Starting at {started_text}", flush=True) print("Running:", flush=True) print(" " + " ".join(cmd), flush=True) print(f"Log: {log_path}", flush=True) with log_path.open("w", encoding="utf-8") as handle: handle.write(f"[{stage_name}] Starting at {started_text}\n") handle.write("Running:\n") handle.write(" " + " ".join(cmd) + "\n") handle.write(f"Log: {log_path}\n") handle.flush() child_env = build_process_env() child_env["PYTHONUNBUFFERED"] = "1" proc = subprocess.Popen( cmd, cwd=cwd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True, bufsize=1, env=child_env, ) assert proc.stdout is not None for line in proc.stdout: sys.stdout.write(line) sys.stdout.flush() handle.write(line) handle.flush() status = proc.wait() elapsed_seconds = time.monotonic() - started_monotonic handle.write(f"[{stage_name}] Finished with exit code {status} after {elapsed_seconds:.1f}s\n") handle.flush() print( f"[{stage_name}] Finished with exit code {status} after {elapsed_seconds:.1f}s", flush=True, ) if status != 0: raise RuntimeError(f"Command failed with exit code {status}: {' '.join(cmd)}") def parse_bbox_arg(value: str | None) -> list[float] | None: if value is None: return None parts = [item.strip() for item in value.split(",")] if len(parts) != 4: raise ValueError("Bounding box must be south,north,west,east") return [float(item) for item in parts] def load_estimated_bbox(topo_xml: Path) -> list[float]: root = ET.fromstring(topo_xml.read_text(encoding="utf-8")) for prop in root.findall("property"): if prop.attrib.get("name") == "estimatedboundingbox": value = prop.findtext("value") if not value: break bbox = ast.literal_eval(value) return [float(item) for item in bbox] raise ValueError(f"estimatedboundingbox not found in {topo_xml}") def expand_bbox(bbox: list[float], margin: float) -> list[float]: south, north, west, east = bbox return [ max(-90.0, south - margin), min(90.0, north + margin), max(-180.0, west - margin), min(180.0, east + margin), ] def resolve_auto_geocode_bbox(work_dir: Path, bbox_margin: float) -> tuple[list[float], list[float]]: topo_xml = work_dir / "PICKLE" / "topo.xml" if not topo_xml.exists(): raise FileNotFoundError("topo step output is missing; cannot resolve the geocode bounding box.") estimated_bbox = load_estimated_bbox(topo_xml) return estimated_bbox, expand_bbox(estimated_bbox, bbox_margin) def ensure_geocode_bbox(work_dir: Path, config: PipelineConfig, bbox_margin: float) -> None: if config.bbox is not None: return estimated_bbox, expanded_bbox = resolve_auto_geocode_bbox(work_dir, bbox_margin) config.bbox = expanded_bbox print(f"Auto bbox from topo: {estimated_bbox}") print(f"Expanded bbox used for geocode: {config.bbox}") def cleanup_geocode_outputs(work_dir: Path, geocode_products: list[str] | None) -> None: if not geocode_products: return removed: list[Path] = [] for product in geocode_products: base_path = work_dir / product for suffix in (".geo", ".geo.xml", ".geo.vrt", ".geo.aux.xml"): candidate = Path(str(base_path) + suffix) if candidate.exists(): candidate.unlink() removed.append(candidate) if removed: print(f"Removed {len(removed)} stale geocode output file(s).") def cleanup_dem_subset_outputs(base_path: Path) -> None: for suffix in ("", ".hdr", ".xml", ".vrt", ".aux.xml"): candidate = Path(str(base_path) + suffix) if candidate.exists(): candidate.unlink() def prepare_geocode_dem_subset(work_dir: Path, source_dem_path: Path, bbox: list[float]) -> Path: import isce # noqa: F401 # Ensures the bundled ISCE packages are initialized on sys.path. import isceobj from osgeo import gdal gdal.UseExceptions() source_xml = Path(str(source_dem_path) + ".xml") if not source_xml.exists(): raise FileNotFoundError(f"Missing DEM XML sidecar: {source_xml}") source_vrt = Path(str(source_dem_path) + ".vrt") source_open_path = source_vrt if source_vrt.exists() else source_dem_path if not source_open_path.exists(): raise FileNotFoundError(f"Missing DEM source for geocode subset: {source_open_path}") subset_base = work_dir / "geocode_dem" cleanup_dem_subset_outputs(subset_base) south, north, west, east = bbox ds = gdal.Translate( subset_base.as_posix(), source_open_path.as_posix(), format="ENVI", projWin=[west, north, east, south], ) if ds is None: raise RuntimeError(f"Failed to crop DEM subset from {source_open_path}") width = int(ds.RasterXSize or 0) length = int(ds.RasterYSize or 0) geotransform = ds.GetGeoTransform(can_return_null=True) ds = None if width <= 0 or length <= 0 or geotransform is None: raise RuntimeError("Cropped DEM subset is empty or missing georeferencing metadata.") source_dem = isceobj.createDemImage() source_dem.load(source_xml.as_posix()) dem_reference = str(source_dem.reference or "").strip() or "UNKNOWN" source_dem.filename = subset_base.as_posix() source_dem.width = width source_dem.length = length source_dem.coord1.coordStart = geotransform[0] source_dem.coord1.coordDelta = geotransform[1] source_dem.coord1.coordSize = width source_dem.coord2.coordStart = geotransform[3] source_dem.coord2.coordDelta = geotransform[5] source_dem.coord2.coordSize = length source_dem.dump(subset_base.as_posix() + ".xml") source_dem.renderVRT() print( "Prepared geocode DEM subset: " f"{subset_base} ({width}x{length}, reference={dem_reference})" ) return subset_base def prepare_geocode_dem(work_dir: Path, config: PipelineConfig) -> None: if config.bbox is None: raise ValueError("Cannot prepare a geocode DEM subset without a resolved bbox.") config.dem_path = prepare_geocode_dem_subset(work_dir, config.dem_path, config.bbox) def should_run_stage(start_stage: str, stage_name: str) -> bool: start_index = PIPELINE_STAGE_ORDER.index(start_stage) stage_index = PIPELINE_STAGE_ORDER.index(stage_name) return stage_index >= start_index def has_pickle_state(work_dir: Path, state_name: str) -> bool: pickle_dir = work_dir / "PICKLE" return (pickle_dir / state_name).exists() and (pickle_dir / f"{state_name}.xml").exists() def copy_pickle_state(work_dir: Path, source_state: str, target_state: str) -> Path: pickle_dir = work_dir / "PICKLE" src = pickle_dir / source_state src_xml = pickle_dir / f"{source_state}.xml" dst = pickle_dir / target_state dst_xml = pickle_dir / f"{target_state}.xml" if not src.exists() or not src_xml.exists(): raise FileNotFoundError( f"PICKLE/{source_state} state is missing; cannot prepare PICKLE/{target_state}." ) shutil.copy2(src, dst) shutil.copy2(src_xml, dst_xml) return dst_xml def prepare_no_ionosphere_snaphu_resume(work_dir: Path, bbox: list[float] | None) -> None: dst_xml = copy_pickle_state(work_dir, "filter", "filter_high_band") 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") print("Prepared no-ionosphere SNAPHU resume state: PICKLE/filter_high_band", flush=True) def prepare_no_ionosphere_geocode_resume(work_dir: Path) -> None: copy_pickle_state(work_dir, "unwrap", "ionosphere") print("Prepared no-ionosphere geocode resume state: PICKLE/ionosphere", flush=True) 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 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." ) 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( task_dir: Path, work_dir: Path, output_dir: Path, config: PipelineConfig, ) -> None: print(f"Task dir: {task_dir}") print(f"Task name: {config.task_name}") print(f"Output prefix:{config.output_prefix}") print(f"Work dir: {work_dir}") print(f"Output dir: {output_dir}") print(f"DEM: {config.dem_path}") print(f"Reference: {config.reference.tiff_path} [{config.reference.satellite}]") print(f"Secondary: {config.secondary.tiff_path} [{config.secondary.satellite}]") print(f"Ref orbit: {config.reference.orbit_xml_path}") print(f"Sec orbit: {config.secondary.orbit_xml_path}") 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: " + ( ", ".join(config.geocode_products) if config.geocode_products else "ISCE2 default" ) ) 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() if not task_dir.exists(): raise FileNotFoundError(f"Task directory not found: {task_dir}") task_name = args.task_name or task_dir.name output_prefix = args.output_prefix or task_name work_root = normalize_linux_path(args.work_root).resolve() work_dir = normalize_linux_path(args.work_dir).resolve() if args.work_dir else work_root / task_name output_dir = normalize_linux_path(args.output_dir).resolve() if args.output_dir else work_dir orbit_root = normalize_linux_path(args.orbit_root).resolve() orbit_out_dir = ( normalize_linux_path(args.orbit_output_dir).resolve() if args.orbit_output_dir else work_dir / "orbits" ) if work_dir.exists(): if resume_from: pass elif args.force: shutil.rmtree(work_dir) else: raise FileExistsError(f"Work directory already exists: {work_dir}. Use --force to recreate it.") elif resume_from: raise FileNotFoundError( f"Resume requested from {resume_from}, but work directory does not exist: {work_dir}" ) work_dir.mkdir(parents=True, exist_ok=True) reference, secondary = resolve_task( task_dir=task_dir, orbit_root=orbit_root, orbit_out_dir=orbit_out_dir, margin_sec=args.orbit_margin_sec, master_dir_name=args.master_dir_name, slave_dir_name=args.slave_dir_name, scene_glob=args.scene_glob, prefer_scene_keyword=args.prefer_scene_keyword, reference_satellite=args.reference_satellite, secondary_satellite=args.secondary_satellite, ) dem_path = resolve_dem(args.dem) bbox = parse_bbox_arg(args.bbox) geo_posting_deg = meters_to_geoposting_degrees(args.target_grid_size_m) config = PipelineConfig( task_name=task_name, output_prefix=output_prefix, dem_path=dem_path, reference=reference, secondary=secondary, 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 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"}: ensure_geocode_bbox(work_dir, config, args.bbox_margin) if should_run_stage(start_stage, "geocode"): prepare_geocode_dem(work_dir, config) print_summary(task_dir=task_dir, work_dir=work_dir, output_dir=output_dir, config=config) xml_path = work_dir / f"{task_name}_stripmap.xml" write_stripmap_xml(xml_path, config) if args.dry_run: print(f"Generated XML: {xml_path}") return 0 app_py = locate_stripmap_app() if should_run_stage(start_stage, "filter"): run_logged( "01_to_filter", [sys.executable, app_py.as_posix(), xml_path.as_posix(), "--steps", "--end=filter"], cwd=work_dir, log_path=work_dir / "01_to_filter.log", ) if should_run_stage(start_stage, "geocode") and config.bbox is None: ensure_geocode_bbox(work_dir, config, args.bbox_margin) prepare_geocode_dem(work_dir, config) write_stripmap_xml(xml_path, config) if should_run_stage(start_stage, "unwrap"): if not config.ionosphere_correction: prepare_no_ionosphere_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( 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 / f"{unwrap_stage_name}.log", ) if should_run_stage(start_stage, "geocode"): if not config.ionosphere_correction: prepare_no_ionosphere_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", f"--start={geocode_start_step}", "--end=geocode", ], cwd=work_dir, log_path=work_dir / "03_geocode.log", ) outputs: dict[str, Path] = {} if should_run_stage(start_stage, "export"): outputs = export_products( work_dir=work_dir, output_dir=output_dir, prefix=output_prefix, wavelength=args.wavelength, 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, ) print("Pipeline finished.") for key, path in outputs.items(): print(f"{key}: {path}") return 0 if __name__ == "__main__": raise SystemExit(main())