diff --git a/README.md b/README.md index e96deb7..a243ae3 100644 --- a/README.md +++ b/README.md @@ -22,4 +22,4 @@ Full documentation, including the core data model, processing registries, and cu ## Contributing -See [CONTRIBUTING.md](CONTRIBUTING.md) for how to set up a development environment and submit changes. +See [CONTRIBUTING.md](CONTRIBUTING.md) for how to set up a development environment and submit changes. \ No newline at end of file diff --git a/docs/concepts.md b/docs/concepts.md index 6763d2c..5b4e0b9 100644 --- a/docs/concepts.md +++ b/docs/concepts.md @@ -35,3 +35,15 @@ regardless of which field an importer populated, `scene.cloud_mask`, `scene.clou `scene.masks["cloud"]`, and `scene.metadata["cloud_mask"]` all agree. Downstream algorithms (cloud masking, cloud filtering) read from this normalised schema rather than product-specific field names. + +## Spectral indices + +`earthrs.indices` provides `ndvi`, `ndwi`, and `evi`, each taking a `Scene` and returning +per-pixel index values computed from its named bands. + +## Sampling + +`earthrs.sampling` provides `sample_points`, `sample_polygons`, and `sample_transects` (plus a +`sample()` dispatcher that infers the target from geometry type), all returning `Samples` — one +row per point, polygon, or transect vertex observation. `Dataset.sample_points` delegates to +`earthrs.sampling.sample_points` per scene and aggregates the results into a single `Samples`. diff --git a/docs/index.md b/docs/index.md index 1ffb9f7..a6e5bb5 100644 --- a/docs/index.md +++ b/docs/index.md @@ -16,9 +16,9 @@ GeoPandas) instead of replacing it. ## Where to go next - [Getting started](getting-started.md) — install the package and run a first example. -- [Core concepts](concepts.md) — the `Scene`, `Dataset`, and `Samples` data model, and cloud - metadata normalisation. -- [Processing](processing.md) — glint removal, depth correction, and cloud masking, and how to - register new methods. +- [Core concepts](concepts.md) — the `Scene`, `Dataset`, and `Samples` data model, cloud + metadata normalisation, spectral indices, and sampling. +- [Processing](processing.md) — glint removal, depth correction, cloud masking, + reprojection/registration, and how to register new methods. - [Project status](status.md) — what's implemented, what's known-limited, and what's scaffolded but not yet built. diff --git a/docs/processing.md b/docs/processing.md index 9fa285a..f230671 100644 --- a/docs/processing.md +++ b/docs/processing.md @@ -118,3 +118,16 @@ depth_correct( None of the three variants take a known `depth` as input — the whole point is to derive depth (or a depth-invariant index) from radiance, unlike `"maritorena"` which requires one. + +## Reprojection and registration + +`earthrs.processing.reproject` resamples a scene onto a new pixel grid — described by an +`(a, b, c, d, e, f)` affine `transform` and a `(rows, cols)` `shape` — using nearest-neighbour or +bilinear resampling. `earthrs.processing.register` resamples a scene onto a reference scene's +grid (transform, shape, and CRS). Both operate on 2D grid band data (a list of rows), matching +the layout `earthrs.sampling` already assumes — see [Project status](status.md) for how this +differs from the flat/nested-list layout accepted elsewhere in `earthrs.processing`. + +Neither function bundles a CRS database: if the destination `crs` differs from `scene.crs`, you +must supply a `transformer` callable mapping `(x, y)` in `scene.crs` to `(x, y)` in the +destination CRS (for example backed by `pyproj`), otherwise a `ValueError` is raised. diff --git a/docs/status.md b/docs/status.md index 20a234c..bc8e82c 100644 --- a/docs/status.md +++ b/docs/status.md @@ -16,4 +16,4 @@ scaffolded but not yet built — see [Getting started](getting-started.md) and - sensor definitions (`earthrs.sensors`) beyond placeholder name-only classes - machine-learning helpers (`earthrs.ml`) -- `register`/`reproject` raster-registration and reprojection helpers +- visualisation helpers (`earthrs.plotting`) — currently an empty package diff --git a/src/earthrs/processing/core.py b/src/earthrs/processing/core.py index e9d1668..7629a1d 100644 --- a/src/earthrs/processing/core.py +++ b/src/earthrs/processing/core.py @@ -2,6 +2,8 @@ from __future__ import annotations +import math +from collections.abc import Callable from typing import Any from earthrs.scene import Scene @@ -27,26 +29,218 @@ def register_cloud_method(name: str, func: CloudProcessor) -> None: _CLOUD_REGISTRY[name.lower()] = func -def register(data: Any, reference: Any, *, method: str = "default") -> Any: - """Register raster data to a reference scene. +_RESAMPLING_METHODS = frozenset({"nearest", "bilinear"}) - Placeholder implementation keeps API stable until full registration support is added. +Transform = tuple[float, float, float, float, float, float] +CoordTransformer = Callable[[float, float], tuple[float, float]] + + +def reproject( + scene: Scene, + *, + transform: Transform, + shape: tuple[int, int], + crs: Any = None, + resampling: str = "nearest", + transformer: CoordTransformer | None = None, +) -> Scene: + """Resample a scene onto a new pixel grid, optionally changing CRS. + + ``transform`` and ``shape`` describe the destination grid: ``transform`` is an + ``(a, b, c, d, e, f)`` affine coefficient tuple (``x = a*col + b*row + c``, + ``y = d*col + e*row + f``), and ``shape`` is ``(rows, cols)``. Each band in + ``scene.data`` must be a 2D grid (a list of rows). + + When ``crs`` differs from ``scene.crs``, a ``transformer`` callable mapping + ``(x, y)`` in ``scene.crs`` to ``(x, y)`` in ``crs`` must be supplied (for + example backed by ``pyproj``). """ - _ = method - _ = reference - return data + if scene.transform is None: + raise ValueError("`scene.transform` is required to reproject.") + target_crs = scene.crs if crs is None else crs + if transformer is None and scene.crs is not None and target_crs != scene.crs: + raise ValueError( + "Reprojecting between different CRS requires a `transformer` callable " + "(for example one backed by pyproj); earthrs does not bundle a CRS database." + ) + return _resample_scene( + scene, + dst_transform=transform, + dst_shape=shape, + crs=target_crs, + resampling=resampling, + transformer=transformer, + history_entry=f"reproject:{resampling}", + ) -def reproject(data: Any, *, crs: Any, resampling: str = "nearest") -> Any: - """Reproject raster data to a target coordinate reference system. +def register(scene: Scene, reference: Scene, *, method: str = "nearest") -> Scene: + """Align ``scene`` onto ``reference``'s pixel grid. - Placeholder implementation keeps API stable until full reprojection support is added. + Resamples ``scene`` so its transform, shape, and CRS match ``reference``. CRS + reconciliation follows the same rules as :func:`reproject`: if ``scene.crs`` + and ``reference.crs`` differ, reproject with an explicit ``transformer`` first. """ - _ = crs - _ = resampling - return data + if reference.transform is None: + raise ValueError("`reference.transform` is required to register a scene.") + if scene.crs is not None and reference.crs is not None and scene.crs != reference.crs: + raise ValueError( + "`register` cannot reconcile differing CRS without a transformer; " + "call `reproject` with an explicit `transformer` first." + ) + return _resample_scene( + scene, + dst_transform=reference.transform, + dst_shape=_scene_shape(reference), + crs=reference.crs, + resampling=method, + transformer=None, + history_entry=f"register:{method}", + ) + + +def _resample_scene( + scene: Scene, + *, + dst_transform: Transform, + dst_shape: tuple[int, int], + crs: Any, + resampling: str, + transformer: CoordTransformer | None, + history_entry: str, +) -> Scene: + if resampling not in _RESAMPLING_METHODS: + raise ValueError(f"Unsupported resampling method '{resampling}'.") + if not isinstance(scene.data, dict): + raise TypeError("Reprojection expects mapping-based scene data.") + src_transform = scene.transform + if src_transform is None: + raise ValueError("`scene.transform` is required to resample a scene.") + + updated_data = { + band: _resample_grid(grid, src_transform, dst_transform, dst_shape, resampling, transformer) + for band, grid in scene.data.items() + } + + updated_mask = None + if scene.cloud_mask is not None: + updated_mask = _resample_grid( + scene.cloud_mask, src_transform, dst_transform, dst_shape, "nearest", transformer + ) + updated_probability = None + if scene.cloud_probability is not None: + updated_probability = _resample_grid( + scene.cloud_probability, + src_transform, + dst_transform, + dst_shape, + resampling, + transformer, + ) + + masks = dict(scene.masks) + metadata = dict(scene.metadata) + if updated_mask is not None: + masks["cloud"] = updated_mask + metadata["cloud_mask"] = updated_mask + if updated_probability is not None: + masks["cloud_probability"] = updated_probability + metadata["cloud_probability"] = updated_probability + + return Scene( + data=updated_data, + crs=crs, + transform=dst_transform, + metadata=metadata, + band_names=scene.band_names, + acquisition_time=scene.acquisition_time, + sensor=scene.sensor, + history=(*scene.history, history_entry), + masks=masks, + cloud_mask=updated_mask if updated_mask is not None else scene.cloud_mask, + cloud_probability=( + updated_probability if updated_probability is not None else scene.cloud_probability + ), + ) + + +def _scene_shape(scene: Scene) -> tuple[int, int]: + if not isinstance(scene.data, dict) or not scene.data: + raise ValueError("Cannot infer grid shape from a scene without band data.") + first_band = next(iter(scene.data.values())) + if not isinstance(first_band, list) or not first_band or not isinstance(first_band[0], list): + raise ValueError("Reprojection expects each band as a 2D grid (a list of rows).") + return len(first_band), len(first_band[0]) + + +def _resample_grid( + grid: Any, + src_transform: Transform, + dst_transform: Transform, + dst_shape: tuple[int, int], + resampling: str, + transformer: CoordTransformer | None, +) -> list[list[Any]]: + if not isinstance(grid, list) or not grid or not isinstance(grid[0], list): + raise TypeError("Reprojection expects each band as a 2D grid (a list of rows).") + src_rows = len(grid) + src_cols = len(grid[0]) + dst_rows, dst_cols = dst_shape + + result: list[list[Any]] = [] + for dst_row in range(dst_rows): + row_values: list[Any] = [] + for dst_col in range(dst_cols): + x, y = _pixel_to_world(dst_col, dst_row, dst_transform) + if transformer is not None: + x, y = transformer(x, y) + src_col, src_row = _world_to_pixel(x, y, src_transform) + if resampling == "nearest": + value = _sample_nearest(grid, src_row, src_col, src_rows, src_cols) + else: + value = _sample_bilinear(grid, src_row, src_col, src_rows, src_cols) + row_values.append(value) + result.append(row_values) + return result + + +def _pixel_to_world(col: float, row: float, transform: Transform) -> tuple[float, float]: + a, b, c, d, e, f = transform + return a * col + b * row + c, d * col + e * row + f + + +def _world_to_pixel(x: float, y: float, transform: Transform) -> tuple[float, float]: + a, b, c, d, e, f = transform + det = a * e - b * d + if det == 0: + raise ValueError("Transform is not invertible.") + col = (e * (x - c) - b * (y - f)) / det + row = (-d * (x - c) + a * (y - f)) / det + return col, row + + +def _sample_nearest(grid: list[list[Any]], row: float, col: float, rows: int, cols: int) -> Any: + nearest_row = round(row) + nearest_col = round(col) + if 0 <= nearest_row < rows and 0 <= nearest_col < cols: + return grid[nearest_row][nearest_col] + return None + + +def _sample_bilinear(grid: list[list[Any]], row: float, col: float, rows: int, cols: int) -> Any: + row0 = math.floor(row) + col0 = math.floor(col) + row1 = row0 + 1 + col1 = col0 + 1 + if row0 < 0 or col0 < 0 or row1 >= rows or col1 >= cols: + return None + frac_row = row - row0 + frac_col = col - col0 + top = grid[row0][col0] * (1 - frac_col) + grid[row0][col1] * frac_col + bottom = grid[row1][col0] * (1 - frac_col) + grid[row1][col1] * frac_col + return top * (1 - frac_row) + bottom * frac_row def _resolve_cloud_method( diff --git a/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py index 3f0e4d2..8785057 100644 --- a/src/earthrs/processing/depth.py +++ b/src/earthrs/processing/depth.py @@ -43,9 +43,7 @@ def _require_mapping_data(scene: Scene, label: str) -> None: raise TypeError(f"{label} correction expects mapping-based scene data.") -def _band_param( - bands: Iterable[str], value: Any, default: float -) -> dict[str, float]: +def _band_param(bands: Iterable[str], value: Any, default: float) -> dict[str, float]: """Normalise a per-band parameter to a ``{band: value}`` mapping. Accepts ``None`` (every band gets `default`), a single scalar (every band gets @@ -112,9 +110,7 @@ def _lyzenga_1978( updated[band] = values continue updated[band] = _log_transform(values, dwr[band], epsilon) - return _finalise_lyzenga( - scene, data=updated, band_names=scene.band_names, method="lyzenga1978" - ) + return _finalise_lyzenga(scene, data=updated, band_names=scene.band_names, method="lyzenga1978") def _lyzenga_1981( diff --git a/tests/test_processing.py b/tests/test_processing.py index d84550c..11fb4c3 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -2,7 +2,14 @@ import pytest -from earthrs.processing import atmospheric_correction, cloud_mask, depth_correct, remove_glint +from earthrs.processing import ( + atmospheric_correction, + cloud_mask, + depth_correct, + register, + remove_glint, + reproject, +) from earthrs.scene import Scene @@ -288,3 +295,141 @@ def test_cloud_mask_explicit_method_overrides_sensor_auto_resolution() -> None: result = cloud_mask(scene, method="probability", probability=[0.9]) assert result.cloud_mask == [True] + + +def test_reproject_requires_scene_transform() -> None: + scene = Scene(data={"nir": [[1, 2], [3, 4]]}, band_names=["nir"]) + + with pytest.raises(ValueError): + reproject(scene, transform=(1, 0, 0, 0, 1, 0), shape=(2, 2)) + + +def test_reproject_identity_transform_reproduces_grid() -> None: + scene = Scene( + data={"nir": [[1, 2], [3, 4]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + crs="EPSG:4326", + ) + + result = reproject(scene, transform=(1, 0, 0, 0, 1, 0), shape=(2, 2)) + + assert result.data["nir"] == [[1, 2], [3, 4]] + assert result.transform == (1, 0, 0, 0, 1, 0) + assert result.history == ("reproject:nearest",) + + +def test_reproject_nearest_downsamples_by_picking_source_pixels() -> None: + grid = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]] + scene = Scene(data={"nir": grid}, band_names=["nir"], transform=(1, 0, 0, 0, 1, 0)) + + result = reproject(scene, transform=(2, 0, 0, 0, 2, 0), shape=(2, 2)) + + assert result.data["nir"] == [[1, 3], [9, 11]] + + +def test_reproject_nearest_marks_out_of_bounds_as_none() -> None: + scene = Scene(data={"nir": [[1, 2], [3, 4]]}, band_names=["nir"], transform=(1, 0, 0, 0, 1, 0)) + + result = reproject(scene, transform=(1, 0, 10, 0, 1, 10), shape=(2, 2)) + + assert result.data["nir"] == [[None, None], [None, None]] + + +def test_reproject_bilinear_interpolates_between_source_pixels() -> None: + scene = Scene( + data={"nir": [[0, 10], [20, 30]]}, band_names=["nir"], transform=(1, 0, 0, 0, 1, 0) + ) + + result = reproject(scene, transform=(1, 0, 0.5, 0, 1, 0.5), shape=(1, 1), resampling="bilinear") + + assert result.data["nir"][0] == pytest.approx([15.0]) + + +def test_reproject_rejects_unsupported_resampling_method() -> None: + scene = Scene(data={"nir": [[1, 2], [3, 4]]}, band_names=["nir"], transform=(1, 0, 0, 0, 1, 0)) + + with pytest.raises(ValueError): + reproject(scene, transform=(1, 0, 0, 0, 1, 0), shape=(2, 2), resampling="cubic") + + +def test_reproject_cross_crs_requires_transformer() -> None: + scene = Scene( + data={"nir": [[1, 2], [3, 4]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + crs="EPSG:4326", + ) + + with pytest.raises(ValueError): + reproject(scene, transform=(1, 0, 0, 0, 1, 0), shape=(2, 2), crs="EPSG:3857") + + +def test_reproject_cross_crs_uses_transformer() -> None: + scene = Scene( + data={"nir": [[1, 2], [3, 4]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + crs="EPSG:4326", + ) + + result = reproject( + scene, + transform=(1, 0, 0, 0, 1, 0), + shape=(2, 2), + crs="EPSG:3857", + transformer=lambda x, y: (x, y), + ) + + assert result.data["nir"] == [[1, 2], [3, 4]] + assert result.crs == "EPSG:3857" + + +def test_reproject_resamples_cloud_mask_with_nearest_regardless_of_method() -> None: + scene = Scene( + data={"nir": [[1.0, 2.0], [3.0, 4.0]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + masks={"cloud": [[True, False], [False, True]]}, + ) + + result = reproject(scene, transform=(1, 0, 0.5, 0, 1, 0.5), shape=(1, 1), resampling="bilinear") + + assert result.cloud_mask == [[True]] + assert result.masks["cloud"] == result.cloud_mask + + +def test_register_aligns_scene_onto_reference_grid() -> None: + scene = Scene( + data={"nir": [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], [13, 14, 15, 16]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + ) + reference = Scene(data={"nir": [[0, 0], [0, 0]]}, transform=(2, 0, 0, 0, 2, 0)) + + result = register(scene, reference) + + assert result.data["nir"] == [[1, 3], [9, 11]] + assert result.transform == (2, 0, 0, 0, 2, 0) + assert result.history == ("register:nearest",) + + +def test_register_requires_reference_transform() -> None: + scene = Scene(data={"nir": [[1, 2], [3, 4]]}, band_names=["nir"], transform=(1, 0, 0, 0, 1, 0)) + reference = Scene(data={"nir": [[0, 0], [0, 0]]}) + + with pytest.raises(ValueError): + register(scene, reference) + + +def test_register_rejects_mismatched_crs_without_transformer() -> None: + scene = Scene( + data={"nir": [[1, 2], [3, 4]]}, + band_names=["nir"], + transform=(1, 0, 0, 0, 1, 0), + crs="EPSG:4326", + ) + reference = Scene(data={"nir": [[0, 0], [0, 0]]}, transform=(1, 0, 0, 0, 1, 0), crs="EPSG:3857") + + with pytest.raises(ValueError): + register(scene, reference)