""" GF3 (GaoFen-3) L1A -> L2 processing service. Pipeline: Extract archive -> Parse XML metadata -> Radiometric calibration -> RPC geometric correction. Pure Python implementation using GDAL/numpy, no ENVI/SARscape dependency. """ from __future__ import annotations import json import logging import math import os import shutil import tarfile import zipfile from typing import Any, Dict, List, Optional, Tuple import numpy as np logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- # Archive extraction # --------------------------------------------------------------------------- def _extract_archive(path: str, dest: str) -> str: """Extract .tar.gz or .zip archive to *dest*, return extracted directory path.""" os.makedirs(dest, exist_ok=True) if tarfile.is_tarfile(path): with tarfile.open(path, "r:*") as tf: # Security: prevent path traversal for member in tf.getmembers(): if member.name.startswith("/") or ".." in member.name: raise ValueError(f"Unsafe path in archive: {member.name}") tf.extractall(dest) elif zipfile.is_zipfile(path): with zipfile.ZipFile(path, "r") as zf: for info in zf.infolist(): if info.filename.startswith("/") or ".." in info.filename: raise ValueError(f"Unsafe path in archive: {info.filename}") zf.extractall(dest) else: raise ValueError(f"Unsupported archive format: {path}") return dest # --------------------------------------------------------------------------- # GF3 XML metadata parser # --------------------------------------------------------------------------- def _parse_gf3_meta(xml_path: str) -> Dict[str, Any]: """Parse GF3 product XML metadata. Returns dict with keys per polarization: polarizations: list of str (e.g. ["HH", "VV"]) calibration: {pol: {"QualifyValue": float, "CalibrationConst": float}} """ import defusedxml.ElementTree as ET tree = ET.parse(xml_path) root = tree.getroot() result: Dict[str, Any] = {"polarizations": [], "calibration": {}} # Find all imageinfo or channel elements # GF3 XML structure varies; search by tag name for robustness def _find_all_recursive(element, tag): found = [] for child in element.iter(): if tag.lower() in child.tag.lower(): found.append(child) return found # Try to find QualifyValue and CalibrationConst # Typical GF3 XML has -> and # per polarization channel qualify_values = {} cal_consts = {} # Search for elements containing polarization-specific calibration info for elem in root.iter(): tag = elem.tag.split("}")[-1] if "}" in elem.tag else elem.tag # strip namespace if tag == "QualifyValue" and elem.text: # Parent should indicate which polarization this belongs to parent = _find_parent(root, elem) pol = _extract_polarization(parent, elem) if pol: try: qualify_values[pol] = float(elem.text.strip()) except ValueError: pass if tag == "CalibrationConst" and elem.text: parent = _find_parent(root, elem) pol = _extract_polarization(parent, elem) if pol: try: cal_consts[pol] = float(elem.text.strip()) except ValueError: pass # If per-polarization search didn't work, try flat extraction if not qualify_values: # Fallback: find all QualifyValue elements in order qv_elems = [e for e in root.iter() if e.tag.split("}")[-1] == "QualifyValue" and e.text] cc_elems = [e for e in root.iter() if e.tag.split("}")[-1] == "CalibrationConst" and e.text] # Find polarization list pol_elems = [e for e in root.iter() if e.tag.split("}")[-1] == "Polarisation" and e.text] if not pol_elems: pol_elems = [e for e in root.iter() if e.tag.split("}")[-1] == "polarization" and e.text] pols = [e.text.strip().upper() for e in pol_elems] for i, pol in enumerate(pols): if i < len(qv_elems): try: qualify_values[pol] = float(qv_elems[i].text.strip()) except ValueError: pass if i < len(cc_elems): try: cal_consts[pol] = float(cc_elems[i].text.strip()) except ValueError: pass polarizations = sorted(set(list(qualify_values.keys()) + list(cal_consts.keys()))) if not polarizations: # Last resort: guess from TIFF filenames in same directory xml_dir = os.path.dirname(xml_path) for f in os.listdir(xml_dir): fl = f.upper() for pol in ("HH", "HV", "VH", "VV"): if pol in fl and f.lower().endswith((".tif", ".tiff")): if pol not in polarizations: polarizations.append(pol) polarizations.sort() calibration = {} for pol in polarizations: calibration[pol] = { "QualifyValue": qualify_values.get(pol, 1.0), "CalibrationConst": cal_consts.get(pol, 0.0), } result["polarizations"] = polarizations result["calibration"] = calibration logger.info("[GF3] Parsed metadata: polarizations=%s, calibration=%s", polarizations, calibration) return result def _find_parent(root, target): """Find the parent element of *target* in the tree.""" for parent in root.iter(): for child in parent: if child is target: return parent return None def _extract_polarization(parent, elem) -> Optional[str]: """Try to extract polarization from context around an XML element.""" if parent is None: return None # Check parent tag or sibling elements tag = parent.tag.split("}")[-1] if "}" in parent.tag else parent.tag for pol in ("HH", "HV", "VH", "VV"): if pol in tag.upper(): return pol # Check sibling text for child in parent: child_tag = child.tag.split("}")[-1] if "}" in child.tag else child.tag if "polaris" in child_tag.lower() and child.text: return child.text.strip().upper() return None # --------------------------------------------------------------------------- # RPC file parser # --------------------------------------------------------------------------- def _read_rpb(rpb_path: str) -> Dict[str, Any]: """Parse a .rpb (Rational Polynomial Coefficients) file into GDAL metadata dict.""" rpc = {} current_key = None values = [] with open(rpb_path, "r") as f: for line in f: line = line.strip().rstrip(";") if "=" in line: if current_key and values: rpc[current_key] = values values = [] key, _, val = line.partition("=") current_key = key.strip() val = val.strip().strip("(").strip(")") if val: for v in val.replace(",", " ").split(): try: values.append(float(v)) except ValueError: values.append(v) else: # continuation of values val = line.strip("()").strip() if val: for v in val.replace(",", " ").split(): try: values.append(float(v)) except ValueError: values.append(v) if current_key and values: rpc[current_key] = values # Map to GDAL RPC metadata keys gdal_rpc = {} key_map = { "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", } for rpb_key, gdal_key in key_map.items(): if rpb_key in rpc: val = rpc[rpb_key] if isinstance(val, list): if len(val) == 1: gdal_rpc[gdal_key] = str(val[0]) else: gdal_rpc[gdal_key] = " ".join(str(v) for v in val) else: gdal_rpc[gdal_key] = str(val) return gdal_rpc # --------------------------------------------------------------------------- # Radiometric calibration # --------------------------------------------------------------------------- def _radiometric_calibration(tiff_path: str, qv: float, cal: float, output_path: str) -> str: """L1A -> L1B single-polarization radiometric calibration. Formula: A = sqrt(I^2 + Q^2), dB = 20*log10(A * QV / 65535) - Cal For amplitude-only TIFF: dB = 20*log10(A * QV / 65535) - Cal """ from osgeo import gdal ds = gdal.Open(tiff_path, gdal.GA_ReadOnly) if ds is None: raise FileNotFoundError(f"Cannot open TIFF: {tiff_path}") n_bands = ds.RasterCount width = ds.RasterXSize height = ds.RasterYSize if n_bands >= 2: # Complex I/Q data band_i = ds.GetRasterBand(1).ReadAsArray().astype(np.float64) band_q = ds.GetRasterBand(2).ReadAsArray().astype(np.float64) amplitude = np.sqrt(band_i ** 2 + band_q ** 2) else: # Amplitude only amplitude = ds.GetRasterBand(1).ReadAsArray().astype(np.float64) # Avoid log of zero amplitude = np.where(amplitude > 0, amplitude, np.nan) db_values = 20.0 * np.log10(amplitude * qv / 65535.0) - cal db_values = np.where(np.isfinite(db_values), db_values, 0).astype(np.float32) # Write calibrated result preserving georeference driver = gdal.GetDriverByName("GTiff") out_ds = driver.Create(output_path, width, height, 1, gdal.GDT_Float32, options=["COMPRESS=DEFLATE"]) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) # Copy RPC metadata if present rpc_md = ds.GetMetadata("RPC") if rpc_md: out_ds.SetMetadata(rpc_md, "RPC") out_ds.GetRasterBand(1).WriteArray(db_values) out_ds.GetRasterBand(1).SetNoDataValue(0) out_ds.FlushCache() out_ds = None ds = None logger.info("[GF3] Calibration done: %s -> %s (QV=%.2f, Cal=%.2f)", tiff_path, output_path, qv, cal) return output_path # --------------------------------------------------------------------------- # Geometric correction # --------------------------------------------------------------------------- def _geometric_correction( l1b_path: str, rpb_path: Optional[str], output_path: str, resolution: float, dem_path: str, ) -> str: """L1B -> L2 RPC geometric correction using GDAL Warp.""" from osgeo import gdal ds = gdal.Open(l1b_path, gdal.GA_ReadOnly) if ds is None: raise FileNotFoundError(f"Cannot open L1B: {l1b_path}") # If RPC not already in dataset, load from .rpb rpc_md = ds.GetMetadata("RPC") if not rpc_md and rpb_path and os.path.isfile(rpb_path): rpc_md = _read_rpb(rpb_path) ds.SetMetadata(rpc_md, "RPC") logger.info("[GF3] Loaded RPC from %s", rpb_path) warp_options = gdal.WarpOptions( dstSRS="EPSG:4326", format="GTiff", xRes=resolution, yRes=resolution, rpc=True, creationOptions=["COMPRESS=DEFLATE"], ) # Use DEM if available if dem_path and os.path.isfile(dem_path): warp_options = gdal.WarpOptions( dstSRS="EPSG:4326", format="GTiff", xRes=resolution, yRes=resolution, rpc=True, transformerOptions=[f"RPC_DEM={dem_path}"], creationOptions=["COMPRESS=DEFLATE"], ) result = gdal.Warp(output_path, ds, options=warp_options) ds = None if result is None: raise RuntimeError(f"GDAL Warp failed for {l1b_path}") result = None logger.info("[GF3] Geometric correction done: %s -> %s", l1b_path, output_path) return output_path # --------------------------------------------------------------------------- # Main pipeline # --------------------------------------------------------------------------- def run_gf3_l1a_to_l2( input_dir: str, output_dir: str, resolution: float = 0.0002, job_id: Optional[str] = None, ) -> Dict[str, Any]: """Run GF3 L1A -> L2 pipeline: extract -> calibrate -> geometric correction. Args: input_dir: Path to GF3 L1A product directory (or archive file) output_dir: Output directory for L2 products resolution: Output resolution in degrees (default 0.0002 ~ 20m) job_id: Optional job ID for progress tracking Returns: dict with keys: ok, l2_paths, polarizations, output_dir """ from ..config import settings dem_path = settings.GF3_GEO_DEM_PATH os.makedirs(output_dir, exist_ok=True) # Step 1: Extract if archive work_dir = input_dir if os.path.isfile(input_dir): logger.info("[GF3] Extracting archive: %s", input_dir) extract_dir = os.path.join(output_dir, "_extracted") _extract_archive(input_dir, extract_dir) # Find actual data directory (may be nested) subdirs = [d for d in os.listdir(extract_dir) if os.path.isdir(os.path.join(extract_dir, d))] work_dir = os.path.join(extract_dir, subdirs[0]) if subdirs else extract_dir # Step 2: Scan for XML + TIFF + RPB files xml_path = None tiff_files: Dict[str, str] = {} # pol -> tiff path rpb_files: Dict[str, str] = {} # pol -> rpb path for f in os.listdir(work_dir): fl = f.lower() fp = os.path.join(work_dir, f) if fl.endswith(".meta.xml") or (fl.endswith(".xml") and "meta" in fl): xml_path = fp elif fl.endswith(".xml") and xml_path is None: xml_path = fp elif fl.endswith((".tif", ".tiff")): for pol in ("HH", "HV", "VH", "VV"): if pol in f.upper(): tiff_files[pol] = fp break elif fl.endswith(".rpb"): for pol in ("HH", "HV", "VH", "VV"): if pol in f.upper(): rpb_files[pol] = fp break if not xml_path: # Try subdirectories for sub in os.listdir(work_dir): sub_path = os.path.join(work_dir, sub) if os.path.isdir(sub_path): for f in os.listdir(sub_path): fl = f.lower() fp = os.path.join(sub_path, f) if fl.endswith(".xml") and not xml_path: xml_path = fp elif fl.endswith((".tif", ".tiff")): for pol in ("HH", "HV", "VH", "VV"): if pol in f.upper(): tiff_files[pol] = fp break elif fl.endswith(".rpb"): for pol in ("HH", "HV", "VH", "VV"): if pol in f.upper(): rpb_files[pol] = fp break if not tiff_files: return {"ok": False, "error": f"No TIFF files found in {work_dir}"} # Step 3: Parse XML metadata meta = {"polarizations": list(tiff_files.keys()), "calibration": {}} if xml_path: try: meta = _parse_gf3_meta(xml_path) except Exception as e: logger.warning("[GF3] Failed to parse XML %s: %s, using defaults", xml_path, e) polarizations = meta.get("polarizations", list(tiff_files.keys())) calibration = meta.get("calibration", {}) # Step 4-5: Process each polarization l2_paths = [] for pol in polarizations: if pol not in tiff_files: logger.warning("[GF3] No TIFF found for polarization %s, skipping", pol) continue tiff_path = tiff_files[pol] cal_info = calibration.get(pol, {"QualifyValue": 1.0, "CalibrationConst": 0.0}) qv = cal_info["QualifyValue"] cal = cal_info["CalibrationConst"] # L1A -> L1B (calibration) l1b_path = os.path.join(output_dir, f"{pol}_L1B.tif") logger.info("[GF3] Calibrating %s (QV=%.2f, Cal=%.2f)", pol, qv, cal) _radiometric_calibration(tiff_path, qv, cal, l1b_path) # L1B -> L2 (geometric correction) l2_path = os.path.join(output_dir, f"{pol}_L2.tif") rpb_path = rpb_files.get(pol) logger.info("[GF3] Geometric correction %s (resolution=%.6f)", pol, resolution) _geometric_correction(l1b_path, rpb_path, l2_path, resolution, dem_path) l2_paths.append(l2_path) # Clean up intermediate L1B try: os.remove(l1b_path) except OSError: pass if not l2_paths: return {"ok": False, "error": "No polarization channels processed successfully"} logger.info("[GF3] Pipeline complete: %d L2 products", len(l2_paths)) return { "ok": True, "l2_paths": l2_paths, "polarizations": polarizations, "output_dir": output_dir, "input_dir_name": os.path.basename(input_dir), } # --------------------------------------------------------------------------- # Auto-register L2 result into radar_data table # --------------------------------------------------------------------------- async def register_l2_to_radar_data(l2_dir: str, input_dir_name: str, polarizations: List[str], db) -> Optional[int]: """Register a GF3 L2 output directory as a radar_data record. Args: l2_dir: Path to the L2 output directory input_dir_name: Original L1A input directory name (for metadata extraction) polarizations: List of polarization channels processed db: AsyncSession Returns: radar_data.id if successfully registered, else None """ from ..utils import parse_gf3_l2_dirname from ..models import RadarDataORM from .data_service import extract_geotiff_bounds from geoalchemy2.shape import from_shape from shapely.geometry import Polygon as ShapelyPolygon from sqlalchemy.future import select # Check if already registered (by file_path) result = await db.execute( select(RadarDataORM).where(RadarDataORM.file_path == l2_dir).limit(1) ) if result.scalar_one_or_none(): logger.info("[GF3] L2 dir already registered: %s", l2_dir) return None # Parse metadata from the input directory name meta = parse_gf3_l2_dirname(input_dir_name) if not meta: # Fallback: minimal metadata meta = { "satellite": "GF3", "imaging_date": None, "polarization": ",".join(polarizations) if polarizations else None, } # Try to extract polygon from the first L2 GeoTIFF polygon = None try: for f in os.listdir(l2_dir): if f.lower().endswith((".tif", ".tiff")) and "L2" in f: tiff_path = os.path.join(l2_dir, f) polygon = extract_geotiff_bounds(tiff_path) if polygon: break except OSError: pass coverage_geom = None if polygon and len(polygon) >= 4: try: shp = ShapelyPolygon(polygon) if shp.is_valid: coverage_geom = from_shape(shp, srid=4326) except Exception: pass radar = RadarDataORM( file_path=l2_dir, satellite=meta.get("satellite", "GF3"), imaging_date=meta.get("imaging_date"), imaging_mode=meta.get("imaging_mode"), polarization=meta.get("polarization") or (",".join(polarizations) if polarizations else None), scene_center_lon=meta.get("scene_center_lon"), scene_center_lat=meta.get("scene_center_lat"), coverage_polygon=coverage_geom, ) db.add(radar) await db.flush() radar_id = radar.id await db.commit() logger.info("[GF3] Registered L2 in radar_data: id=%s, path=%s", radar_id, l2_dir) return radar_id