Add reprojection and registration

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.

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`.

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.

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. 

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

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(

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(