diff --git a/pixi.toml b/pixi.toml index b67e92ef..1eecc38d 100644 --- a/pixi.toml +++ b/pixi.toml @@ -117,7 +117,7 @@ lint = "ruff check python/siac tests" format-check = "ruff format --check python/siac tests" format = "ruff format python/siac tests" typecheck = "mypy python/siac" -typecheck-scoped = "mypy --follow-imports=skip python/siac/config/algorithms.py python/siac/config/load.py python/siac/config/providers.py python/siac/config/public.py python/siac/config/request.py python/siac/config/resolved.py python/siac/config/snapshot.py python/siac/config/system.py python/siac/config/types.py python/siac/adapters/auth.py python/siac/adapters/bestpixel.py python/siac/adapters/bestpixel_surface_prior.py python/siac/adapters/atmo/cams.py python/siac/adapters/atmo/maiac_day_aod.py python/siac/adapters/atmo/mcd19_earthaccess.py python/siac/adapters/atmo/merra2.py python/siac/adapters/brdf/mcd43_earthaccess.py python/siac/adapters/brdf/_mcd43_kernels.py python/siac/adapters/brdf/_mcd43_qa.py python/siac/adapters/data/copernicus_dataspace.py python/siac/adapters/data/gcs_sentinel2.py python/siac/adapters/data/earthaccess_catalog.py python/siac/adapters/earthdata.py python/siac/adapters/satellite/sentinel2.py python/siac/algorithms/brdf/kernels.py python/siac/algorithms/cloud/providers/omnicloudmask.py python/siac/algorithms/rt/lut/backend.py python/siac/algorithms/rt/lut/_compact.py python/siac/algorithms/rt/lut/_spectral.py python/siac/algorithms/rt/lut/http_zip_store.py python/siac/algorithms/rt/lut/store.py python/siac/algorithms/rt/direct/libradtran.py python/siac/algorithms/rt/direct/libradtran_runner.py python/siac/algorithms/rt/direct/libradtran_build.py python/siac/algorithms/rt/direct/_build_common.py python/siac/algorithms/rt/direct/_common.py python/siac/algorithms/rt/_run_cache.py python/siac/workflows/scene_setup.py python/siac/algorithms/surface/_swir_refine_resample.py python/siac/algorithms/surface/brdf_monthly_database.py python/siac/algorithms/surface/brdf_whittaker.py python/siac/algorithms/surface/kernel_model.py python/siac/algorithms/grid/toa_psf.py python/siac/algorithms/surface/spectral_mapping.py python/siac/algorithms/surface/swir_refine.py python/siac/storage/writers.py python/siac/storage/readers.py python/siac/geo/reprojection.py python/siac/api/public.py python/siac/app/_assembly_correction.py python/siac/app/_assembly_io.py python/siac/app/_assembly_preprocessor.py python/siac/app/_assembly_providers.py python/siac/app/_assembly_rt.py python/siac/app/_assembly_solver.py python/siac/app/_assembly_surface.py python/siac/app/planning.py python/siac/app/s2_backend.py python/siac/app/sentinel2.py python/siac/workflows/_pipeline_config.py python/siac/workflows/scene.py python/siac/workflows/sentinel2.py" +typecheck-scoped = "mypy --follow-imports=skip python/siac/config/algorithms.py python/siac/config/load.py python/siac/config/providers.py python/siac/config/public.py python/siac/config/request.py python/siac/config/resolved.py python/siac/config/snapshot.py python/siac/config/system.py python/siac/config/types.py python/siac/adapters/auth.py python/siac/adapters/bestpixel.py python/siac/adapters/bestpixel_surface_prior.py python/siac/adapters/atmo/cams.py python/siac/adapters/atmo/maiac_day_aod.py python/siac/adapters/atmo/mcd19_earthaccess.py python/siac/adapters/atmo/merra2.py python/siac/adapters/brdf/mcd43_earthaccess.py python/siac/adapters/brdf/_mcd43_kernels.py python/siac/adapters/brdf/_mcd43_qa.py python/siac/adapters/data/copernicus_dataspace.py python/siac/adapters/data/gcs_sentinel2.py python/siac/adapters/data/earthaccess_catalog.py python/siac/adapters/earthdata.py python/siac/adapters/satellite/sentinel2.py python/siac/adapters/satellite/landsat.py python/siac/algorithms/brdf/kernels.py python/siac/algorithms/cloud/providers/omnicloudmask.py python/siac/algorithms/rt/lut/backend.py python/siac/algorithms/rt/lut/_compact.py python/siac/algorithms/rt/lut/_spectral.py python/siac/algorithms/rt/lut/http_zip_store.py python/siac/algorithms/rt/lut/store.py python/siac/algorithms/rt/direct/libradtran.py python/siac/algorithms/rt/direct/libradtran_runner.py python/siac/algorithms/rt/direct/libradtran_build.py python/siac/algorithms/rt/direct/_build_common.py python/siac/algorithms/rt/direct/_common.py python/siac/algorithms/rt/_run_cache.py python/siac/workflows/scene_setup.py python/siac/algorithms/surface/_swir_refine_resample.py python/siac/algorithms/surface/brdf_monthly_database.py python/siac/algorithms/surface/brdf_whittaker.py python/siac/algorithms/surface/kernel_model.py python/siac/algorithms/grid/toa_psf.py python/siac/algorithms/surface/spectral_mapping.py python/siac/algorithms/surface/swir_refine.py python/siac/storage/writers.py python/siac/storage/readers.py python/siac/geo/reprojection.py python/siac/api/public.py python/siac/app/_assembly_correction.py python/siac/app/_assembly_io.py python/siac/app/_assembly_preprocessor.py python/siac/app/_assembly_providers.py python/siac/app/_assembly_rt.py python/siac/app/_assembly_solver.py python/siac/app/_assembly_surface.py python/siac/app/planning.py python/siac/app/s2_backend.py python/siac/app/sentinel2.py python/siac/workflows/_pipeline_config.py python/siac/workflows/scene.py python/siac/workflows/sentinel2.py" build-rust = "maturin develop --release --skip-install --manifest-path src/siac_rs/Cargo.toml" rust-test = "cargo test --manifest-path src/siac_rs/Cargo.toml" # Invoke cargo-clippy directly (resolved via PATH -> the pixi env's clippy) diff --git a/python/siac/adapters/satellite/__init__.py b/python/siac/adapters/satellite/__init__.py index 3aa5eb0a..2c0dfa66 100644 --- a/python/siac/adapters/satellite/__init__.py +++ b/python/siac/adapters/satellite/__init__.py @@ -47,4 +47,5 @@ ) # Import sensor-specific preprocessors to trigger registration +from siac.adapters.satellite.landsat import LandsatPreprocessor as LandsatPreprocessor from siac.adapters.satellite.sentinel2 import Sentinel2Preprocessor as Sentinel2Preprocessor diff --git a/python/siac/adapters/satellite/base.py b/python/siac/adapters/satellite/base.py index 335267c4..aea337af 100644 --- a/python/siac/adapters/satellite/base.py +++ b/python/siac/adapters/satellite/base.py @@ -397,11 +397,20 @@ def detect_sensor(input_path: str | Path) -> str: return "s2" # Check for Landsat MTL - mtl_files = list(input_path.glob("*_MTL.txt")) + list(input_path.glob("*_MTL.xml")) + if input_path.is_file() and input_path.name.upper().endswith(("_MTL.TXT", "_MTL.JSON")): + mtl_files = [input_path] + else: + mtl_files = [ + *input_path.glob("*_MTL.txt"), + *input_path.glob("*_MTL.json"), + *input_path.glob("*_MTL.xml"), + ] if mtl_files: # Check if Landsat 8/9 or older mtl_content = mtl_files[0].read_text() - if "LANDSAT_8" in mtl_content or "LANDSAT_9" in mtl_content: + if "LANDSAT_9" in mtl_content: + return "l9" + if "LANDSAT_8" in mtl_content: return "l8" elif "LANDSAT_5" in mtl_content or "LANDSAT_7" in mtl_content: return "l5" diff --git a/python/siac/adapters/satellite/landsat.py b/python/siac/adapters/satellite/landsat.py new file mode 100644 index 00000000..2b519aed --- /dev/null +++ b/python/siac/adapters/satellite/landsat.py @@ -0,0 +1,718 @@ +""" +Landsat 8/9 OLI Collection-2 Level-1 preprocessor. + +Reads an unpacked USGS Collection-2 Level-1 scene (``*_MTL.txt`` or +``*_MTL.json`` plus the per-band GeoTIFFs) into the M1 contract: + +- TOA reflectance for the reflective OLI bands ``B1``-``B7``, rescaled with the + per-band ``REFLECTANCE_MULT/ADD`` coefficients from the MTL and corrected for + the per-pixel solar zenith angle from the Collection-2 ``SZA`` angle band + (scene-centre ``SUN_ELEVATION`` when the angle band is absent). +- Sun and view geometry from the Collection-2 ``SZA``/``SAA``/``VZA``/``VAA`` + angle bands, block-averaged to a coarse grid (the Sentinel-2 path keeps the + 5 km metadata grid; every downstream consumer resamples to its own grid). +- Standardized cloud classes, either from the configured cloud-mask provider or, + with ``provider = "product_qa"``, from the product's own ``QA_PIXEL`` band. +""" + +from __future__ import annotations + +import json +import logging +import re +from datetime import datetime +from pathlib import Path +from typing import TYPE_CHECKING, Any, cast + +import numpy as np +import xarray as xr + +from siac.adapters.rsrf import load_sensor_config_with_rsrf +from siac.adapters.satellite.base import ( + BaseSatellitePreprocessor, + degrees_to_radians, + register_preprocessor, +) +from siac.algorithms.cloud import build_cloud_classes, classes_to_bool_mask +from siac.geo import reproject_match +from siac.runtime import GeometryAngles +from siac.storage.readers import read_raster, read_raster_window + +if TYPE_CHECKING: + from siac.domain import SensorConfig + +logger = logging.getLogger(__name__) +_TOA_BAND_LOADER_ATTR = "_siac_toa_band_loader" + +_SATELLITE_BY_SPACECRAFT = {"LANDSAT_8": "L8", "LANDSAT_9": "L9"} + +# Collection-2 angle bands are int16 in hundredths of a degree. +_ANGLE_SCALE = 0.01 +_ANGLE_SUFFIXES = { + "sza": ("FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4", "_SZA.TIF"), + "saa": ("FILE_NAME_ANGLE_SOLAR_AZIMUTH_BAND_4", "_SAA.TIF"), + "vza": ("FILE_NAME_ANGLE_SENSOR_ZENITH_BAND_4", "_VZA.TIF"), + "vaa": ("FILE_NAME_ANGLE_SENSOR_AZIMUTH_BAND_4", "_VAA.TIF"), +} +_QA_PIXEL_KEY = ("FILE_NAME_QUALITY_L1_PIXEL", "_QA_PIXEL.TIF") + +# Collection-2 QA_PIXEL bit positions. +_QA_FILL = 0 +_QA_DILATED_CLOUD = 1 +_QA_CIRRUS = 2 +_QA_CLOUD = 3 +_QA_CLOUD_SHADOW = 4 + +_DEFAULT_ANGLE_BLOCK_PIXELS = 100 # 3 km blocks at 30 m + + +def parse_mtl_text(text: str) -> dict[str, Any]: + """Parse a Landsat ``*_MTL.txt`` ODL document into nested dictionaries. + + Quoted values stay strings; unquoted values are returned as the raw token, + leaving numeric conversion to the caller (MTL mixes dates, times, integers + and floats under the same syntax). + """ + root: dict[str, Any] = {} + stack: list[dict[str, Any]] = [root] + for raw_line in text.splitlines(): + line = raw_line.strip() + if not line or line == "END": + continue + if "=" not in line: + continue + key, value = (part.strip() for part in line.split("=", 1)) + if key == "GROUP": + group: dict[str, Any] = {} + stack[-1][value] = group + stack.append(group) + elif key == "END_GROUP": + if len(stack) > 1: + stack.pop() + else: + if len(value) >= 2 and value[0] == value[-1] == '"': + value = value[1:-1] + stack[-1][key] = value + return root + + +def flatten_mtl(tree: dict[str, Any]) -> dict[str, str]: + """Flatten nested MTL groups to a ``{KEY: value}`` map (first occurrence wins).""" + flat: dict[str, str] = {} + + def _walk(node: dict[str, Any]) -> None: + for key, value in node.items(): + if isinstance(value, dict): + _walk(value) + elif key not in flat: + flat[key] = str(value) + + _walk(tree) + return flat + + +def qa_pixel_to_cloud_classes(qa: np.ndarray) -> np.ndarray: + """Map Collection-2 ``QA_PIXEL`` bits to SIAC classes {0,1,2,3}. + + Fill -> 0 (missing); cloud, dilated cloud or high-confidence cirrus -> 2; + cloud shadow -> 3; everything else (clear, water, snow) -> 1. + """ + bits = np.asarray(qa).astype(np.uint16) + + def _bit(position: int) -> np.ndarray: + return cast("np.ndarray", ((bits >> np.uint16(position)) & np.uint16(1)).astype(bool)) + + classes: np.ndarray = np.ones(bits.shape, dtype=np.uint8) + classes[_bit(_QA_CLOUD_SHADOW)] = 3 + cloud = _bit(_QA_CLOUD) | _bit(_QA_DILATED_CLOUD) | _bit(_QA_CIRRUS) + classes[cloud] = 2 + classes[_bit(_QA_FILL)] = 0 + return classes + + +@register_preprocessor("landsat") +@register_preprocessor("l9") +@register_preprocessor("l8") +class LandsatPreprocessor(BaseSatellitePreprocessor): + """Preprocessor for Landsat 8/9 OLI Collection-2 Level-1 scenes.""" + + def __init__(self, config: dict[str, Any] | None = None): + super().__init__(config) + self._satellite_id: str | None = None + self._sensor_config_cache: SensorConfig | None = None + self._sensor_config_cache_satellite_id: str | None = None + self._scene_dir: Path | None = None + self._mtl: dict[str, str] | None = None + self._resolved_input_path: Path | None = None + self._reference_grid: xr.DataArray | None = None + self._cos_sza: xr.DataArray | None = None + self._last_cloud_classes: xr.DataArray | None = None + self._subset_bounds: tuple[float, float, float, float] | None = None + self._subset_crs: str | None = None + + # ------------------------------------------------------------------ + # Spatial subset + # ------------------------------------------------------------------ + + def set_spatial_subset( + self, + bounds: tuple[float, float, float, float], + *, + crs: str, + ) -> None: + """Restrict raster reads to an AOI window.""" + xmin, ymin, xmax, ymax = bounds + self._subset_bounds = (float(xmin), float(ymin), float(xmax), float(ymax)) + self._subset_crs = str(crs) + self._reference_grid = None + self._cos_sza = None + + def clear_spatial_subset(self) -> None: + """Clear any active AOI read window.""" + self._subset_bounds = None + self._subset_crs = None + self._reference_grid = None + self._cos_sza = None + + # ------------------------------------------------------------------ + # Sensor configuration + # ------------------------------------------------------------------ + + @property + def sensor_config(self) -> SensorConfig: + """Return the RSRF-backed OLI configuration for the resolved platform.""" + satellite_id = self._satellite_id or "L8" + if ( + self._sensor_config_cache is not None + and self._sensor_config_cache_satellite_id == satellite_id + ): + return self._sensor_config_cache + try: + sensor_config = load_sensor_config_with_rsrf( + "OLI", + satellite_id, + rsrf_root=self.config.get("rsrf_root"), + ) + except Exception as exc: + raise RuntimeError( + f"Unable to load Landsat RSRF metadata for {satellite_id}. Configure " + f"paths.rsrf_root to a catalog containing {satellite_id} OLI band responses." + ) from exc + self._sensor_config_cache = sensor_config + self._sensor_config_cache_satellite_id = satellite_id + return sensor_config + + @property + def band_names(self) -> tuple[str, ...]: + return tuple(band.name for band in self.sensor_config.bands) + + # ------------------------------------------------------------------ + # Metadata + # ------------------------------------------------------------------ + + def get_metadata(self, input_path: str | Path) -> dict[str, Any]: + """Extract observation metadata from the scene MTL.""" + input_path = Path(input_path).expanduser() + mtl = self._resolve_scene(input_path) + + metadata: dict[str, Any] = { + "sensor": "OLI", + "satellite": self._satellite_id, + "input_path": str(input_path), + "product_id": mtl.get("LANDSAT_PRODUCT_ID"), + "processing_level": mtl.get("PROCESSING_LEVEL"), + "collection_number": mtl.get("COLLECTION_NUMBER"), + "collection_category": mtl.get("COLLECTION_CATEGORY"), + } + path, row = mtl.get("WRS_PATH"), mtl.get("WRS_ROW") + if path is not None and row is not None: + metadata["wrs_path"] = int(path) + metadata["wrs_row"] = int(row) + metadata["tile_id"] = f"{int(path):03d}{int(row):03d}" + for key, name in ( + ("SUN_AZIMUTH", "sun_azimuth"), + ("SUN_ELEVATION", "sun_elevation"), + ("EARTH_SUN_DISTANCE", "earth_sun_distance"), + ("CLOUD_COVER", "cloud_cover"), + ): + if key in mtl: + metadata[name] = float(mtl[key]) + + observation_time = self._parse_observation_time(mtl) + if observation_time is None: + observation_time = self._parse_observation_time_from_name( + mtl.get("LANDSAT_PRODUCT_ID") or input_path.name + ) + if observation_time is not None: + metadata["observation_time"] = observation_time + + metadata["reflectance_mult"] = {} + metadata["reflectance_add"] = {} + for band_name in self.band_names: + number = self._band_number(band_name) + mult = mtl.get(f"REFLECTANCE_MULT_BAND_{number}") + add = mtl.get(f"REFLECTANCE_ADD_BAND_{number}") + if mult is not None and add is not None: + metadata["reflectance_mult"][band_name] = float(mult) + metadata["reflectance_add"][band_name] = float(add) + return metadata + + @staticmethod + def _parse_observation_time(mtl: dict[str, str]) -> datetime | None: + date_text = mtl.get("DATE_ACQUIRED") + time_text = mtl.get("SCENE_CENTER_TIME") + if not date_text: + return None + if not time_text: + return datetime.strptime(date_text, "%Y-%m-%d") + # SCENE_CENTER_TIME carries up to 7 fractional digits; datetime takes 6. + match = re.match(r"(\d{2}):(\d{2}):(\d{2})(?:\.(\d+))?", time_text.strip().rstrip("Z")) + if match is None: + logger.warning("Could not parse Landsat SCENE_CENTER_TIME %r", time_text) + return datetime.strptime(date_text, "%Y-%m-%d") + hour, minute, second, fraction = match.groups() + microsecond = int((fraction or "0")[:6].ljust(6, "0")) + day = datetime.strptime(date_text, "%Y-%m-%d") + return day.replace( + hour=int(hour), minute=int(minute), second=int(second), microsecond=microsecond + ) + + @staticmethod + def _parse_observation_time_from_name(name: str) -> datetime | None: + """Infer the acquisition date from a Collection-2 product id.""" + match = re.search(r"_\d{6}_(\d{8})_\d{8}_", name) + if match is None: + return None + try: + return datetime.strptime(match.group(1), "%Y%m%d") + except ValueError: + return None + + # ------------------------------------------------------------------ + # TOA reflectance + # ------------------------------------------------------------------ + + def load_toa( + self, + input_path: str | Path, + *, + band_names: tuple[str, ...] | list[str] | None = None, + ) -> xr.Dataset: + """Load TOA reflectance for the reflective OLI bands.""" + input_path = Path(input_path).expanduser() + metadata = self.get_metadata(input_path) + requested = self._resolve_requested_band_names(band_names) + + reference_name = "B4" if "B4" in self.band_names else self.band_names[0] + reference = self._load_band(reference_name, metadata=metadata, reference_grid=None) + self._reference_grid = reference + + loaded: dict[str, xr.DataArray] = {} + if reference_name in requested: + loaded[reference_name] = reference + for band_name in requested: + if band_name not in loaded: + loaded[band_name] = self._load_band( + band_name, metadata=metadata, reference_grid=reference + ) + + ds = xr.Dataset(loaded) + ds.attrs[_TOA_BAND_LOADER_ATTR] = self._make_band_loader( + metadata=metadata, + reference_grid=reference, + loaded_band_names=tuple(loaded), + ) + ds.attrs["loaded_toa_bands"] = tuple(loaded) + ds.attrs["sensor"] = "OLI" + ds.attrs["satellite"] = self._satellite_id + obs_time = metadata.get("observation_time") + ds.attrs["observation_time"] = obs_time.isoformat() if obs_time else "" + return ds + + def _resolve_requested_band_names( + self, band_names: tuple[str, ...] | list[str] | None + ) -> tuple[str, ...]: + available = self.band_names + requested = tuple(str(name) for name in band_names) if band_names else available + resolved: list[str] = [] + for name in requested: + if name not in available: + logger.warning("Ignoring unsupported Landsat band request %s", name) + continue + if self._band_file(name) is None: + logger.warning("Band %s not found", name) + continue + if name not in resolved: + resolved.append(name) + if not resolved: + raise RuntimeError("No requested Landsat bands were available for TOA loading.") + return tuple(resolved) + + def _load_band( + self, + band_name: str, + *, + metadata: dict[str, Any], + reference_grid: xr.DataArray | None, + ) -> xr.DataArray: + band_file = self._band_file(band_name) + if band_file is None: + raise KeyError(f"Band {band_name!r} is not available in the Landsat scene.") + mult = metadata["reflectance_mult"].get(band_name) + add = metadata["reflectance_add"].get(band_name) + if mult is None or add is None: + raise ValueError( + f"Landsat MTL lacks REFLECTANCE_MULT/ADD for {band_name}; " + "Collection-2 Level-1 metadata is required." + ) + dn = self._read(band_file) + if reference_grid is not None and not self._coords_match(dn, reference_grid): + dn = reproject_match(dn, reference_grid, resampling="nearest") + reflectance = self._calibrate(dn, mult=mult, add=add, metadata=metadata) + reflectance.name = band_name + return reflectance + + def _calibrate( + self, + dn: xr.DataArray, + *, + mult: float, + add: float, + metadata: dict[str, Any], + ) -> xr.DataArray: + """DN -> TOA reflectance: (M*Q + A) / cos(SZA); DN 0 is fill.""" + values = dn.astype(np.float32) + fill = ~np.isfinite(values) | (values == 0) + rho_prime = values * np.float32(mult) + np.float32(add) + cos_sza = self._cos_solar_zenith(dn, metadata=metadata) + rho = (rho_prime / cos_sza).clip(np.float32(0.0), np.float32(1.5)) + return cast("xr.DataArray", rho.where(~fill).astype(np.float32)) + + def _cos_solar_zenith( + self, template: xr.DataArray, *, metadata: dict[str, Any] + ) -> xr.DataArray | np.float32: + """Per-pixel cos(SZA) from the angle band, else scene-centre sun elevation.""" + cached = self._cos_sza + if cached is not None and self._coords_match(cached, template): + return cached + sza_file = self._auxiliary_file(*_ANGLE_SUFFIXES["sza"]) + if sza_file is not None: + sza = self._read(sza_file).astype(np.float32) * np.float32(_ANGLE_SCALE) + if not self._coords_match(sza, template): + sza = reproject_match(sza, template, resampling="bilinear") + cos_sza = np.cos(np.deg2rad(sza)).astype(np.float32) + cos_sza = cast("xr.DataArray", cos_sza.where(cos_sza > 0)) + self._cos_sza = cos_sza + return cos_sza + elevation = metadata.get("sun_elevation") + if elevation is None: + raise ValueError("Landsat scene has neither an SZA angle band nor SUN_ELEVATION.") + logger.warning( + "Landsat SZA angle band not found; using scene-centre SUN_ELEVATION for the " + "TOA solar-zenith correction." + ) + return np.float32(np.sin(np.deg2rad(float(elevation))).item()) + + def _make_band_loader( + self, + *, + metadata: dict[str, Any], + reference_grid: xr.DataArray, + loaded_band_names: tuple[str, ...], + ) -> Any: + cache: dict[tuple[str, bool], xr.DataArray] = {} + loaded = set(loaded_band_names) + + def _load(band_name: str, *, native: bool = False) -> xr.DataArray: + key = (band_name, native) + if key in cache: + return cache[key] + if band_name in loaded and not native: + raise KeyError(f"Band {band_name!r} is already present in the loaded TOA dataset.") + band = self._load_band( + band_name, + metadata=metadata, + reference_grid=None if native else reference_grid, + ) + cache[key] = band + return band + + return _load + + # ------------------------------------------------------------------ + # Geometry + # ------------------------------------------------------------------ + + def extract_geometry(self, input_path: str | Path) -> GeometryAngles: + """Block-average the Collection-2 angle bands to a coarse geometry grid.""" + input_path = Path(input_path).expanduser() + self._resolve_scene(input_path) + block = int(self.config.get("angle_block_pixels", _DEFAULT_ANGLE_BLOCK_PIXELS)) + if block < 1: + raise ValueError(f"angle_block_pixels must be >= 1, got {block}") + + grids: dict[str, xr.DataArray] = {} + for name, (mtl_key, suffix) in _ANGLE_SUFFIXES.items(): + path = self._auxiliary_file(mtl_key, suffix) + if path is None: + raise FileNotFoundError( + f"Landsat angle band {suffix} not found for {input_path}. " + "Collection-2 Level-1 scenes ship SZA/SAA/VZA/VAA bands." + ) + raw = self._read(path).astype(np.float32) + degrees = raw * np.float32(_ANGLE_SCALE) + # Angle bands carry 0 outside the imaged footprint; the zero test is + # safe for SZA but not for VZA/VAA, so footprint comes from QA/TOA. + degrees = degrees.where(self._footprint(raw)) + grids[name] = self._block_reduce(degrees, block, circular=name in ("saa", "vaa")) + + return GeometryAngles( + sza=degrees_to_radians(grids["sza"]), + saa=degrees_to_radians(grids["saa"]), + vza=degrees_to_radians(grids["vza"]), + vaa=degrees_to_radians(grids["vaa"]), + ) + + def _footprint(self, template: xr.DataArray) -> xr.DataArray: + qa_file = self._auxiliary_file(*_QA_PIXEL_KEY) + if qa_file is not None: + qa = self._read(qa_file) + if not self._coords_match(qa, template): + qa = reproject_match(qa, template, resampling="nearest") + valid = np.isfinite(qa) & (qa_pixel_to_cloud_classes(np.nan_to_num(qa.values)) != 0) + return xr.DataArray(np.asarray(valid), dims=template.dims, coords=template.coords) + if self._reference_grid is not None and self._coords_match(self._reference_grid, template): + return cast("xr.DataArray", np.isfinite(self._reference_grid)) + return cast("xr.DataArray", np.isfinite(template) & (template != 0)) + + @staticmethod + def _block_reduce(da: xr.DataArray, block: int, *, circular: bool) -> xr.DataArray: + """NaN-aware block mean (circular for azimuths), then fill empty blocks.""" + from scipy.ndimage import distance_transform_edt + + values = np.asarray(da.values, dtype=np.float64) + height, width = values.shape + rows, cols = -(-height // block), -(-width // block) + padded = np.full((rows * block, cols * block), np.nan) + padded[:height, :width] = values + tiles = padded.reshape(rows, block, cols, block).transpose(0, 2, 1, 3) + tiles = tiles.reshape(rows, cols, block * block) + finite = np.isfinite(tiles) + count = finite.sum(axis=-1) + if circular: + radians = np.deg2rad(np.where(finite, tiles, 0.0)) + sine = np.where(finite, np.sin(radians), 0.0).sum(axis=-1) + cosine = np.where(finite, np.cos(radians), 0.0).sum(axis=-1) + reduced = np.mod(np.rad2deg(np.arctan2(sine, cosine)), 360.0) + else: + total = np.where(finite, tiles, 0.0).sum(axis=-1) + reduced = np.divide(total, np.maximum(count, 1)) + reduced = np.where(count > 0, reduced, np.nan) + if not np.isfinite(reduced).any(): + raise ValueError("Landsat angle band has no valid pixels inside the scene footprint.") + empty = ~np.isfinite(reduced) + if empty.any(): + _, (iy, ix) = distance_transform_edt(empty, return_indices=True) + reduced = reduced[iy, ix] + if circular: + unwrapped = np.unwrap(np.unwrap(np.deg2rad(reduced), axis=1), axis=0) + reduced = np.rad2deg(unwrapped) + + x = np.asarray(da.coords["x"].values, dtype=np.float64) + y = np.asarray(da.coords["y"].values, dtype=np.float64) + dx = float(x[1] - x[0]) if x.size > 1 else 30.0 + dy = float(y[1] - y[0]) if y.size > 1 else -30.0 + offset = (block - 1) / 2.0 + x_coarse = x[0] + (np.arange(cols) * block + offset) * dx + y_coarse = y[0] + (np.arange(rows) * block + offset) * dy + out = xr.DataArray( + reduced.astype(np.float32), + dims=("y", "x"), + coords={"y": y_coarse, "x": x_coarse}, + ) + crs = da.rio.crs + if crs is not None: + out = out.rio.write_crs(crs) + return out + + # ------------------------------------------------------------------ + # Cloud mask + # ------------------------------------------------------------------ + + def extract_cloud_mask( + self, + input_path: str | Path, + toa: xr.Dataset | None = None, + ) -> xr.DataArray: + """Build standardized cloud classes and return the boolean mask.""" + input_path = Path(input_path).expanduser() + if toa is None: + toa = self.load_toa(input_path) + settings = self._cloud_mask_settings(input_path) + + if settings["mode"] == "auto" and settings["provider"] == "product_qa": + cloud_classes = self._product_qa_classes(toa) + else: + cloud_classes = build_cloud_classes( + toa, + self.sensor_config, + mode=settings["mode"], + provider=settings["provider"], + class_mapping=settings["class_mapping"], + external_mask_path=settings["external_mask_path"], + user_callable=settings["user_callable"], + target_resolution_m=settings["target_resolution_m"], + resolution_policy=settings["resolution_policy"], + allow_upsample_to_target=settings["allow_upsample_to_target"], + unmapped_to_missing=settings["unmapped_to_missing"], + cloud_cache_dir=settings.get("cache_dir"), + inference_device=settings.get("inference_device"), + ) + self._last_cloud_classes = cloud_classes + return classes_to_bool_mask(cloud_classes) + + def _product_qa_classes(self, toa: xr.Dataset) -> xr.DataArray: + qa_file = self._auxiliary_file(*_QA_PIXEL_KEY) + if qa_file is None: + raise FileNotFoundError( + "cloud_mask.provider='product_qa' needs the Collection-2 QA_PIXEL band." + ) + reference = next(iter(toa.data_vars.values())) + qa = self._read(qa_file) + if not self._coords_match(qa, reference): + qa = reproject_match(qa, reference, resampling="nearest") + classes = qa_pixel_to_cloud_classes(np.nan_to_num(np.asarray(qa.values), nan=1.0)) + classes[~np.isfinite(np.asarray(reference.values))] = 0 + out = xr.DataArray( + classes, dims=reference.dims, coords=reference.coords, name="cloud_classes" + ) + out.attrs["class_values"] = {0: "missing", 1: "clear", 2: "cloud", 3: "shadow"} + out.attrs["source"] = "landsat_c2_qa_pixel" + return out + + def _cloud_mask_settings(self, input_path: Path) -> dict[str, Any]: + defaults: dict[str, Any] = { + "mode": "auto", + "provider": "omnicloudmask", + "class_mapping": None, + "external_mask_path": None, + "user_callable": None, + "target_resolution_m": 30.0, + "resolution_policy": "auto", + "allow_upsample_to_target": False, + "unmapped_to_missing": True, + "inference_device": "cpu", + } + cloud_mask_config = ( + self.config.get("cloud_mask", {}) if isinstance(self.config, dict) else {} + ) + if isinstance(cloud_mask_config, dict): + defaults.update(cloud_mask_config) + for key in ("mode", "provider", "resolution_policy"): + defaults[key] = str(getattr(defaults[key], "value", defaults[key])) + external = defaults.get("external_mask_path") + if isinstance(external, (str, Path)) and str(external): + p = Path(external).expanduser() + if not p.is_absolute(): + p = (self._scene_dir or input_path) / p + defaults["external_mask_path"] = p + return defaults + + # ------------------------------------------------------------------ + # Scene resolution and I/O + # ------------------------------------------------------------------ + + def _resolve_scene(self, input_path: Path) -> dict[str, str]: + if self._resolved_input_path == input_path and self._mtl is not None: + return self._mtl + if input_path.is_file(): + mtl_path = input_path + scene_dir = input_path.parent + else: + scene_dir = input_path + candidates = sorted(scene_dir.glob("*_MTL.txt")) or sorted(scene_dir.glob("*_MTL.json")) + if not candidates: + raise FileNotFoundError(f"No Landsat *_MTL.txt or *_MTL.json under {scene_dir}") + if len(candidates) > 1: + raise ValueError( + f"Multiple Landsat MTL files under {scene_dir}: " + + ", ".join(path.name for path in candidates) + ) + mtl_path = candidates[0] + + text = mtl_path.read_text(encoding="utf-8") + tree = json.loads(text) if mtl_path.suffix.lower() == ".json" else parse_mtl_text(text) + mtl = flatten_mtl(tree) + + spacecraft = mtl.get("SPACECRAFT_ID", "").upper() + satellite_id = _SATELLITE_BY_SPACECRAFT.get(spacecraft) + if satellite_id is None: + raise ValueError( + f"Unsupported Landsat spacecraft {spacecraft!r}; only Landsat 8/9 OLI " + "Collection-2 Level-1 scenes are supported." + ) + level = mtl.get("PROCESSING_LEVEL", "") + if level and not level.startswith("L1"): + raise ValueError( + f"Landsat product processing level is {level!r}; SIAC expects a Level-1 " + "(L1TP/L1GT/L1GS) scene, not surface reflectance." + ) + + self._satellite_id = satellite_id + self._scene_dir = scene_dir + self._mtl = mtl + self._resolved_input_path = input_path + self._reference_grid = None + self._cos_sza = None + return mtl + + def _band_file(self, band_name: str) -> Path | None: + number = self._band_number(band_name) + return self._auxiliary_file(f"FILE_NAME_BAND_{number}", f"_B{number}.TIF") + + def _auxiliary_file(self, mtl_key: str, suffix: str) -> Path | None: + if self._scene_dir is None or self._mtl is None: + raise RuntimeError("Landsat scene is not resolved; call get_metadata() first.") + named = self._mtl.get(mtl_key) + if named: + candidate = self._scene_dir / named + if candidate.exists(): + return candidate + matches = [ + path + for path in sorted(self._scene_dir.iterdir()) + if path.name.upper().endswith(suffix.upper()) + ] + return matches[0] if matches else None + + @staticmethod + def _band_number(band_name: str) -> int: + match = re.fullmatch(r"B(\d{1,2})", band_name) + if match is None: + raise ValueError(f"Not a Landsat OLI band name: {band_name!r}") + return int(match.group(1)) + + def _read(self, path: Path) -> xr.DataArray: + if self._subset_bounds is None or self._subset_crs is None: + da = read_raster(path, masked=False) + else: + da = read_raster_window( + path, + bounds=self._subset_bounds, + bounds_crs=self._subset_crs, + masked=False, + ) + if "band" in da.dims: + da = da.squeeze("band", drop=True) + return da + + @staticmethod + def _coords_match(a: xr.DataArray, b: xr.DataArray) -> bool: + if a.shape != b.shape: + return False + for dim in ("x", "y"): + if dim not in a.coords or dim not in b.coords: + return False + if not np.allclose(a.coords[dim].values, b.coords[dim].values): + return False + return True diff --git a/python/siac/app/_assembly_preprocessor.py b/python/siac/app/_assembly_preprocessor.py index a3b463d1..6186d15a 100644 --- a/python/siac/app/_assembly_preprocessor.py +++ b/python/siac/app/_assembly_preprocessor.py @@ -128,7 +128,7 @@ def build_preprocessor_runtime( except KeyError as exc: raise SensorNotSupportedError( f"Sensor {sensor_name!r} has no registered scene preprocessor. " - "Sentinel-2 ('s2') is the only end-to-end processing sensor currently available." + "Registered sensors: Sentinel-2 ('s2') and Landsat 8/9 OLI ('l8', 'l9')." ) from exc if input_path is not None and hasattr(preprocessor_obj, "get_metadata"): diff --git a/python/siac/catalog/sensors/landsat.py b/python/siac/catalog/sensors/landsat.py index 17755e3b..87deec8d 100644 --- a/python/siac/catalog/sensors/landsat.py +++ b/python/siac/catalog/sensors/landsat.py @@ -9,6 +9,10 @@ # specifically for Landsat. See REVIEW.md ยง3.2. _OLI_AEROSOL_BANDS = ("B1", "B2") +# Collection-2 Level-1 TOA reflectance rescaling (REFLECTANCE_MULT/ADD_BAND_n). +# The scene preprocessor reads the per-band values from the MTL; these are the +# nominal defaults. (2.75e-5 / -0.2 is the Level-2 surface-reflectance scale.) + LANDSAT8_OLI_CONFIG = SensorConfig( sensor_id="OLI", satellite_id="L8", @@ -21,8 +25,8 @@ rsrf_band("B6", 1608.5, 85.0, 30.0, 5, sensor_unit_id="landsat-8_oli"), rsrf_band("B7", 2200.5, 187.0, 30.0, 6, sensor_unit_id="landsat-8_oli"), ), - default_ref_scale=2.75e-5, - default_ref_offset=-0.2, + default_ref_scale=2.0e-5, + default_ref_offset=-0.1, aerosol_solver_band_names=_OLI_AEROSOL_BANDS, ) @@ -38,7 +42,7 @@ rsrf_band("B6", 1608.5, 85.0, 30.0, 5, sensor_unit_id="landsat-9_oli2"), rsrf_band("B7", 2200.5, 187.0, 30.0, 6, sensor_unit_id="landsat-9_oli2"), ), - default_ref_scale=2.75e-5, - default_ref_offset=-0.2, + default_ref_scale=2.0e-5, + default_ref_offset=-0.1, aerosol_solver_band_names=_OLI_AEROSOL_BANDS, ) diff --git a/python/siac/config/types.py b/python/siac/config/types.py index 27d7423c..8e7ea3d2 100644 --- a/python/siac/config/types.py +++ b/python/siac/config/types.py @@ -26,6 +26,8 @@ def __str__(self) -> str: class SensorName(SIACEnum): S2 = "s2" + L8 = "l8" + L9 = "l9" AUTO = "auto" @@ -258,6 +260,9 @@ class CloudMaskMode(SIACEnum): class CloudMaskProvider(SIACEnum): OMNICLOUDMASK = "omnicloudmask" + # Cloud/shadow classification shipped with the Level-1 product + # (Landsat Collection-2 QA_PIXEL). + PRODUCT_QA = "product_qa" class ResolutionPolicy(SIACEnum): diff --git a/tests/integration/test_e2e_synthetic_landsat.py b/tests/integration/test_e2e_synthetic_landsat.py new file mode 100644 index 00000000..8ba02b66 --- /dev/null +++ b/tests/integration/test_e2e_synthetic_landsat.py @@ -0,0 +1,259 @@ +"""Synthetic OLI M1-M6 integration and Collection-2 preprocessing seam tests. + +Like the Sentinel-style synthetic E2E test, M2/M3/M5/M6 and RT are injected; +M4 is real. No emulator forward pass, BRDF kernels, or network are required. +""" + +from __future__ import annotations + +from datetime import datetime +from pathlib import Path +from typing import TYPE_CHECKING, Any, cast + +import numpy as np +import pytest +import xarray as xr + +from siac.adapters.auth import CredentialManager +from siac.adapters.satellite.landsat import LandsatPreprocessor +from siac.adapters.satellite.observation import preprocess_to_observation_bundle +from siac.algorithms.grid.assembler import assemble_grids +from siac.app.planning import build_execution_plan +from siac.app.requests import SceneProcessRequest +from siac.catalog.sensors import get_sensor_config +from siac.config import SIACConfig +from siac.runtime import ( + AtmosphericState, + CorrectionResult, + GeometryAngles, + ObservationBundle, + SolvedAtmosphere, + SolverInputBundle, + SurfacePrior, +) +from siac.runtime.validation import ( + validate_atmospheric_state, + validate_observation_bundle, + validate_solver_input_bundle, + validate_surface_prior, +) +from siac.workflows.pipeline import run_pipeline +from siac.workflows.scene import execute_plan +from tests.integration.test_e2e_synthetic import _make_atmo, _make_surface +from tests.landsat_scene import CRS as SCENE_CRS +from tests.landsat_scene import HEIGHT, ORIGIN_X, ORIGIN_Y, WIDTH, _make_scene + +if TYPE_CHECKING: + from siac.adapters.satellite.observation import SatellitePreprocessorLike + from siac.domain import SensorConfig + +pytestmark = pytest.mark.integration +SHAPE = (64, 64) +BANDS = {f"B{n}" for n in range(1, 8)} +CRS = "EPSG:32632" +BOUNDS = (300000.0, 5500000.0, 301920.0, 5501920.0) + + +@pytest.fixture(autouse=True) +def _forbid_network(monkeypatch: pytest.MonkeyPatch) -> None: + def fail_request(*args: Any, **kwargs: Any) -> Any: + pytest.fail("Synthetic Landsat tests must not make network requests") + + monkeypatch.setattr("requests.sessions.Session.request", fail_request) + + +def _make_observation(config: SensorConfig) -> ObservationBundle: + rng = np.random.RandomState(99) + coords = { + "x": BOUNDS[0] + (np.arange(SHAPE[1]) + 0.5) * 30.0, + "y": BOUNDS[3] - (np.arange(SHAPE[0]) + 0.5) * 30.0, + } + + def raster(values: Any) -> xr.DataArray: + return cast( + "xr.DataArray", xr.DataArray(values, dims=("y", "x"), coords=coords).rio.write_crs(CRS) + ) + + toa = xr.Dataset( + { + band.name: raster(rng.uniform(0.05, 0.4, SHAPE).astype(np.float32)) + for band in config.bands + } + ) + cloud = np.zeros(SHAPE, dtype=bool) + cloud[20:40, 20:40] = True + return ObservationBundle( + toa=toa, + geometry=GeometryAngles( + sza=raster(np.full(SHAPE, 0.5)), + saa=raster(np.full(SHAPE, 2.5)), + vza=raster(np.full(SHAPE, 0.1)), + vaa=raster(np.full(SHAPE, 1.5)), + ), + cloud_mask=raster(cloud), + sensor_config=config, + metadata={"observation_time": datetime(2023, 7, 15, 10, 30)}, + crs=CRS, + bounds=BOUNDS, + ) + + +def _injected_stages(rt_model: Any, grids: list[SolverInputBundle]) -> dict[str, Any]: + atmo = _make_atmo() + surface = _make_surface() + validate_atmospheric_state(atmo) + validate_surface_prior(surface) + + def atmo_provider( + bounds: tuple[float, float, float, float], crs: str, obs_time: datetime, resolution: float + ) -> AtmosphericState: + return atmo + + def surface_provider( + observation: ObservationBundle, + atmo_state: AtmosphericState, + rt_model: Any, + resolution: float, + ) -> SurfacePrior: + return surface + + def grid_assembler( + observation: ObservationBundle, + atmo_state: AtmosphericState, + surface_prior: SurfacePrior, + rt: Any, + **kwargs: Any, + ) -> SolverInputBundle: + # The default pipeline water mask is remote; synthetic tests omit that optional input. + kwargs["water_mask_path"] = None + # PSF setup uses native BRDF kernels; omit it as in the original synthetic E2E. + kwargs["toa_psf_config"] = None + bundle = assemble_grids(observation, atmo_state, surface_prior, rt, **kwargs) + validate_solver_input_bundle(bundle) + grids.append(bundle) + return bundle + + def solver(inputs: SolverInputBundle, config: Any) -> SolvedAtmosphere: + return SolvedAtmosphere( + atmo_state=inputs.atmo_prior, + aot=inputs.atmo_prior.aot, + tcwv=inputs.atmo_prior.tcwv, + aot_unc=inputs.atmo_prior.aot_unc, + tcwv_unc=inputs.atmo_prior.tcwv_unc, + cost_final=0.001, + n_iterations=3, + converged=True, + ) + + def corrector( + observation: ObservationBundle, solved: SolvedAtmosphere, rt: Any + ) -> CorrectionResult: + return CorrectionResult( + boa=observation.toa, + boa_unc=None, + aot=solved.aot, + tcwv=solved.tcwv, + cloud_mask=observation.cloud_mask, + metadata={"corrector": "synthetic_identity"}, + ) + + return { + "atmo_provider": atmo_provider, + "surface_prior_provider": surface_provider, + "grid_assembler": grid_assembler, + "solver": solver, + "corrector": corrector, + "rt_model": rt_model, + } + + +def _assert_grid(bundle: SolverInputBundle, config: SensorConfig) -> None: + assert bundle.sensor_config == config + assert set(bundle.toa.band.values) == { + band.name for band in config.default_aerosol_solver_bands() + } + assert bundle.cloud_mask.values.any() + assert not bundle.cloud_mask.values.all() + + +@pytest.mark.parametrize("satellite", ["L8", "L9"]) +def test_full_oli_pipeline_with_injected_stages(satellite: str, mock_rt_model: Any) -> None: + config = get_sensor_config("OLI", satellite) + assert {band.name for band in config.bands} == BANDS + assert all(band.resolution == 30.0 for band in config.bands) + observation = _make_observation(config) + validate_observation_bundle(observation) + grids: list[SolverInputBundle] = [] + + def preprocessor(path: Path, aoi: Any = None) -> ObservationBundle: + return observation + + result = run_pipeline( + Path("synthetic-oli"), + None, + None, + preprocessor=preprocessor, + **_injected_stages(mock_rt_model, grids), + ) + + assert isinstance(result, CorrectionResult) + assert set(result.boa.data_vars) == BANDS + assert len(grids) == 1 + _assert_grid(grids[0], config) + assert grids[0].toa.shape[1:] == (16, 16) # 64 native 30 m pixels -> 120 m grid + assert str(grids[0].toa.rio.crs) == CRS + assert grids[0].toa.rio.bounds() == pytest.approx(BOUNDS) + assert np.all(grids[0].cloud_mask.values[5:10, 5:10]) + for name in BANDS: + assert result.boa[name].shape == SHAPE + assert str(result.boa[name].rio.crs) == observation.crs + assert result.boa[name].rio.bounds() == pytest.approx(observation.bounds) + assert result.boa[name].rio.resolution() == pytest.approx((30.0, -30.0)) + xr.testing.assert_identical(result.cloud_mask, observation.cloud_mask) + + +@pytest.mark.parametrize("spacecraft", ["LANDSAT_8", "LANDSAT_9"]) +def test_collection2_preprocessing_and_execution_plan( + tmp_path: Path, spacecraft: str, mock_rt_model: Any +) -> None: + scene = _make_scene(tmp_path, spacecraft=spacecraft) + observation = preprocess_to_observation_bundle( + # The adapter only reads sensor_config, which Landsat exposes as a property. + cast( + "SatellitePreprocessorLike", + LandsatPreprocessor(config={"cloud_mask": {"provider": "product_qa"}}), + ), + scene, + ) + validate_observation_bundle(observation) + assert set(observation.toa.data_vars) == BANDS + assert observation.crs == SCENE_CRS + assert observation.bounds == pytest.approx( + (ORIGIN_X, ORIGIN_Y - HEIGHT * 30.0, ORIGIN_X + WIDTH * 30.0, ORIGIN_Y) + ) + assert bool(observation.cloud_mask.values[12, 12]) + assert not bool(observation.cloud_mask.values[7, 0]) + + grids: list[SolverInputBundle] = [] + config = SIACConfig.model_validate( + {"sensor": "auto", "algorithms": {"cloud_mask": {"provider": "product_qa"}}} + ) + plan = build_execution_plan( + SceneProcessRequest(config=config, input_path=scene, auth=CredentialManager()), + **_injected_stages(mock_rt_model, grids), + ) + planned_observation = plan.preprocessor(scene, None) + validate_observation_bundle(planned_observation) + assert planned_observation.sensor_config.sensor_id == "OLI" + assert planned_observation.sensor_config.satellite_id == spacecraft.replace("LANDSAT_", "L") + assert planned_observation.bounds == observation.bounds + xr.testing.assert_identical(planned_observation.cloud_mask, observation.cloud_mask) + + # Planning eagerly loads cloud/reference bands; real M4 loads remaining solver bands lazily. + assert "B4" in planned_observation.toa + result = execute_plan(plan) + assert len(grids) == 1 + _assert_grid(grids[0], planned_observation.sensor_config) + assert str(result.boa.rio.crs) == SCENE_CRS + assert result.boa.rio.bounds() == pytest.approx(observation.bounds) + xr.testing.assert_identical(result.cloud_mask, observation.cloud_mask) diff --git a/tests/landsat_scene.py b/tests/landsat_scene.py new file mode 100644 index 00000000..261a8f8a --- /dev/null +++ b/tests/landsat_scene.py @@ -0,0 +1,134 @@ +"""Synthetic Landsat Collection-2 Level-1 files shared by unit and seam tests.""" + +from __future__ import annotations + +import json +from typing import TYPE_CHECKING + +import numpy as np +import rasterio +from rasterio.transform import from_origin + +from siac.adapters.satellite.landsat import parse_mtl_text + +if TYPE_CHECKING: + from pathlib import Path + +PRODUCT_ID = "LC09_L1TP_199024_20240615_20240615_02_T1" +HEIGHT, WIDTH = 40, 50 +ORIGIN_X, ORIGIN_Y = 600000.0, 5700000.0 +CRS = "EPSG:32631" +MULT, ADD = 2.0e-5, -0.1 +SZA_DEG = 35.0 + +# QA_PIXEL test values (Collection-2 bit layout). +QA_FILL = 1 << 0 +QA_CLEAR = 1 << 6 +QA_CLOUD = (1 << 3) | (1 << 1) +QA_SHADOW = 1 << 4 +QA_CIRRUS = 1 << 2 + + +def _write(path: Path, data: np.ndarray, dtype: str) -> None: + with rasterio.open( + path, + "w", + driver="GTiff", + height=data.shape[0], + width=data.shape[1], + count=1, + dtype=dtype, + crs=CRS, + transform=from_origin(ORIGIN_X, ORIGIN_Y, 30.0, 30.0), + ) as dst: + dst.write(data.astype(dtype), 1) + + +def _mtl_text(spacecraft: str, level: str, *, with_angle_names: bool) -> str: + files = "\n".join(f' FILE_NAME_BAND_{n} = "{PRODUCT_ID}_B{n}.TIF"' for n in range(1, 8)) + angles = "" + if with_angle_names: + angles = ( + f' FILE_NAME_ANGLE_SOLAR_ZENITH_BAND_4 = "{PRODUCT_ID}_SZA.TIF"\n' + f' FILE_NAME_ANGLE_SOLAR_AZIMUTH_BAND_4 = "{PRODUCT_ID}_SAA.TIF"\n' + f' FILE_NAME_ANGLE_SENSOR_ZENITH_BAND_4 = "{PRODUCT_ID}_VZA.TIF"\n' + f' FILE_NAME_ANGLE_SENSOR_AZIMUTH_BAND_4 = "{PRODUCT_ID}_VAA.TIF"\n' + f' FILE_NAME_QUALITY_L1_PIXEL = "{PRODUCT_ID}_QA_PIXEL.TIF"\n' + ) + rescale = "\n".join( + f" REFLECTANCE_MULT_BAND_{n} = {MULT:.4E}\n REFLECTANCE_ADD_BAND_{n} = {ADD:.6f}" + for n in range(1, 10) + ) + return f"""GROUP = LANDSAT_METADATA_FILE + GROUP = PRODUCT_CONTENTS + LANDSAT_PRODUCT_ID = "{PRODUCT_ID}" + COLLECTION_NUMBER = 02 + COLLECTION_CATEGORY = "T1" + PROCESSING_LEVEL = "{level}" +{files} +{angles} END_GROUP = PRODUCT_CONTENTS + GROUP = IMAGE_ATTRIBUTES + SPACECRAFT_ID = "{spacecraft}" + SENSOR_ID = "OLI_TIRS" + WRS_PATH = 199 + WRS_ROW = 24 + DATE_ACQUIRED = 2024-06-15 + SCENE_CENTER_TIME = "10:47:12.3456789Z" + CLOUD_COVER = 12.5 + SUN_AZIMUTH = 150.0 + SUN_ELEVATION = {90.0 - SZA_DEG} + EARTH_SUN_DISTANCE = 1.0158 + END_GROUP = IMAGE_ATTRIBUTES + GROUP = LEVEL1_RADIOMETRIC_RESCALING +{rescale} + END_GROUP = LEVEL1_RADIOMETRIC_RESCALING +END_GROUP = LANDSAT_METADATA_FILE +END +""" + + +def _qa_grid() -> np.ndarray: + qa = np.full((HEIGHT, WIDTH), QA_CLEAR, dtype=np.uint16) + qa[:5, :] = QA_FILL + qa[10:15, 10:20] = QA_CLOUD + qa[20:25, 10:20] = QA_SHADOW + qa[30:35, 30:40] = QA_CIRRUS + return qa + + +def _make_scene( + root: Path, + *, + spacecraft: str = "LANDSAT_9", + level: str = "L1TP", + with_angles: bool = True, + with_angle_names: bool = True, + mtl_format: str = "txt", +) -> Path: + scene = root / PRODUCT_ID + scene.mkdir() + fill_rows = slice(0, 5) + for n in range(1, 8): + dn = np.full((HEIGHT, WIDTH), 10000 + 1000 * n, dtype=np.uint16) + dn[fill_rows, :] = 0 + _write(scene / f"{PRODUCT_ID}_B{n}.TIF", dn, "uint16") + if with_angles: + sza = np.full((HEIGHT, WIDTH), SZA_DEG * 100) + saa = np.full((HEIGHT, WIDTH), 150.0 * 100) + vza = np.tile(np.linspace(-700, 700, WIDTH), (HEIGHT, 1)) + vza = np.abs(vza) + # Azimuths are stored in [-180, 180] deg; alternate -179 / +179 deg inside each block. + vaa = np.where(np.arange(WIDTH)[None, :] % 2 == 0, -17900, 17900) * np.ones((HEIGHT, 1)) + for name, grid in (("SZA", sza), ("SAA", saa), ("VZA", vza), ("VAA", vaa)): + grid = grid.copy() + grid[fill_rows, :] = 0 + _write(scene / f"{PRODUCT_ID}_{name}.TIF", grid, "int16") + _write(scene / f"{PRODUCT_ID}_QA_PIXEL.TIF", _qa_grid(), "uint16") + + text = _mtl_text(spacecraft, level, with_angle_names=with_angle_names and with_angles) + if mtl_format == "json": + payload = json.dumps(parse_mtl_text(text)) + (scene / f"{PRODUCT_ID}_MTL.json").write_text(payload, encoding="utf-8") + else: + (scene / f"{PRODUCT_ID}_MTL.txt").write_text(text, encoding="utf-8") + return scene diff --git a/tests/unit/test_config.py b/tests/unit/test_config.py index 8639f97d..269098f6 100644 --- a/tests/unit/test_config.py +++ b/tests/unit/test_config.py @@ -243,7 +243,11 @@ def test_with_overrides(self): def test_unsupported_processing_sensor_is_rejected(self): with pytest.raises(ValidationError): - SIACConfig(sensor="l8") + SIACConfig(sensor="l5") + + @pytest.mark.parametrize("sensor", ["l8", "l9"]) + def test_landsat_processing_sensors_are_accepted(self, sensor): + assert SIACConfig(sensor=sensor).sensor.value == sensor def test_resolve_uses_python_only_sensor_and_aoi_defaults(self): config = SIACConfig(sensor="s2", aoi="/tmp/aoi.geojson") diff --git a/tests/unit/test_landsat_preprocessor.py b/tests/unit/test_landsat_preprocessor.py new file mode 100644 index 00000000..718241fc --- /dev/null +++ b/tests/unit/test_landsat_preprocessor.py @@ -0,0 +1,249 @@ +"""Unit tests for the Landsat 8/9 OLI Collection-2 Level-1 preprocessor.""" + +from __future__ import annotations + +from datetime import datetime +from typing import TYPE_CHECKING + +import numpy as np +import pytest + +from siac.adapters.satellite import detect_sensor, get_preprocessor +from siac.adapters.satellite.landsat import ( + LandsatPreprocessor, + flatten_mtl, + parse_mtl_text, + qa_pixel_to_cloud_classes, +) +from siac.adapters.satellite.observation import preprocess_to_observation_bundle +from siac.algorithms.rt.emulator.two_nn import TwoLayerNNEmulator +from siac.runtime.validation import validate_observation_bundle +from tests.landsat_scene import ( + ADD, + CRS, + HEIGHT, + MULT, + ORIGIN_X, + PRODUCT_ID, + QA_CIRRUS, + QA_CLEAR, + QA_CLOUD, + QA_FILL, + QA_SHADOW, + SZA_DEG, + WIDTH, + _make_scene, + _mtl_text, + _qa_grid, +) + +if TYPE_CHECKING: + from pathlib import Path + + +def _expected_reflectance(n: int, sza_deg: float = SZA_DEG) -> float: + return ((10000 + 1000 * n) * MULT + ADD) / np.cos(np.deg2rad(sza_deg)) + + +# --------------------------------------------------------------------------- +# MTL parsing +# --------------------------------------------------------------------------- + + +def test_parse_mtl_text_nests_groups_and_strips_quotes() -> None: + tree = parse_mtl_text(_mtl_text("LANDSAT_8", "L1TP", with_angle_names=False)) + contents = tree["LANDSAT_METADATA_FILE"]["PRODUCT_CONTENTS"] + assert contents["LANDSAT_PRODUCT_ID"] == PRODUCT_ID + assert contents["COLLECTION_NUMBER"] == "02" + flat = flatten_mtl(tree) + assert flat["SPACECRAFT_ID"] == "LANDSAT_8" + assert float(flat["REFLECTANCE_MULT_BAND_4"]) == pytest.approx(MULT) + + +# --------------------------------------------------------------------------- +# Registration and detection +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize(("spacecraft", "expected"), [("LANDSAT_8", "l8"), ("LANDSAT_9", "l9")]) +def test_detect_sensor_distinguishes_landsat_8_and_9( + tmp_path: Path, spacecraft: str, expected: str +) -> None: + scene = _make_scene(tmp_path, spacecraft=spacecraft) + assert detect_sensor(scene) == expected + assert detect_sensor(scene / f"{PRODUCT_ID}_MTL.txt") == expected + + +def test_detect_sensor_reads_json_only_scene(tmp_path: Path) -> None: + scene = _make_scene(tmp_path, mtl_format="json") + assert detect_sensor(scene) == "l9" + + +@pytest.mark.parametrize("sensor_id", ["l8", "l9", "landsat"]) +def test_registry_resolves_landsat_preprocessor(sensor_id: str) -> None: + assert isinstance(get_preprocessor(sensor_id), LandsatPreprocessor) + + +def test_platform_comes_from_mtl_not_registry_key(tmp_path: Path) -> None: + scene = _make_scene(tmp_path, spacecraft="LANDSAT_9") + pre = get_preprocessor("l8") + pre.get_metadata(scene) + assert pre.sensor_config.sensor_id == "OLI" + assert pre.sensor_config.satellite_id == "L9" + + +@pytest.mark.parametrize( + ("spacecraft", "level", "message"), + [("LANDSAT_7", "L1TP", "Unsupported Landsat spacecraft"), ("LANDSAT_8", "L2SP", "Level-1")], +) +def test_rejects_unsupported_scenes( + tmp_path: Path, spacecraft: str, level: str, message: str +) -> None: + scene = _make_scene(tmp_path, spacecraft=spacecraft, level=level) + with pytest.raises(ValueError, match=message): + LandsatPreprocessor().get_metadata(scene) + + +# --------------------------------------------------------------------------- +# Metadata +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("mtl_format", ["txt", "json"]) +def test_metadata_from_mtl(tmp_path: Path, mtl_format: str) -> None: + scene = _make_scene(tmp_path, mtl_format=mtl_format) + metadata = LandsatPreprocessor().get_metadata(scene) + assert metadata["observation_time"] == datetime(2024, 6, 15, 10, 47, 12, 345678) + assert metadata["satellite"] == "L9" + assert metadata["tile_id"] == "199024" + assert metadata["sun_elevation"] == pytest.approx(90.0 - SZA_DEG) + assert metadata["reflectance_mult"]["B7"] == pytest.approx(MULT) + assert set(metadata["reflectance_add"]) == {f"B{n}" for n in range(1, 8)} + + +# --------------------------------------------------------------------------- +# TOA reflectance +# --------------------------------------------------------------------------- + + +def test_load_toa_applies_rescaling_and_per_pixel_solar_zenith(tmp_path: Path) -> None: + scene = _make_scene(tmp_path) + toa = LandsatPreprocessor().load_toa(scene) + assert set(toa.data_vars) == {f"B{n}" for n in range(1, 8)} + for n in range(1, 8): + band = toa[f"B{n}"] + assert band.dtype == np.float32 + assert np.all(np.isnan(band.values[:5])) + np.testing.assert_allclose(band.values[5:], _expected_reflectance(n), rtol=1e-5) + assert toa.attrs["sensor"] == "OLI" + assert toa.attrs["satellite"] == "L9" + + +def test_load_toa_falls_back_to_scene_sun_elevation(tmp_path: Path) -> None: + scene = _make_scene(tmp_path, with_angles=False) + toa = LandsatPreprocessor().load_toa(scene, band_names=("B4",)) + assert list(toa.data_vars) == ["B4"] + np.testing.assert_allclose(toa["B4"].values[5:], _expected_reflectance(4), rtol=1e-5) + + +def test_band_loader_serves_bands_not_preloaded(tmp_path: Path) -> None: + scene = _make_scene(tmp_path) + toa = LandsatPreprocessor().load_toa(scene, band_names=("B4",)) + loader = toa.attrs["_siac_toa_band_loader"] + b6 = loader("B6") + np.testing.assert_allclose(b6.values[5:], _expected_reflectance(6), rtol=1e-5) + with pytest.raises(KeyError): + loader("B4") + + +def test_band_files_found_by_suffix_when_mtl_names_absent(tmp_path: Path) -> None: + scene = _make_scene(tmp_path, with_angle_names=False) + geometry = LandsatPreprocessor().extract_geometry(scene) + assert float(geometry.sza.mean()) == pytest.approx(np.deg2rad(SZA_DEG), rel=1e-5) + + +# --------------------------------------------------------------------------- +# Geometry +# --------------------------------------------------------------------------- + + +def test_extract_geometry_block_averages_angle_bands(tmp_path: Path) -> None: + scene = _make_scene(tmp_path) + pre = LandsatPreprocessor(config={"angle_block_pixels": 10}) + geometry = pre.extract_geometry(scene) + assert geometry.sza.shape == (4, 5) + # Fill rows are excluded from the block mean, and no block is left empty. + np.testing.assert_allclose(geometry.sza.values, np.deg2rad(SZA_DEG), rtol=1e-5) + np.testing.assert_allclose(geometry.saa.values, np.deg2rad(150.0), rtol=1e-5) + assert np.all(np.isfinite(geometry.vza.values)) + # -179 deg and +179 deg average circularly to 180 deg (an arithmetic mean gives 0). + vaa_deg = np.mod(np.rad2deg(geometry.vaa.values), 360.0) + np.testing.assert_allclose(vaa_deg, 180.0, atol=1e-3) + # Block centres sit on the native grid. + assert float(geometry.sza.x[0]) == pytest.approx(ORIGIN_X + 15.0 + 4.5 * 30.0) + assert str(geometry.sza.rio.crs) == CRS + + +def test_extract_geometry_requires_angle_bands(tmp_path: Path) -> None: + scene = _make_scene(tmp_path, with_angles=False) + with pytest.raises(FileNotFoundError, match="angle band"): + LandsatPreprocessor().extract_geometry(scene) + + +# --------------------------------------------------------------------------- +# Cloud classes +# --------------------------------------------------------------------------- + + +def test_qa_pixel_bits_map_to_standard_classes() -> None: + qa = np.array([QA_FILL, QA_CLEAR, QA_CLOUD, QA_SHADOW, QA_CIRRUS, QA_CLOUD | QA_SHADOW]) + np.testing.assert_array_equal(qa_pixel_to_cloud_classes(qa), [0, 1, 2, 3, 2, 2]) + + +def test_product_qa_cloud_mask(tmp_path: Path) -> None: + scene = _make_scene(tmp_path) + pre = LandsatPreprocessor(config={"cloud_mask": {"provider": "product_qa"}}) + mask = pre.extract_cloud_mask(scene) + classes = pre._last_cloud_classes + assert classes is not None + np.testing.assert_array_equal(classes.values, qa_pixel_to_cloud_classes(_qa_grid())) + assert mask.dtype == bool + assert bool(mask.values[0, 0]) # fill is not usable + assert not bool(mask.values[7, 0]) # clear + assert bool(mask.values[12, 12]) # cloud + assert bool(mask.values[22, 12]) # shadow + + +# --------------------------------------------------------------------------- +# M1 contract +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + "cloud_mask", [{"provider": "product_qa"}, {"mode": "none"}], ids=["product_qa", "none"] +) +def test_preprocess_produces_valid_observation_bundle(tmp_path: Path, cloud_mask: dict) -> None: + scene = _make_scene(tmp_path) + pre = LandsatPreprocessor(config={"cloud_mask": cloud_mask, "angle_block_pixels": 10}) + bundle = preprocess_to_observation_bundle(pre, scene) + validate_observation_bundle(bundle) + assert bundle.sensor_config.satellite_id == "L9" + assert bundle.crs == CRS + assert bundle.cloud_mask.shape == (HEIGHT, WIDTH) + assert set(bundle.toa.data_vars) == {f"B{n}" for n in range(1, 8)} + + +# --------------------------------------------------------------------------- +# RT emulator platform fallback +# --------------------------------------------------------------------------- + + +def test_landsat9_emulator_falls_back_to_landsat8_weights(tmp_path: Path) -> None: + for n in range(1, 8): + for part in ("xap", "xbp", "xcp"): + (tmp_path / f"isotropic_OLI_L8_B{n}_{part}.npz").touch() + emulator = TwoLayerNNEmulator(emulator_dir=tmp_path, sensor_id="OLI", satellite_id="L9") + assert emulator.available_bands == [f"B{n}" for n in range(1, 8)] + path = emulator._get_emulator_path("B4") + assert path is not None + assert "_L8_B4_" in path.name