Refactor GF3 L1A to L2 pipeline
This commit is contained in:
@@ -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)
|
||||
@@ -0,0 +1,5 @@
|
||||
from gf3_l1a2l2.cli import main
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
raise SystemExit(main())
|
||||
@@ -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)
|
||||
@@ -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
|
||||
@@ -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)),
|
||||
)
|
||||
@@ -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
|
||||
@@ -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}")
|
||||
@@ -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
|
||||
@@ -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."""
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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)
|
||||
@@ -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
|
||||
@@ -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
|
||||
@@ -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())
|
||||
Reference in New Issue
Block a user