From 46474274c2b9d1f5afdc0719cead60700be40b63 Mon Sep 17 00:00:00 2001 From: Harmon Date: Sun, 3 May 2026 15:41:08 +0800 Subject: [PATCH] Refactor GF3 L1A to L2 pipeline --- .gitignore | 9 ++ Decompression.py | 30 ++++++ README.md | 59 +++++++++++ environment.yml | 10 ++ gf3_L1A_To_L2.py | 16 +++ pyproject.toml | 24 +++++ src/gf3_l1a2l2/__init__.py | 16 +++ src/gf3_l1a2l2/__main__.py | 5 + src/gf3_l1a2l2/calibration.py | 88 +++++++++++++++ src/gf3_l1a2l2/cli.py | 107 +++++++++++++++++++ src/gf3_l1a2l2/config.py | 41 +++++++ src/gf3_l1a2l2/constants.py | 6 ++ src/gf3_l1a2l2/decompress.py | 112 ++++++++++++++++++++ src/gf3_l1a2l2/discovery.py | 90 ++++++++++++++++ src/gf3_l1a2l2/errors.py | 18 ++++ src/gf3_l1a2l2/geocorrection.py | 78 ++++++++++++++ src/gf3_l1a2l2/log.py | 34 ++++++ src/gf3_l1a2l2/metadata.py | 145 +++++++++++++++++++++++++ src/gf3_l1a2l2/models.py | 66 ++++++++++++ src/gf3_l1a2l2/paths.py | 43 ++++++++ src/gf3_l1a2l2/pipeline.py | 182 ++++++++++++++++++++++++++++++++ src/gf3_l1a2l2/rpc.py | 70 ++++++++++++ tests/fixtures/sample.meta.xml | 8 ++ tests/fixtures/sample.rpb | 16 +++ tests/test_metadata.py | 19 ++++ tests/test_rpc.py | 18 ++++ 26 files changed, 1310 insertions(+) create mode 100644 .gitignore create mode 100644 Decompression.py create mode 100644 README.md create mode 100644 environment.yml create mode 100644 gf3_L1A_To_L2.py create mode 100644 pyproject.toml create mode 100644 src/gf3_l1a2l2/__init__.py create mode 100644 src/gf3_l1a2l2/__main__.py create mode 100644 src/gf3_l1a2l2/calibration.py create mode 100644 src/gf3_l1a2l2/cli.py create mode 100644 src/gf3_l1a2l2/config.py create mode 100644 src/gf3_l1a2l2/constants.py create mode 100644 src/gf3_l1a2l2/decompress.py create mode 100644 src/gf3_l1a2l2/discovery.py create mode 100644 src/gf3_l1a2l2/errors.py create mode 100644 src/gf3_l1a2l2/geocorrection.py create mode 100644 src/gf3_l1a2l2/log.py create mode 100644 src/gf3_l1a2l2/metadata.py create mode 100644 src/gf3_l1a2l2/models.py create mode 100644 src/gf3_l1a2l2/paths.py create mode 100644 src/gf3_l1a2l2/pipeline.py create mode 100644 src/gf3_l1a2l2/rpc.py create mode 100644 tests/fixtures/sample.meta.xml create mode 100644 tests/fixtures/sample.rpb create mode 100644 tests/test_metadata.py create mode 100644 tests/test_rpc.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..37d61be --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +.venv/ +__pycache__/ +*.py[cod] +.pytest_cache/ +.mypy_cache/ +.ruff_cache/ +/.tmp_tests/ +logs/ +GMTED2010.jp2 diff --git a/Decompression.py b/Decompression.py new file mode 100644 index 0000000..6863954 --- /dev/null +++ b/Decompression.py @@ -0,0 +1,30 @@ +"""Compatibility wrapper for archive extraction.""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +_SRC = Path(__file__).resolve().parent / "src" +if _SRC.is_dir(): + sys.path.insert(0, str(_SRC)) + +from gf3_l1a2l2.decompress import extract_compressed_files + + +def main() -> int: + parser = argparse.ArgumentParser(description="Extract GF3 archives") + parser.add_argument("input", type=Path, help="folder containing .tar.gz or .zip archives") + parser.add_argument("--overwrite", action="store_true", help="overwrite existing extracted folders") + args = parser.parse_args() + + results = extract_compressed_files(args.input, overwrite=args.overwrite) + failed = [item for item in results if item.status == "failed"] + for item in results: + print(f"{item.status}\t{item.archive}\t{item.output_dir}\t{item.message}") + return 1 if failed else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/README.md b/README.md new file mode 100644 index 0000000..800ed51 --- /dev/null +++ b/README.md @@ -0,0 +1,59 @@ +# GF3 L1A to L2 Pipeline + +This repository contains a Python package and CLI for converting GF3 L1A +products to L2 GeoTIFF outputs. + +The heavy work is still handled by GDAL, Rasterio, and NumPy. The package adds +pipeline structure around the original script: explicit configuration, safe path +handling, logging, restart-friendly skipping, per-scene error isolation, and +block-based radiometric calibration. + +## Install + +Use an environment that already has a working GDAL/OSGeo Python binding. On +Windows, Conda/Miniforge is usually the least fragile route. + +```powershell +conda env create -f environment.yml +conda activate gf3-l1a2l2 +pip install -e . --no-deps +``` + +## Run + +```powershell +gf3-l1a2l2 run --input D:\GF3\L1A_BATCH --output E:\GF3\L2_Image --dem .\GMTED2010.jp2 +``` + +Common options: + +```powershell +gf3-l1a2l2 run ` + --input D:\GF3\L1A_BATCH ` + --output E:\GF3\L2_Image ` + --dem .\GMTED2010.jp2 ` + --workers 1 ` + --x-res 0.0002 ` + --y-res 0.0002 ` + --dst-srs EPSG:4326 +``` + +By default, output files are written under one subdirectory per discovered +scene. Use `--flat-output` to write all outputs directly into the output +directory. + +The legacy files are now compatibility wrappers: + +```powershell +python gf3_L1A_To_L2.py run --input D:\GF3\L1A_BATCH --output E:\GF3\L2_Image --dem .\GMTED2010.jp2 +python Decompression.py D:\GF3\L1A_BATCH +``` + +## Notes + +- The CLI extracts `.tar.gz` and `.zip` products before processing unless + `--skip-decompress` is used. +- Existing L2 outputs are skipped unless `--overwrite` is used. +- Logs are written to `\logs`. +- Without representative GF3 products, verification is limited to parser unit + tests and Python syntax/import checks. diff --git a/environment.yml b/environment.yml new file mode 100644 index 0000000..10b28ef --- /dev/null +++ b/environment.yml @@ -0,0 +1,10 @@ +name: gf3-l1a2l2 +channels: + - conda-forge +dependencies: + - python=3.11 + - gdal + - rasterio + - numpy + - tqdm + - pip diff --git a/gf3_L1A_To_L2.py b/gf3_L1A_To_L2.py new file mode 100644 index 0000000..24b6c2d --- /dev/null +++ b/gf3_L1A_To_L2.py @@ -0,0 +1,16 @@ +"""Compatibility wrapper for the packaged GF3 L1A to L2 pipeline.""" + +from __future__ import annotations + +import sys +from pathlib import Path + +_SRC = Path(__file__).resolve().parent / "src" +if _SRC.is_dir(): + sys.path.insert(0, str(_SRC)) + +from gf3_l1a2l2.cli import main + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..bad80b9 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,24 @@ +[build-system] +requires = ["setuptools>=68", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "gf3-l1a2l2-pipeline" +version = "0.1.0" +description = "GF3 L1A to L2 processing pipeline" +readme = "README.md" +requires-python = ">=3.10" +dependencies = [ + "numpy>=1.22", + "rasterio>=1.3", + "tqdm>=4.64", +] + +[project.scripts] +gf3-l1a2l2 = "gf3_l1a2l2.cli:main" + +[tool.setuptools] +package-dir = {"" = "src"} + +[tool.setuptools.packages.find] +where = ["src"] diff --git a/src/gf3_l1a2l2/__init__.py b/src/gf3_l1a2l2/__init__.py new file mode 100644 index 0000000..283dafb --- /dev/null +++ b/src/gf3_l1a2l2/__init__.py @@ -0,0 +1,16 @@ +"""GF3 L1A to L2 processing pipeline.""" + +__all__ = ["PipelineConfig", "run_pipeline"] +__version__ = "0.1.0" + + +def __getattr__(name: str): + if name == "PipelineConfig": + from gf3_l1a2l2.config import PipelineConfig + + return PipelineConfig + if name == "run_pipeline": + from gf3_l1a2l2.pipeline import run_pipeline + + return run_pipeline + raise AttributeError(name) diff --git a/src/gf3_l1a2l2/__main__.py b/src/gf3_l1a2l2/__main__.py new file mode 100644 index 0000000..1ae10bd --- /dev/null +++ b/src/gf3_l1a2l2/__main__.py @@ -0,0 +1,5 @@ +from gf3_l1a2l2.cli import main + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/gf3_l1a2l2/calibration.py b/src/gf3_l1a2l2/calibration.py new file mode 100644 index 0000000..0483931 --- /dev/null +++ b/src/gf3_l1a2l2/calibration.py @@ -0,0 +1,88 @@ +from __future__ import annotations + +import logging +import math +from pathlib import Path +from typing import Iterator, Optional + +import numpy as np +import rasterio +from rasterio.windows import Window + +from gf3_l1a2l2.errors import ProcessingError +from gf3_l1a2l2.models import CalibrationValue +from gf3_l1a2l2.paths import output_name + + +def radiometric_calibrate_l1a_to_l1b( + image_path: Path | str, + output_dir: Path | str, + calibration: CalibrationValue, + *, + block_size: int, + overwrite: bool = False, + logger: Optional[logging.Logger] = None, +) -> Path: + source = Path(image_path) + out_dir = Path(output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + l1b_path = out_dir / output_name(source, "L1A", "L1B") + + if l1b_path.exists() and not overwrite: + if logger: + logger.info("Skip existing L1B: %s", l1b_path) + return l1b_path + + if not math.isfinite(calibration.qualify_value): + raise ProcessingError(f"Invalid QualifyValue for {source.name}: {calibration.qualify_value}") + if not math.isfinite(calibration.calibration_const): + raise ProcessingError( + f"Invalid CalibrationConst for {source.name}: {calibration.calibration_const}" + ) + + if logger: + logger.info("Calibrating %s -> %s", source, l1b_path) + + try: + with rasterio.open(source) as src: + if src.count < 2: + raise ProcessingError(f"Expected at least 2 bands in L1A image: {source}") + + profile = src.profile.copy() + profile.update( + driver="GTiff", + count=1, + dtype="float32", + compress="lzw", + BIGTIFF="YES", + tiled=True, + ) + + with rasterio.open(l1b_path, "w", **profile) as dst: + for window in _iter_windows(src.width, src.height, block_size): + real = src.read(1, window=window, out_dtype="float32") + imag = src.read(2, window=window, out_dtype="float32") + db = _calibrate_window(real, imag, calibration) + dst.write(db, 1, window=window) + except ProcessingError: + raise + except Exception as exc: + raise ProcessingError(f"Failed to calibrate {source}") from exc + + return l1b_path + + +def _calibrate_window(real: np.ndarray, imag: np.ndarray, calibration: CalibrationValue) -> np.ndarray: + amplitude = np.hypot(real, imag) + scaled = (amplitude / np.float32(32767.0)) * np.float32(calibration.qualify_value) + with np.errstate(divide="ignore", invalid="ignore"): + db = np.float32(20.0) * np.log10(scaled) - np.float32(calibration.calibration_const) + return db.astype("float32", copy=False) + + +def _iter_windows(width: int, height: int, block_size: int) -> Iterator[Window]: + for row_off in range(0, height, block_size): + window_height = min(block_size, height - row_off) + for col_off in range(0, width, block_size): + window_width = min(block_size, width - col_off) + yield Window(col_off, row_off, window_width, window_height) diff --git a/src/gf3_l1a2l2/cli.py b/src/gf3_l1a2l2/cli.py new file mode 100644 index 0000000..21f28b9 --- /dev/null +++ b/src/gf3_l1a2l2/cli.py @@ -0,0 +1,107 @@ +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +from gf3_l1a2l2.config import PipelineConfig +from gf3_l1a2l2.decompress import extract_compressed_files +from gf3_l1a2l2.discovery import discover_scenes +from gf3_l1a2l2.errors import GF3PipelineError + + +def main(argv: list[str] | None = None) -> int: + parser = build_parser() + args = parser.parse_args(argv) + + if not hasattr(args, "handler"): + parser.print_help() + return 2 + + try: + return args.handler(args) + except GF3PipelineError as exc: + print(f"ERROR: {exc}", file=sys.stderr) + return 1 + + +def build_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser(prog="gf3-l1a2l2", description="GF3 L1A to L2 pipeline") + subparsers = parser.add_subparsers(dest="command") + + run_parser = subparsers.add_parser("run", help="extract and process GF3 L1A products") + run_parser.add_argument("--input", required=True, type=Path, help="input folder containing archives or scenes") + run_parser.add_argument("--output", required=True, type=Path, help="output folder") + run_parser.add_argument("--dem", type=Path, default=None, help="DEM file for RPC orthorectification") + run_parser.add_argument("--x-res", type=float, default=0.0002, help="output x resolution") + run_parser.add_argument("--y-res", type=float, default=0.0002, help="output y resolution") + run_parser.add_argument("--dst-srs", default="EPSG:4326", help="destination spatial reference") + run_parser.add_argument("--block-size", type=int, default=2048, help="calibration block size in pixels") + run_parser.add_argument("--workers", type=int, default=1, help="number of scenes to process concurrently") + run_parser.add_argument("--skip-decompress", action="store_true", help="do not extract archives before processing") + run_parser.add_argument("--overwrite", action="store_true", help="overwrite existing extraction and output files") + run_parser.add_argument("--flat-output", action="store_true", help="write all outputs directly to output folder") + run_parser.add_argument("--fail-fast", action="store_true", help="stop on the first failed scene or polarization") + run_parser.add_argument("--log-level", default="INFO", choices=["DEBUG", "INFO", "WARNING", "ERROR"]) + run_parser.set_defaults(handler=_run) + + decompress_parser = subparsers.add_parser("decompress", help="extract .tar.gz and .zip archives") + decompress_parser.add_argument("input", type=Path, help="folder containing archives") + decompress_parser.add_argument("--overwrite", action="store_true", help="overwrite existing extracted folders") + decompress_parser.set_defaults(handler=_decompress) + + discover_parser = subparsers.add_parser("discover", help="list discovered GF3 scenes") + discover_parser.add_argument("input", type=Path, help="input folder") + discover_parser.set_defaults(handler=_discover) + + return parser + + +def _run(args: argparse.Namespace) -> int: + from gf3_l1a2l2.pipeline import run_pipeline + + dem_path = args.dem or _default_dem_path() + config = PipelineConfig( + input_dir=args.input, + output_dir=args.output, + dem_path=dem_path, + x_res=args.x_res, + y_res=args.y_res, + dst_srs=args.dst_srs, + block_size=args.block_size, + workers=args.workers, + skip_decompress=args.skip_decompress, + overwrite=args.overwrite, + flat_output=args.flat_output, + fail_fast=args.fail_fast, + log_level=args.log_level, + ) + summary = run_pipeline(config) + return 1 if summary.failed else 0 + + +def _decompress(args: argparse.Namespace) -> int: + results = extract_compressed_files(args.input, overwrite=args.overwrite) + failed = [item for item in results if item.status == "failed"] + for item in results: + print(f"{item.status}\t{item.archive}\t{item.output_dir}\t{item.message}") + return 1 if failed else 0 + + +def _discover(args: argparse.Namespace) -> int: + scenes = discover_scenes(args.input) + for scene in scenes: + polarizations = ",".join(scene.available_polarizations) + print(f"{scene.root}\tmetadata={scene.metadata_xml.name}\tpol={polarizations}") + return 0 + + +def _default_dem_path() -> Path | None: + candidates = [ + Path.cwd() / "GMTED2010.jp2", + Path(__file__).resolve().parents[2] / "GMTED2010.jp2", + ] + for candidate in candidates: + if candidate.is_file(): + return candidate + return None diff --git a/src/gf3_l1a2l2/config.py b/src/gf3_l1a2l2/config.py new file mode 100644 index 0000000..87e190e --- /dev/null +++ b/src/gf3_l1a2l2/config.py @@ -0,0 +1,41 @@ +from __future__ import annotations + +from dataclasses import dataclass, replace +from pathlib import Path +from typing import Optional + +from gf3_l1a2l2.constants import ( + DEFAULT_BLOCK_SIZE, + DEFAULT_DST_SRS, + DEFAULT_X_RES, + DEFAULT_Y_RES, +) + + +@dataclass(frozen=True) +class PipelineConfig: + input_dir: Path + output_dir: Path + dem_path: Optional[Path] = None + x_res: float = DEFAULT_X_RES + y_res: float = DEFAULT_Y_RES + dst_srs: str = DEFAULT_DST_SRS + block_size: int = DEFAULT_BLOCK_SIZE + workers: int = 1 + skip_decompress: bool = False + overwrite: bool = False + flat_output: bool = False + keep_l1b: bool = True + fail_fast: bool = False + log_level: str = "INFO" + + def resolved(self) -> "PipelineConfig": + dem_path = self.dem_path.resolve() if self.dem_path else None + return replace( + self, + input_dir=self.input_dir.resolve(), + output_dir=self.output_dir.resolve(), + dem_path=dem_path, + workers=max(1, int(self.workers)), + block_size=max(128, int(self.block_size)), + ) diff --git a/src/gf3_l1a2l2/constants.py b/src/gf3_l1a2l2/constants.py new file mode 100644 index 0000000..55f8adf --- /dev/null +++ b/src/gf3_l1a2l2/constants.py @@ -0,0 +1,6 @@ +POLARIZATIONS = ("HH", "HV", "VH", "VV") + +DEFAULT_X_RES = 0.0002 +DEFAULT_Y_RES = 0.0002 +DEFAULT_DST_SRS = "EPSG:4326" +DEFAULT_BLOCK_SIZE = 2048 diff --git a/src/gf3_l1a2l2/decompress.py b/src/gf3_l1a2l2/decompress.py new file mode 100644 index 0000000..c9a43bd --- /dev/null +++ b/src/gf3_l1a2l2/decompress.py @@ -0,0 +1,112 @@ +from __future__ import annotations + +import logging +import shutil +import tarfile +import zipfile +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable, Optional + +try: + from tqdm import tqdm +except Exception: # pragma: no cover - only used when tqdm is unavailable + def tqdm(iterable: Iterable, **_: object) -> Iterable: + return iterable + +from gf3_l1a2l2.errors import ProcessingError + + +@dataclass(frozen=True) +class ExtractionResult: + archive: Path + output_dir: Path + status: str + message: str = "" + + +def extract_compressed_files( + folder_path: Path | str, + *, + overwrite: bool = False, + logger: Optional[logging.Logger] = None, +) -> list[ExtractionResult]: + folder = Path(folder_path).resolve() + if not folder.is_dir(): + raise ProcessingError(f"Input folder does not exist: {folder}") + + archives = sorted( + path + for path in folder.iterdir() + if path.is_file() + and (path.name.lower().endswith(".tar.gz") or path.suffix.lower() == ".zip") + ) + + results: list[ExtractionResult] = [] + if not archives: + if logger: + logger.info("No .tar.gz or .zip archives found in %s", folder) + return results + + for archive in tqdm(archives, desc="Extract archives", unit="file"): + output_dir = _archive_output_dir(archive) + try: + if output_dir.exists(): + if not overwrite: + results.append(ExtractionResult(archive, output_dir, "skipped", "output exists")) + if logger: + logger.info("Skip existing extraction: %s", output_dir) + continue + shutil.rmtree(output_dir) + + output_dir.mkdir(parents=True, exist_ok=True) + if archive.name.lower().endswith(".tar.gz"): + _extract_tar_gz(archive, output_dir) + elif archive.suffix.lower() == ".zip": + _extract_zip(archive, output_dir) + else: + continue + + results.append(ExtractionResult(archive, output_dir, "done")) + if logger: + logger.info("Extracted %s -> %s", archive, output_dir) + except Exception as exc: + results.append(ExtractionResult(archive, output_dir, "failed", str(exc))) + if logger: + logger.exception("Failed to extract %s", archive) + + return results + + +def _archive_output_dir(archive: Path) -> Path: + if archive.name.lower().endswith(".tar.gz"): + return archive.with_name(archive.name[:-7]) + return archive.with_suffix("") + + +def _extract_tar_gz(archive: Path, output_dir: Path) -> None: + base = output_dir.resolve() + with tarfile.open(archive, "r:gz") as tar: + members = tar.getmembers() + for member in tqdm(members, desc=f"Extract {archive.name}", leave=False, unit="file"): + if member.issym() or member.islnk(): + raise ProcessingError(f"Refusing to extract link from archive: {member.name}") + _ensure_safe_target(base, member.name) + tar.extract(member, path=base) + + +def _extract_zip(archive: Path, output_dir: Path) -> None: + base = output_dir.resolve() + with zipfile.ZipFile(archive, "r") as zip_ref: + names = zip_ref.namelist() + for name in tqdm(names, desc=f"Extract {archive.name}", leave=False, unit="file"): + _ensure_safe_target(base, name) + zip_ref.extract(name, base) + + +def _ensure_safe_target(base: Path, member_name: str) -> None: + target = (base / member_name).resolve() + if target == base: + return + if base not in target.parents: + raise ProcessingError(f"Unsafe archive path: {member_name}") diff --git a/src/gf3_l1a2l2/discovery.py b/src/gf3_l1a2l2/discovery.py new file mode 100644 index 0000000..264dc8c --- /dev/null +++ b/src/gf3_l1a2l2/discovery.py @@ -0,0 +1,90 @@ +from __future__ import annotations + +import logging +from pathlib import Path +from typing import Iterable, Optional + +from gf3_l1a2l2.errors import DiscoveryError +from gf3_l1a2l2.models import SceneFiles +from gf3_l1a2l2.paths import detect_polarization + +SKIP_DIR_NAMES = {".git", ".venv", ".tmp_tests", "__pycache__", ".pytest_cache", "logs"} + + +def discover_scenes(input_dir: Path | str, logger: Optional[logging.Logger] = None) -> list[SceneFiles]: + root_dir = Path(input_dir).resolve() + if not root_dir.is_dir(): + raise DiscoveryError(f"Input directory does not exist: {root_dir}") + + scenes: list[SceneFiles] = [] + for directory in _walk_directories(root_dir): + try: + files = [path for path in directory.iterdir() if path.is_file()] + except OSError as exc: + if logger: + logger.warning("Skipping unreadable directory %s: %s", directory, exc) + continue + + metadata_files = sorted(path for path in files if path.name.lower().endswith(".meta.xml")) + image_files = sorted( + path + for path in files + if path.suffix.lower() in (".tif", ".tiff") and detect_polarization(path) + ) + + if not metadata_files or not image_files: + continue + + metadata_xml = metadata_files[0] + if len(metadata_files) > 1 and logger: + logger.warning("Multiple metadata files in %s; using %s", directory, metadata_xml.name) + + images = _collect_by_polarization(image_files, logger) + rpcs = _collect_by_polarization( + sorted(path for path in files if path.suffix.lower() in (".rpc", ".rpb")), + logger, + ) + + scenes.append(SceneFiles(root=directory, metadata_xml=metadata_xml, images=images, rpcs=rpcs)) + + scenes.sort(key=lambda scene: str(scene.root).lower()) + if logger: + logger.info("Discovered %s scene(s) under %s", len(scenes), root_dir) + return scenes + + +def _walk_directories(root_dir: Path) -> Iterable[Path]: + stack = [root_dir] + while stack: + directory = stack.pop() + if directory.name in SKIP_DIR_NAMES: + continue + yield directory + try: + children = list(directory.iterdir()) + except OSError: + continue + stack.extend(path for path in children if path.is_dir() and path.name not in SKIP_DIR_NAMES) + + +def _collect_by_polarization( + files: Iterable[Path], + logger: Optional[logging.Logger] = None, +) -> dict[str, Path]: + collected: dict[str, Path] = {} + for path in files: + polarization = detect_polarization(path) + if not polarization: + continue + if polarization in collected: + if logger: + logger.warning( + "Duplicate %s file in %s; using %s and ignoring %s", + polarization, + path.parent, + collected[polarization].name, + path.name, + ) + continue + collected[polarization] = path + return collected diff --git a/src/gf3_l1a2l2/errors.py b/src/gf3_l1a2l2/errors.py new file mode 100644 index 0000000..291c1da --- /dev/null +++ b/src/gf3_l1a2l2/errors.py @@ -0,0 +1,18 @@ +class GF3PipelineError(Exception): + """Base error for pipeline failures.""" + + +class DiscoveryError(GF3PipelineError): + """Raised when input products cannot be discovered.""" + + +class MetadataError(GF3PipelineError): + """Raised when GF3 metadata cannot be parsed.""" + + +class RpcError(GF3PipelineError): + """Raised when RPC metadata cannot be parsed.""" + + +class ProcessingError(GF3PipelineError): + """Raised when image processing fails.""" diff --git a/src/gf3_l1a2l2/geocorrection.py b/src/gf3_l1a2l2/geocorrection.py new file mode 100644 index 0000000..9514098 --- /dev/null +++ b/src/gf3_l1a2l2/geocorrection.py @@ -0,0 +1,78 @@ +from __future__ import annotations + +import logging +from pathlib import Path +from typing import Optional + +from gf3_l1a2l2.errors import ProcessingError +from gf3_l1a2l2.paths import output_name +from gf3_l1a2l2.rpc import read_rpc_metadata + + +def geometric_correction( + l1b_path: Path | str, + rpc_path: Path | str, + output_dir: Path | str, + *, + dem_path: Optional[Path], + x_res: float, + y_res: float, + dst_srs: str, + overwrite: bool = False, + logger: Optional[logging.Logger] = None, +) -> Path: + source = Path(l1b_path) + out_dir = Path(output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + l2_path = out_dir / output_name(source, "L1B", "L2") + + if l2_path.exists() and not overwrite: + if logger: + logger.info("Skip existing L2: %s", l2_path) + return l2_path + + if dem_path and not Path(dem_path).is_file(): + raise ProcessingError(f"DEM file does not exist: {dem_path}") + + rpc_metadata = read_rpc_metadata(rpc_path) + + try: + from osgeo import gdal + except Exception as exc: + raise ProcessingError("Cannot import osgeo.gdal. Install GDAL Python bindings.") from exc + + if logger: + logger.info("Geocorrecting %s -> %s", source, l2_path) + + dataset = gdal.Open(str(source), gdal.GA_Update) + if dataset is None: + raise ProcessingError(f"GDAL cannot open L1B image: {source}") + + try: + dataset.SetMetadata(rpc_metadata, "RPC") + transformer_options = [] + if dem_path: + transformer_options.append(f"RPC_DEM={Path(dem_path)}") + + warp_options = gdal.WarpOptions( + dstSRS=dst_srs, + xRes=x_res, + yRes=y_res, + rpc=True, + multithread=True, + transformerOptions=transformer_options, + creationOptions=["COMPRESS=LZW", "TILED=YES", "BIGTIFF=YES"], + ) + result = gdal.Warp(str(l2_path), dataset, options=warp_options) + if result is None: + raise ProcessingError(f"GDAL Warp failed for {source}") + result.FlushCache() + result = None + except ProcessingError: + raise + except Exception as exc: + raise ProcessingError(f"Failed to geocorrect {source}") from exc + finally: + dataset = None + + return l2_path diff --git a/src/gf3_l1a2l2/log.py b/src/gf3_l1a2l2/log.py new file mode 100644 index 0000000..40616d8 --- /dev/null +++ b/src/gf3_l1a2l2/log.py @@ -0,0 +1,34 @@ +from __future__ import annotations + +import logging +from datetime import datetime +from pathlib import Path + + +def configure_logging(output_dir: Path, level: str = "INFO") -> logging.Logger: + output_dir.mkdir(parents=True, exist_ok=True) + log_dir = output_dir / "logs" + log_dir.mkdir(parents=True, exist_ok=True) + + logger = logging.getLogger("gf3_l1a2l2") + logger.setLevel(getattr(logging, level.upper(), logging.INFO)) + logger.propagate = False + logger.handlers.clear() + + formatter = logging.Formatter( + fmt="%(asctime)s %(levelname)s %(message)s", + datefmt="%Y-%m-%d %H:%M:%S", + ) + + console = logging.StreamHandler() + console.setFormatter(formatter) + console.setLevel(getattr(logging, level.upper(), logging.INFO)) + logger.addHandler(console) + + timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") + file_handler = logging.FileHandler(log_dir / f"gf3_l1a2l2_{timestamp}.log", encoding="utf-8") + file_handler.setFormatter(formatter) + file_handler.setLevel(logging.DEBUG) + logger.addHandler(file_handler) + + return logger diff --git a/src/gf3_l1a2l2/metadata.py b/src/gf3_l1a2l2/metadata.py new file mode 100644 index 0000000..31e8334 --- /dev/null +++ b/src/gf3_l1a2l2/metadata.py @@ -0,0 +1,145 @@ +from __future__ import annotations + +import math +import xml.etree.ElementTree as ET +from pathlib import Path +from typing import Iterable, Optional + +from gf3_l1a2l2.constants import POLARIZATIONS +from gf3_l1a2l2.errors import MetadataError +from gf3_l1a2l2.models import CalibrationMetadata, CalibrationValue +from gf3_l1a2l2.paths import detect_polarization + + +QUALIFY_KEYWORDS = ("qualifyvalue", "qualityvalue") +CALIBRATION_KEYWORDS = ("calibrationconst", "calibrationconstant") + + +def read_calibration_metadata(xml_path: Path | str) -> CalibrationMetadata: + path = Path(xml_path) + try: + root = ET.parse(path).getroot() + except ET.ParseError as exc: + raise MetadataError(f"Invalid metadata XML: {path}") from exc + except OSError as exc: + raise MetadataError(f"Cannot read metadata XML: {path}") from exc + + qualify_values = _find_values(root, QUALIFY_KEYWORDS) or _legacy_values(root, (17, 13)) + calibration_values = _find_values(root, CALIBRATION_KEYWORDS) or _legacy_values(root, (18, 3)) + + if not qualify_values: + raise MetadataError(f"Cannot find QualifyValue values in {path}") + if not calibration_values: + raise MetadataError(f"Cannot find CalibrationConst values in {path}") + + values: dict[str, CalibrationValue] = {} + for polarization in POLARIZATIONS: + if polarization not in qualify_values or polarization not in calibration_values: + continue + values[polarization] = CalibrationValue( + qualify_value=qualify_values[polarization], + calibration_const=calibration_values[polarization], + ) + + if not values: + raise MetadataError(f"No usable calibration values found in {path}") + + return CalibrationMetadata(values=values) + + +def _find_values(root: ET.Element, keywords: Iterable[str]) -> Optional[dict[str, float]]: + lowered = tuple(keyword.lower() for keyword in keywords) + for element in root.iter(): + name = _local_name(element.tag).lower() + if not any(keyword in name for keyword in lowered): + continue + + values = _values_by_child_polarization(element) + if values: + return values + + values = _values_by_child_order(element) + if values: + return values + + values = _values_from_text(element.text) + if values: + return values + + return None + + +def _legacy_values(root: ET.Element, indexes: tuple[int, int]) -> Optional[dict[str, float]]: + try: + node = root + for index in indexes: + node = list(node)[index] + values = [_parse_float(child.text) for child in list(node)[:4]] + except (IndexError, TypeError, ValueError): + return None + + if len(values) < 4: + return None + return dict(zip(POLARIZATIONS, values)) + + +def _values_by_child_polarization(element: ET.Element) -> Optional[dict[str, float]]: + values: dict[str, float] = {} + for child in list(element): + polarization = _element_polarization(child) + if not polarization: + continue + try: + values[polarization] = _parse_float(child.text) + except ValueError: + continue + return values if values else None + + +def _values_by_child_order(element: ET.Element) -> Optional[dict[str, float]]: + numeric_values: list[float] = [] + for child in list(element): + try: + numeric_values.append(_parse_float(child.text)) + except ValueError: + continue + if len(numeric_values) < 4: + return None + return dict(zip(POLARIZATIONS, numeric_values[:4])) + + +def _values_from_text(text: Optional[str]) -> Optional[dict[str, float]]: + if not text: + return None + normalized = text.replace(",", " ").replace(";", " ") + parts = [part for part in normalized.split() if part] + if len(parts) < 4: + return None + try: + values = [_parse_float(part) for part in parts[:4]] + except ValueError: + return None + return dict(zip(POLARIZATIONS, values)) + + +def _element_polarization(element: ET.Element) -> Optional[str]: + for value in element.attrib.values(): + polarization = detect_polarization(value) + if polarization: + return polarization + return detect_polarization(_local_name(element.tag)) + + +def _parse_float(text: Optional[str]) -> float: + if text is None: + raise ValueError("missing value") + value = text.strip().strip("=") + if not value or value.upper() == "NULL": + return math.nan + return float(value) + + +def _local_name(tag: str) -> str: + if "}" in tag: + return tag.rsplit("}", 1)[1] + return tag diff --git a/src/gf3_l1a2l2/models.py b/src/gf3_l1a2l2/models.py new file mode 100644 index 0000000..5295117 --- /dev/null +++ b/src/gf3_l1a2l2/models.py @@ -0,0 +1,66 @@ +from __future__ import annotations + +from dataclasses import dataclass, field +from pathlib import Path +from typing import Mapping, Optional + + +@dataclass(frozen=True) +class CalibrationValue: + qualify_value: float + calibration_const: float + + +@dataclass(frozen=True) +class CalibrationMetadata: + values: Mapping[str, CalibrationValue] + + def for_polarization(self, polarization: str) -> CalibrationValue: + return self.values[polarization] + + +@dataclass(frozen=True) +class SceneFiles: + root: Path + metadata_xml: Path + images: Mapping[str, Path] + rpcs: Mapping[str, Path] + + @property + def name(self) -> str: + return self.root.name + + @property + def available_polarizations(self) -> tuple[str, ...]: + return tuple(sorted(self.images)) + + +@dataclass(frozen=True) +class TaskResult: + scene: str + polarization: str + status: str + l1b_path: Optional[Path] = None + l2_path: Optional[Path] = None + message: str = "" + + +@dataclass +class PipelineSummary: + results: list[TaskResult] = field(default_factory=list) + + @property + def succeeded(self) -> int: + return sum(1 for item in self.results if item.status == "done") + + @property + def skipped(self) -> int: + return sum(1 for item in self.results if item.status == "skipped") + + @property + def failed(self) -> int: + return sum(1 for item in self.results if item.status == "failed") + + @property + def total(self) -> int: + return len(self.results) diff --git a/src/gf3_l1a2l2/paths.py b/src/gf3_l1a2l2/paths.py new file mode 100644 index 0000000..37d7ba2 --- /dev/null +++ b/src/gf3_l1a2l2/paths.py @@ -0,0 +1,43 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Optional + +from gf3_l1a2l2.constants import POLARIZATIONS + + +def detect_polarization(path_or_name: Path | str) -> Optional[str]: + name = Path(path_or_name).name.upper() + for polarization in POLARIZATIONS: + if polarization in name: + return polarization + return None + + +def output_name(source: Path, source_level: str, target_level: str) -> str: + name = source.name + index = name.lower().find(source_level.lower()) + if index >= 0: + return name[:index] + target_level + name[index + len(source_level) :] + + stem = source.stem + suffix = "".join(source.suffixes) + if suffix and stem.endswith(suffix): + stem = stem[: -len(suffix)] + return f"{stem}_{target_level}{source.suffix}" + + +def scene_output_dir(input_dir: Path, output_dir: Path, scene_dir: Path, flat: bool) -> Path: + if flat: + return output_dir + + try: + relative = scene_dir.relative_to(input_dir) + except ValueError: + relative = Path(scene_dir.name) + + if str(relative) in ("", "."): + relative = Path(scene_dir.name) + + safe_name = "_".join(part for part in relative.parts if part not in ("", ".")) + return output_dir / safe_name diff --git a/src/gf3_l1a2l2/pipeline.py b/src/gf3_l1a2l2/pipeline.py new file mode 100644 index 0000000..3a53f58 --- /dev/null +++ b/src/gf3_l1a2l2/pipeline.py @@ -0,0 +1,182 @@ +from __future__ import annotations + +import concurrent.futures +import logging +from pathlib import Path +from typing import Iterable + +from gf3_l1a2l2.calibration import radiometric_calibrate_l1a_to_l1b +from gf3_l1a2l2.config import PipelineConfig +from gf3_l1a2l2.constants import POLARIZATIONS +from gf3_l1a2l2.decompress import extract_compressed_files +from gf3_l1a2l2.discovery import discover_scenes +from gf3_l1a2l2.errors import GF3PipelineError +from gf3_l1a2l2.geocorrection import geometric_correction +from gf3_l1a2l2.log import configure_logging +from gf3_l1a2l2.metadata import read_calibration_metadata +from gf3_l1a2l2.models import PipelineSummary, SceneFiles, TaskResult +from gf3_l1a2l2.paths import scene_output_dir + + +def run_pipeline(config: PipelineConfig) -> PipelineSummary: + config = config.resolved() + logger = configure_logging(config.output_dir, config.log_level) + _validate_config(config) + + logger.info("Input: %s", config.input_dir) + logger.info("Output: %s", config.output_dir) + if config.dem_path: + logger.info("DEM: %s", config.dem_path) + + if not config.skip_decompress: + extract_compressed_files(config.input_dir, overwrite=config.overwrite, logger=logger) + + scenes = discover_scenes(config.input_dir, logger=logger) + if not scenes: + logger.warning("No GF3 scenes found under %s", config.input_dir) + return PipelineSummary() + + summary = PipelineSummary() + for results in _process_scenes(scenes, config, logger): + summary.results.extend(results) + + logger.info( + "Pipeline finished: total=%s done=%s skipped=%s failed=%s", + summary.total, + summary.succeeded, + summary.skipped, + summary.failed, + ) + return summary + + +def _validate_config(config: PipelineConfig) -> None: + if not config.input_dir.is_dir(): + raise GF3PipelineError(f"Input directory does not exist: {config.input_dir}") + if config.dem_path and not config.dem_path.is_file(): + raise GF3PipelineError(f"DEM file does not exist: {config.dem_path}") + if config.x_res <= 0 or config.y_res <= 0: + raise GF3PipelineError("x_res and y_res must be positive") + + +def _process_scenes( + scenes: list[SceneFiles], + config: PipelineConfig, + logger: logging.Logger, +) -> Iterable[list[TaskResult]]: + if config.workers == 1: + for scene in scenes: + yield _process_scene(scene, config, logger) + return + + with concurrent.futures.ThreadPoolExecutor(max_workers=config.workers) as executor: + future_map = { + executor.submit(_process_scene, scene, config, logger): scene + for scene in scenes + } + for future in concurrent.futures.as_completed(future_map): + scene = future_map[future] + try: + yield future.result() + except Exception as exc: + logger.exception("Scene failed: %s", scene.root) + if config.fail_fast: + raise + yield [ + TaskResult( + scene=scene.name, + polarization=polarization, + status="failed", + message=str(exc), + ) + for polarization in scene.available_polarizations + ] + + +def _process_scene( + scene: SceneFiles, + config: PipelineConfig, + logger: logging.Logger, +) -> list[TaskResult]: + logger.info("Processing scene: %s", scene.root) + out_dir = scene_output_dir(config.input_dir, config.output_dir, scene.root, config.flat_output) + out_dir.mkdir(parents=True, exist_ok=True) + + try: + metadata = read_calibration_metadata(scene.metadata_xml) + except Exception as exc: + logger.exception("Failed to parse metadata for %s", scene.root) + if config.fail_fast: + raise + return [ + TaskResult(scene=scene.name, polarization=pol, status="failed", message=str(exc)) + for pol in scene.available_polarizations + ] + + results: list[TaskResult] = [] + for polarization in POLARIZATIONS: + image_path = scene.images.get(polarization) + if not image_path: + continue + + rpc_path = scene.rpcs.get(polarization) + if not rpc_path: + message = f"Missing RPC/RPB file for {polarization}" + logger.error("%s in %s", message, scene.root) + results.append(TaskResult(scene=scene.name, polarization=polarization, status="failed", message=message)) + if config.fail_fast: + raise GF3PipelineError(message) + continue + + if polarization not in metadata.values: + message = f"Missing calibration metadata for {polarization}" + logger.error("%s in %s", message, scene.metadata_xml) + results.append(TaskResult(scene=scene.name, polarization=polarization, status="failed", message=message)) + if config.fail_fast: + raise GF3PipelineError(message) + continue + + try: + calibration = metadata.for_polarization(polarization) + l1b_path = radiometric_calibrate_l1a_to_l1b( + image_path, + out_dir, + calibration, + block_size=config.block_size, + overwrite=config.overwrite, + logger=logger, + ) + l2_path = geometric_correction( + l1b_path, + rpc_path, + out_dir, + dem_path=config.dem_path, + x_res=config.x_res, + y_res=config.y_res, + dst_srs=config.dst_srs, + overwrite=config.overwrite, + logger=logger, + ) + results.append( + TaskResult( + scene=scene.name, + polarization=polarization, + status="done", + l1b_path=l1b_path, + l2_path=l2_path, + ) + ) + except Exception as exc: + logger.exception("Failed %s %s", scene.name, polarization) + results.append( + TaskResult( + scene=scene.name, + polarization=polarization, + status="failed", + message=str(exc), + ) + ) + if config.fail_fast: + raise + + return results diff --git a/src/gf3_l1a2l2/rpc.py b/src/gf3_l1a2l2/rpc.py new file mode 100644 index 0000000..e0b0ab2 --- /dev/null +++ b/src/gf3_l1a2l2/rpc.py @@ -0,0 +1,70 @@ +from __future__ import annotations + +import re +from pathlib import Path + +from gf3_l1a2l2.errors import RpcError + + +RPC_KEY_MAP = { + "errBias": "ERR_BIAS", + "errRand": "ERR_RAND", + "lineOffset": "LINE_OFF", + "sampOffset": "SAMP_OFF", + "latOffset": "LAT_OFF", + "longOffset": "LONG_OFF", + "heightOffset": "HEIGHT_OFF", + "lineScale": "LINE_SCALE", + "sampScale": "SAMP_SCALE", + "latScale": "LAT_SCALE", + "longScale": "LONG_SCALE", + "heightScale": "HEIGHT_SCALE", + "lineNumCoef": "LINE_NUM_COEFF", + "lineDenCoef": "LINE_DEN_COEFF", + "sampNumCoef": "SAMP_NUM_COEFF", + "sampDenCoef": "SAMP_DEN_COEFF", +} + +COEFFICIENT_KEYS = { + "lineNumCoef", + "lineDenCoef", + "sampNumCoef", + "sampDenCoef", +} + + +def read_rpc_metadata(rpc_path: Path | str) -> dict[str, str]: + path = Path(rpc_path) + try: + text = path.read_text(encoding="utf-8", errors="ignore") + except OSError as exc: + raise RpcError(f"Cannot read RPC file: {path}") from exc + + metadata: dict[str, str] = {} + missing: list[str] = [] + for source_key, target_key in RPC_KEY_MAP.items(): + value = _find_rpc_value(text, source_key) + if value is None: + missing.append(source_key) + continue + metadata[target_key] = _normalize_rpc_value(source_key, value) + + if missing: + raise RpcError(f"RPC file {path} is missing fields: {', '.join(missing)}") + + return metadata + + +def _find_rpc_value(text: str, key: str) -> str | None: + pattern = re.compile(rf"\b{re.escape(key)}\b\s*[:=]?\s*(.*?);", re.IGNORECASE | re.DOTALL) + match = pattern.search(text) + if not match: + return None + return match.group(1) + + +def _normalize_rpc_value(key: str, raw_value: str) -> str: + value = raw_value.strip().strip("=") + if key in COEFFICIENT_KEYS: + value = value.replace("(", " ").replace(")", " ").replace(",", " ") + return " ".join(value.split()) diff --git a/tests/fixtures/sample.meta.xml b/tests/fixtures/sample.meta.xml new file mode 100644 index 0000000..fe77b3d --- /dev/null +++ b/tests/fixtures/sample.meta.xml @@ -0,0 +1,8 @@ + + + 1.02.0NULL4.0 + + + 10.020.030.040.0 + + diff --git a/tests/fixtures/sample.rpb b/tests/fixtures/sample.rpb new file mode 100644 index 0000000..ccb7b40 --- /dev/null +++ b/tests/fixtures/sample.rpb @@ -0,0 +1,16 @@ +errBias = -1; +errRand = -1; +lineOffset = 10; +sampOffset = 20; +latOffset = 30; +longOffset = 40; +heightOffset = 50; +lineScale = 60; +sampScale = 70; +latScale = 80; +longScale = 90; +heightScale = 100; +lineNumCoef = (1, 2, 3); +lineDenCoef = (4, 5, 6); +sampNumCoef = (7, 8, 9); +sampDenCoef = (10, 11, 12); diff --git a/tests/test_metadata.py b/tests/test_metadata.py new file mode 100644 index 0000000..91b81cb --- /dev/null +++ b/tests/test_metadata.py @@ -0,0 +1,19 @@ +import math +import unittest +from pathlib import Path + +from gf3_l1a2l2.metadata import read_calibration_metadata + + +class MetadataTests(unittest.TestCase): + def test_reads_named_calibration_values(self): + path = Path(__file__).parent / "fixtures" / "sample.meta.xml" + metadata = read_calibration_metadata(path) + + self.assertEqual(metadata.for_polarization("HH").qualify_value, 1.0) + self.assertEqual(metadata.for_polarization("HV").calibration_const, 20.0) + self.assertTrue(math.isnan(metadata.for_polarization("VH").qualify_value)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_rpc.py b/tests/test_rpc.py new file mode 100644 index 0000000..04c1527 --- /dev/null +++ b/tests/test_rpc.py @@ -0,0 +1,18 @@ +import unittest +from pathlib import Path + +from gf3_l1a2l2.rpc import read_rpc_metadata + + +class RpcTests(unittest.TestCase): + def test_reads_rpc_metadata(self): + path = Path(__file__).parent / "fixtures" / "sample.rpb" + metadata = read_rpc_metadata(path) + + self.assertEqual(metadata["ERR_BIAS"], "-1") + self.assertEqual(metadata["LINE_OFF"], "10") + self.assertEqual(metadata["LINE_NUM_COEFF"], "1 2 3") + + +if __name__ == "__main__": + unittest.main()