Georeferencing (GCPs & RPCs)

Hold "Ctrl" to enable pan & zoom

flowchart LR
    GE(("Georef<br/>ds.georef"))
    GE --> GC["<b>ground-control points</b><br/>gcps · gcp_count · gcp_projection<br/>has_gcps · set_gcps"]
    GE --> RP["<b>rational-polynomial coeffs</b><br/>rpcs · has_rpcs · set_rpcs"]
    GE --> WA["<b>warp to a map grid</b><br/>orthorectify · georeference"]

Read, attach, and warp-from ground-control points and rational-polynomial coefficients to georeference raw / un-orthorectified imagery — the one case where a raster is not described by a simple affine geotransform. Accessed as ds.georef, with same-named facades on Dataset (ds.gcps, ds.set_gcps, ds.georeference, ds.rpcs, ds.set_rpcs, ds.orthorectify).

Everything routes through GDAL — pyramids does not implement any sensor model itself.

pyramids.dataset.engines.Georef

Bases: _Engine['Dataset']

Ground-control-point and RPC georeferencing for a :class:Dataset.

A normal raster is georeferenced by an affine geotransform. Raw imagery (scanned maps, drone mosaics, un-orthorectified satellite scenes) instead carries ground-control points (pixel↔map tie points) or rational polynomial coefficients (a vendor sensor model). This engine reads and attaches both, and warps from them into an affine-geotransform raster.

Source code in src/pyramids/dataset/engines/georef.py

class Georef(_Engine["Dataset"]):
    """Ground-control-point and RPC georeferencing for a :class:`Dataset`.

    A normal raster is georeferenced by an affine geotransform. Raw imagery
    (scanned maps, drone mosaics, un-orthorectified satellite scenes) instead
    carries **ground-control points** (pixel↔map tie points) or **rational
    polynomial coefficients** (a vendor sensor model). This engine reads and
    attaches both, and warps from them into an affine-geotransform raster.
    """

    @property
    def gcps(self: Georef) -> list[GroundControlPoint]:
        """The dataset's ground-control points (empty list when it has none).

        Returns:
            list[GroundControlPoint]: one per attached GCP, in GDAL order.

        Examples:
            - Read back the points attached with :meth:`set_gcps`:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from pyramids.dataset._gcp import GroundControlPoint
                >>> ds = Dataset.create_from_array(
                ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
                ... )
                >>> ds.set_gcps([GroundControlPoint(row=0, col=0, x=10.0, y=50.0)], 4326)
                >>> (ds.gcps[0].col, ds.gcps[0].x, ds.gcps[0].y)
                (0.0, 10.0, 50.0)

                ```
        """
        return [GroundControlPoint.from_gdal(gcp) for gcp in self._ds.raster.GetGCPs()]

    @property
    def gcp_count(self: Georef) -> int:
        """Number of ground-control points attached (``0`` when none).

        Returns:
            int: the GCP count.
        """
        return cast(int, self._ds.raster.GetGCPCount())

    @property
    def gcp_projection(self: Georef) -> str | None:
        """WKT of the GCPs' CRS, or ``None`` when the dataset has no GCPs.

        Returns:
            str | None: the GCP-projection WKT, else ``None``.
        """
        wkt = self._ds.raster.GetGCPProjection()
        return wkt if self._ds.raster.GetGCPCount() else None

    @property
    def has_gcps(self: Georef) -> bool:
        """``True`` when the dataset carries at least one ground-control point.

        Returns:
            bool: whether any GCP is attached.
        """
        return cast(int, self._ds.raster.GetGCPCount()) > 0

    @property
    def rpcs(self: Georef) -> dict[str, str] | None:
        """The dataset's rational-polynomial coefficients, or ``None`` if absent.

        RPCs live in GDAL's ``"RPC"`` metadata domain — a vendor sensor model of
        ~90 string coefficients (``LINE_NUM_COEFF``, ``HEIGHT_OFF``, ...) shipped
        with raw high-resolution satellite imagery.

        Returns:
            dict[str, str] | None: the RPC coefficient mapping, or ``None`` when
            the dataset has no RPC metadata.

        Examples:
            - A plain raster has no RPCs:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> ds = Dataset.create_from_array(
                ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
                ... )
                >>> ds.rpcs is None
                True

                ```
        """
        metadata = self._ds.raster.GetMetadata("RPC")
        return metadata or None

    @property
    def has_rpcs(self: Georef) -> bool:
        """``True`` when the dataset carries RPC metadata.

        Returns:
            bool: whether an RPC sensor model is attached.
        """
        return bool(self._ds.raster.GetMetadata("RPC"))

    def set_rpcs(self: Georef, rpc: Mapping[str, str | float]) -> None:
        """Attach rational-polynomial coefficients (an RPC sensor model).

        Values are stringified before being written to GDAL's ``"RPC"`` metadata
        domain. The dataset must be writable.

        Args:
            rpc: The RPC coefficient mapping. Must contain every required key:
                the five ``*_OFF`` offsets, the five ``*_SCALE`` scales, and the
                four coefficient vectors (``LINE_NUM_COEFF``, ``LINE_DEN_COEFF``,
                ``SAMP_NUM_COEFF``, ``SAMP_DEN_COEFF``).

        Raises:
            ReadOnlyError: The dataset is opened read-only.
            ValueError: One or more required RPC keys are missing (the error
                lists them).

        Examples:
            - Round-trip a coefficient set through the RPC domain:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> rpc = {k: "0" for k in (
                ...     "LINE_OFF", "SAMP_OFF", "LAT_OFF", "LONG_OFF", "HEIGHT_OFF",
                ...     "LINE_SCALE", "SAMP_SCALE", "LAT_SCALE", "LONG_SCALE",
                ...     "HEIGHT_SCALE", "LINE_NUM_COEFF", "LINE_DEN_COEFF",
                ...     "SAMP_NUM_COEFF", "SAMP_DEN_COEFF",
                ... )}
                >>> ds = Dataset.create_from_array(
                ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
                ... )
                >>> ds.set_rpcs(rpc)
                >>> ds.rpcs["HEIGHT_OFF"]
                '0'

                ```
        """
        if self._ds.access == "read_only":
            raise ReadOnlyError(
                "The Dataset is opened read-only. Please read the dataset using "
                "read_only=False to attach RPC metadata."
            )
        missing = _REQUIRED_RPC_KEYS - set(rpc)
        if missing:
            raise ValueError(
                f"RPC metadata is missing required keys: {sorted(missing)}"
            )
        stringified = {key: str(value) for key, value in rpc.items()}
        self._ds.raster.SetMetadata(stringified, "RPC")

    def orthorectify(
        self: Georef,
        *,
        dem: str | Path | Dataset | None = None,
        to_epsg: int | str | None = None,
        method: str = "bilinear",
        rpc_height: float | None = None,
        cell_size: float | None = None,
        lazy: bool = False,
    ) -> Dataset:
        """Orthorectify the dataset from its RPC sensor model onto a map grid.

        Uses the attached rational-polynomial coefficients (attach them first
        with :meth:`set_rpcs`) and, ideally, a digital elevation model to remove
        terrain-induced distortion and produce a map-projected raster.

        Args:
            dem: Elevation model used to evaluate the RPCs — a raster path or a
                :class:`Dataset` (an in-memory dataset is staged to ``/vsimem/``).
                ``None`` falls back to a constant height.
            to_epsg: Target CRS. ``None`` keeps the RPCs' native geographic CRS.
            method: Resampling method. Default ``"bilinear"``.
            rpc_height: Constant elevation (map units) to use when no ``dem`` is
                given. If both ``dem`` and ``rpc_height`` are ``None`` GDAL uses
                height ``0`` and a ``logging.WARNING`` is emitted.
            cell_size: Output pixel size in target-CRS units. ``None`` lets GDAL
                pick.
            lazy: ``True`` returns a VRT-backed view; ``False`` (default)
                materialises the result.

        Returns:
            Dataset: the orthorectified raster.

        Raises:
            ValueError: the dataset has no RPC metadata.
            RuntimeError: GDAL could not build the warp.
        """
        if not self._ds.raster.GetMetadata("RPC"):
            raise ValueError("dataset has no RPC metadata; call set_rpcs first.")
        transformer_options = ["METHOD=RPC"]
        dem_path = self._resolve_dem_path(dem)
        if dem_path is not None:
            transformer_options.append(f"RPC_DEM={dem_path}")
        elif rpc_height is not None:
            transformer_options.append(f"RPC_HEIGHT={rpc_height}")
        else:
            logger.warning(
                "orthorectify: no DEM and no rpc_height given; GDAL will use "
                "height 0, which is rarely correct over real terrain."
            )
        warp_kwargs: dict = {
            "format": "VRT" if lazy else "MEM",
            "resampleAlg": resolve_resampling(method),
            "xRes": cell_size,
            "yRes": cell_size,
            "transformerOptions": transformer_options,
        }
        if to_epsg is not None:
            warp_kwargs["dstSRS"] = _dst_srs_arg(sr_from_user_input(to_epsg))
        dst = gdal.Warp("", self._ds.raster, options=gdal.WarpOptions(**warp_kwargs))
        if dst is None:
            raise RuntimeError("GDAL could not orthorectify the dataset.")
        result = self._ds.__class__(dst, access="read_only")
        if lazy:
            result._warp_source = self._ds
        elif dem_path is not None and dem_path.startswith("/vsimem/orthorectify_dem_"):
            # A materialised result no longer references the staged DEM, so free
            # the /vsimem copy we made from a MEM Dataset (a lazy result keeps it).
            gdal.Unlink(dem_path)
        return result

    @staticmethod
    def _resolve_dem_path(dem: str | Path | Dataset | None) -> str | None:
        """Resolve a DEM argument to a GDAL-openable path.

        Args:
            dem: ``None``, a filesystem path, or a Dataset. A file-backed Dataset
                yields its path; an in-memory one is staged to ``/vsimem/``.

        Returns:
            str | None: the path GDAL's ``RPC_DEM`` should open, or ``None``.
        """
        if dem is None:
            result = None
        elif isinstance(dem, (str, Path)):
            result = str(dem)
        else:
            description = dem.raster.GetDescription()
            if description:
                result = description
            else:
                result = f"/vsimem/orthorectify_dem_{uuid4().hex}.tif"
                gdal.Translate(result, dem.raster)
        return result

    def set_gcps(
        self: Georef,
        gcps: Sequence[GroundControlPoint],
        projection: int | str,
    ) -> None:
        """Attach ground-control points (and their CRS) to the dataset.

        Replaces any existing GCPs. The dataset must be opened writable
        (``read_only=False``); a MEM-backed dataset (e.g. from
        :meth:`Dataset.create_from_array`) is always writable.

        Args:
            gcps: One or more :class:`GroundControlPoint` tie points.
            projection: The GCPs' CRS, in any form
                :func:`pyramids.base.crs.sr_from_user_input` accepts (EPSG int,
                ``"EPSG:4326"``, WKT, PROJ4, ...).

        Raises:
            ReadOnlyError: The dataset is opened read-only.
            ValueError: ``gcps`` is empty.

        Examples:
            - Attach four corner points in EPSG:4326 to an in-memory raster:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from pyramids.dataset._gcp import GroundControlPoint
                >>> ds = Dataset.create_from_array(
                ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
                ... )
                >>> pts = [
                ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
                ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
                ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
                ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
                ... ]
                >>> ds.set_gcps(pts, 4326)
                >>> ds.raster.GetGCPCount()
                4

                ```
        """
        if self._ds.access == "read_only":
            raise ReadOnlyError(
                "The Dataset is opened read-only. Please read the dataset using "
                "read_only=False to attach GCPs."
            )
        gcp_list = list(gcps)
        if not gcp_list:
            raise ValueError("set_gcps requires at least one GroundControlPoint.")
        wkt = sr_from_user_input(projection).ExportToWkt()
        self._ds.raster.SetGCPs([point.to_gdal() for point in gcp_list], wkt)

    def georeference(
        self: Georef,
        *,
        to_epsg: int | str | None = None,
        method: str = "nearest neighbor",
        transform: str = "polynomial",
        order: int = 1,
        cell_size: float | None = None,
        lazy: bool = False,
    ) -> Dataset:
        """Warp the dataset **from its GCPs** into an affine-geotransform raster.

        Fits a transform through the attached ground-control points and resamples
        the pixels onto a regular grid — the step that turns a GCP-tagged scan or
        mosaic into a normal georeferenced raster. The GCPs are read from the
        source automatically (attach them first with :meth:`set_gcps`).

        Args:
            to_epsg: Target CRS for the output. ``None`` warps into the GCPs' own
                CRS; otherwise reprojects to ``to_epsg`` in the same pass (any
                form :func:`pyramids.base.crs.sr_from_user_input` accepts).
            method: Resampling method (see :meth:`Spatial.to_crs`). Default
                ``"nearest neighbor"``.
            transform: ``"polynomial"`` (default) fits a polynomial of degree
                ``order``; ``"tps"`` fits a thin-plate spline (good for many,
                irregularly-spaced points / local warps).
            order: Polynomial degree, one of ``1``/``2``/``3``. Ignored when
                ``transform="tps"``.
            cell_size: Output pixel size in target-CRS units (both axes). ``None``
                lets GDAL pick a size that preserves the source resolution.
            lazy: ``True`` returns a VRT-backed view (pixels warped per read);
                ``False`` (default) materialises the result in memory.

        Returns:
            Dataset: the georeferenced raster.

        Raises:
            ValueError: the dataset has no GCPs, or ``transform``/``order`` is
                invalid.
            RuntimeError: GDAL could not build the warp.

        Examples:
            - Georeference an 8x8 image from four corner GCPs in EPSG:4326:
                ```python
                >>> import numpy as np
                >>> from pyramids.dataset import Dataset
                >>> from pyramids.dataset._gcp import GroundControlPoint
                >>> ds = Dataset.create_from_array(
                ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
                ... )
                >>> ds.set_gcps([
                ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
                ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
                ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
                ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
                ... ], 4326)
                >>> out = ds.georeference()
                >>> out.epsg
                4326

                ```
        """
        if self._ds.raster.GetGCPCount() == 0:
            raise ValueError("dataset has no GCPs; call set_gcps first.")
        if transform not in {"polynomial", "tps"}:
            raise ValueError(
                f"transform must be 'polynomial' or 'tps', got {transform!r}."
            )
        if transform == "polynomial" and order not in (1, 2, 3):
            raise ValueError(f"order must be 1, 2, or 3, got {order!r}.")
        # Force the GCP transformer: GDAL otherwise prefers an affine geotransform
        # when the source carries one, ignoring the GCPs. METHOD=GCP_* makes the
        # warp fit the points regardless.
        if transform == "tps":
            transformer_options = ["METHOD=GCP_TPS"]
        else:
            transformer_options = ["METHOD=GCP_POLYNOMIAL", f"MAX_GCP_ORDER={order}"]
        warp_kwargs: dict = {
            "format": "VRT" if lazy else "MEM",
            "resampleAlg": resolve_resampling(method),
            "xRes": cell_size,
            "yRes": cell_size,
            "transformerOptions": transformer_options,
        }
        if to_epsg is not None:
            warp_kwargs["dstSRS"] = _dst_srs_arg(sr_from_user_input(to_epsg))
        dst = gdal.Warp("", self._ds.raster, options=gdal.WarpOptions(**warp_kwargs))
        if dst is None:
            raise RuntimeError("GDAL could not warp the dataset from its GCPs.")
        result = self._ds.__class__(dst, access="read_only")
        if lazy:
            # The VRT references the source GDAL handle; pin the source Dataset so
            # it cannot be garbage-collected underneath the view.
            result._warp_source = self._ds
        return result

gcp_count property

Number of ground-control points attached (0 when none).

Returns:

Name	Type	Description
int	int	the GCP count.

gcp_projection property

WKT of the GCPs' CRS, or None when the dataset has no GCPs.

Returns:

Type	Description
str | None	str | None: the GCP-projection WKT, else None.

gcps property

The dataset's ground-control points (empty list when it has none).

Returns:

Type	Description
list[GroundControlPoint]	list[GroundControlPoint]: one per attached GCP, in GDAL order.

Examples:

• Read back the points attached with :meth:set_gcps:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> ds = Dataset.create_from_array(
  ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
  ... )
  >>> ds.set_gcps([GroundControlPoint(row=0, col=0, x=10.0, y=50.0)], 4326)
  >>> (ds.gcps[0].col, ds.gcps[0].x, ds.gcps[0].y)
  (0.0, 10.0, 50.0)
  

has_gcps property

True when the dataset carries at least one ground-control point.

Returns:

Name	Type	Description
bool	bool	whether any GCP is attached.

has_rpcs property

True when the dataset carries RPC metadata.

Returns:

Name	Type	Description
bool	bool	whether an RPC sensor model is attached.

rpcs property

The dataset's rational-polynomial coefficients, or None if absent.

RPCs live in GDAL's "RPC" metadata domain — a vendor sensor model of ~90 string coefficients (LINE_NUM_COEFF, HEIGHT_OFF, ...) shipped with raw high-resolution satellite imagery.

Returns:

Type	Description
dict[str, str] | None	dict[str, str] | None: the RPC coefficient mapping, or None when
dict[str, str] | None	the dataset has no RPC metadata.

Examples:

• A plain raster has no RPCs:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> ds = Dataset.create_from_array(
  ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
  ... )
  >>> ds.rpcs is None
  True
  

georeference(*, to_epsg=None, method='nearest neighbor', transform='polynomial', order=1, cell_size=None, lazy=False)

Warp the dataset from its GCPs into an affine-geotransform raster.

Fits a transform through the attached ground-control points and resamples the pixels onto a regular grid — the step that turns a GCP-tagged scan or mosaic into a normal georeferenced raster. The GCPs are read from the source automatically (attach them first with :meth:set_gcps).

Parameters:

Name	Type	Description	Default
to_epsg	int | str | None	Target CRS for the output. None warps into the GCPs' own CRS; otherwise reprojects to to_epsg in the same pass (any form :func:pyramids.base.crs.sr_from_user_input accepts).	None
method	str	Resampling method (see :meth:Spatial.to_crs). Default "nearest neighbor".	'nearest neighbor'
transform	str	"polynomial" (default) fits a polynomial of degree order; "tps" fits a thin-plate spline (good for many, irregularly-spaced points / local warps).	'polynomial'
order	int	Polynomial degree, one of 1/2/3. Ignored when transform="tps".	1
cell_size	float | None	Output pixel size in target-CRS units (both axes). None lets GDAL pick a size that preserves the source resolution.	None
lazy	bool	True returns a VRT-backed view (pixels warped per read); False (default) materialises the result in memory.	False

Returns:

Name	Type	Description
Dataset	Dataset	the georeferenced raster.

Raises:

Type	Description
ValueError	the dataset has no GCPs, or transform/order is invalid.
RuntimeError	GDAL could not build the warp.

Examples:

• Georeference an 8x8 image from four corner GCPs in EPSG:4326:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> ds = Dataset.create_from_array(
  ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
  ... )
  >>> ds.set_gcps([
  ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
  ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
  ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
  ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
  ... ], 4326)
  >>> out = ds.georeference()
  >>> out.epsg
  4326
  

Source code in src/pyramids/dataset/engines/georef.py

def georeference(
    self: Georef,
    *,
    to_epsg: int | str | None = None,
    method: str = "nearest neighbor",
    transform: str = "polynomial",
    order: int = 1,
    cell_size: float | None = None,
    lazy: bool = False,
) -> Dataset:
    """Warp the dataset **from its GCPs** into an affine-geotransform raster.

    Fits a transform through the attached ground-control points and resamples
    the pixels onto a regular grid — the step that turns a GCP-tagged scan or
    mosaic into a normal georeferenced raster. The GCPs are read from the
    source automatically (attach them first with :meth:`set_gcps`).

    Args:
        to_epsg: Target CRS for the output. ``None`` warps into the GCPs' own
            CRS; otherwise reprojects to ``to_epsg`` in the same pass (any
            form :func:`pyramids.base.crs.sr_from_user_input` accepts).
        method: Resampling method (see :meth:`Spatial.to_crs`). Default
            ``"nearest neighbor"``.
        transform: ``"polynomial"`` (default) fits a polynomial of degree
            ``order``; ``"tps"`` fits a thin-plate spline (good for many,
            irregularly-spaced points / local warps).
        order: Polynomial degree, one of ``1``/``2``/``3``. Ignored when
            ``transform="tps"``.
        cell_size: Output pixel size in target-CRS units (both axes). ``None``
            lets GDAL pick a size that preserves the source resolution.
        lazy: ``True`` returns a VRT-backed view (pixels warped per read);
            ``False`` (default) materialises the result in memory.

    Returns:
        Dataset: the georeferenced raster.

    Raises:
        ValueError: the dataset has no GCPs, or ``transform``/``order`` is
            invalid.
        RuntimeError: GDAL could not build the warp.

    Examples:
        - Georeference an 8x8 image from four corner GCPs in EPSG:4326:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> ds = Dataset.create_from_array(
            ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
            ... )
            >>> ds.set_gcps([
            ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
            ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
            ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
            ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
            ... ], 4326)
            >>> out = ds.georeference()
            >>> out.epsg
            4326

            ```
    """
    if self._ds.raster.GetGCPCount() == 0:
        raise ValueError("dataset has no GCPs; call set_gcps first.")
    if transform not in {"polynomial", "tps"}:
        raise ValueError(
            f"transform must be 'polynomial' or 'tps', got {transform!r}."
        )
    if transform == "polynomial" and order not in (1, 2, 3):
        raise ValueError(f"order must be 1, 2, or 3, got {order!r}.")
    # Force the GCP transformer: GDAL otherwise prefers an affine geotransform
    # when the source carries one, ignoring the GCPs. METHOD=GCP_* makes the
    # warp fit the points regardless.
    if transform == "tps":
        transformer_options = ["METHOD=GCP_TPS"]
    else:
        transformer_options = ["METHOD=GCP_POLYNOMIAL", f"MAX_GCP_ORDER={order}"]
    warp_kwargs: dict = {
        "format": "VRT" if lazy else "MEM",
        "resampleAlg": resolve_resampling(method),
        "xRes": cell_size,
        "yRes": cell_size,
        "transformerOptions": transformer_options,
    }
    if to_epsg is not None:
        warp_kwargs["dstSRS"] = _dst_srs_arg(sr_from_user_input(to_epsg))
    dst = gdal.Warp("", self._ds.raster, options=gdal.WarpOptions(**warp_kwargs))
    if dst is None:
        raise RuntimeError("GDAL could not warp the dataset from its GCPs.")
    result = self._ds.__class__(dst, access="read_only")
    if lazy:
        # The VRT references the source GDAL handle; pin the source Dataset so
        # it cannot be garbage-collected underneath the view.
        result._warp_source = self._ds
    return result

orthorectify(*, dem=None, to_epsg=None, method='bilinear', rpc_height=None, cell_size=None, lazy=False)

Orthorectify the dataset from its RPC sensor model onto a map grid.

Uses the attached rational-polynomial coefficients (attach them first with :meth:set_rpcs) and, ideally, a digital elevation model to remove terrain-induced distortion and produce a map-projected raster.

Parameters:

Name	Type	Description	Default
dem	str | Path | Dataset | None	Elevation model used to evaluate the RPCs — a raster path or a :class:Dataset (an in-memory dataset is staged to /vsimem/). None falls back to a constant height.	None
to_epsg	int | str | None	Target CRS. None keeps the RPCs' native geographic CRS.	None
method	str	Resampling method. Default "bilinear".	'bilinear'
rpc_height	float | None	Constant elevation (map units) to use when no dem is given. If both dem and rpc_height are None GDAL uses height 0 and a logging.WARNING is emitted.	None
cell_size	float | None	Output pixel size in target-CRS units. None lets GDAL pick.	None
lazy	bool	True returns a VRT-backed view; False (default) materialises the result.	False

Returns:

Name	Type	Description
Dataset	Dataset	the orthorectified raster.

Raises:

Type	Description
ValueError	the dataset has no RPC metadata.
RuntimeError	GDAL could not build the warp.

Source code in src/pyramids/dataset/engines/georef.py

def orthorectify(
    self: Georef,
    *,
    dem: str | Path | Dataset | None = None,
    to_epsg: int | str | None = None,
    method: str = "bilinear",
    rpc_height: float | None = None,
    cell_size: float | None = None,
    lazy: bool = False,
) -> Dataset:
    """Orthorectify the dataset from its RPC sensor model onto a map grid.

    Uses the attached rational-polynomial coefficients (attach them first
    with :meth:`set_rpcs`) and, ideally, a digital elevation model to remove
    terrain-induced distortion and produce a map-projected raster.

    Args:
        dem: Elevation model used to evaluate the RPCs — a raster path or a
            :class:`Dataset` (an in-memory dataset is staged to ``/vsimem/``).
            ``None`` falls back to a constant height.
        to_epsg: Target CRS. ``None`` keeps the RPCs' native geographic CRS.
        method: Resampling method. Default ``"bilinear"``.
        rpc_height: Constant elevation (map units) to use when no ``dem`` is
            given. If both ``dem`` and ``rpc_height`` are ``None`` GDAL uses
            height ``0`` and a ``logging.WARNING`` is emitted.
        cell_size: Output pixel size in target-CRS units. ``None`` lets GDAL
            pick.
        lazy: ``True`` returns a VRT-backed view; ``False`` (default)
            materialises the result.

    Returns:
        Dataset: the orthorectified raster.

    Raises:
        ValueError: the dataset has no RPC metadata.
        RuntimeError: GDAL could not build the warp.
    """
    if not self._ds.raster.GetMetadata("RPC"):
        raise ValueError("dataset has no RPC metadata; call set_rpcs first.")
    transformer_options = ["METHOD=RPC"]
    dem_path = self._resolve_dem_path(dem)
    if dem_path is not None:
        transformer_options.append(f"RPC_DEM={dem_path}")
    elif rpc_height is not None:
        transformer_options.append(f"RPC_HEIGHT={rpc_height}")
    else:
        logger.warning(
            "orthorectify: no DEM and no rpc_height given; GDAL will use "
            "height 0, which is rarely correct over real terrain."
        )
    warp_kwargs: dict = {
        "format": "VRT" if lazy else "MEM",
        "resampleAlg": resolve_resampling(method),
        "xRes": cell_size,
        "yRes": cell_size,
        "transformerOptions": transformer_options,
    }
    if to_epsg is not None:
        warp_kwargs["dstSRS"] = _dst_srs_arg(sr_from_user_input(to_epsg))
    dst = gdal.Warp("", self._ds.raster, options=gdal.WarpOptions(**warp_kwargs))
    if dst is None:
        raise RuntimeError("GDAL could not orthorectify the dataset.")
    result = self._ds.__class__(dst, access="read_only")
    if lazy:
        result._warp_source = self._ds
    elif dem_path is not None and dem_path.startswith("/vsimem/orthorectify_dem_"):
        # A materialised result no longer references the staged DEM, so free
        # the /vsimem copy we made from a MEM Dataset (a lazy result keeps it).
        gdal.Unlink(dem_path)
    return result

set_gcps(gcps, projection)

Attach ground-control points (and their CRS) to the dataset.

Replaces any existing GCPs. The dataset must be opened writable (read_only=False); a MEM-backed dataset (e.g. from :meth:Dataset.create_from_array) is always writable.

Parameters:

Name	Type	Description	Default
gcps	Sequence[GroundControlPoint]	One or more :class:GroundControlPoint tie points.	required
projection	int | str	The GCPs' CRS, in any form :func:pyramids.base.crs.sr_from_user_input accepts (EPSG int, "EPSG:4326", WKT, PROJ4, ...).	required

Raises:

Type	Description
ReadOnlyError	The dataset is opened read-only.
ValueError	gcps is empty.

Examples:

• Attach four corner points in EPSG:4326 to an in-memory raster:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> ds = Dataset.create_from_array(
  ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
  ... )
  >>> pts = [
  ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
  ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
  ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
  ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
  ... ]
  >>> ds.set_gcps(pts, 4326)
  >>> ds.raster.GetGCPCount()
  4
  

Source code in src/pyramids/dataset/engines/georef.py

def set_gcps(
    self: Georef,
    gcps: Sequence[GroundControlPoint],
    projection: int | str,
) -> None:
    """Attach ground-control points (and their CRS) to the dataset.

    Replaces any existing GCPs. The dataset must be opened writable
    (``read_only=False``); a MEM-backed dataset (e.g. from
    :meth:`Dataset.create_from_array`) is always writable.

    Args:
        gcps: One or more :class:`GroundControlPoint` tie points.
        projection: The GCPs' CRS, in any form
            :func:`pyramids.base.crs.sr_from_user_input` accepts (EPSG int,
            ``"EPSG:4326"``, WKT, PROJ4, ...).

    Raises:
        ReadOnlyError: The dataset is opened read-only.
        ValueError: ``gcps`` is empty.

    Examples:
        - Attach four corner points in EPSG:4326 to an in-memory raster:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> ds = Dataset.create_from_array(
            ...     np.ones((8, 8), "float32"), top_left_corner=(0.0, 8.0), cell_size=1.0
            ... )
            >>> pts = [
            ...     GroundControlPoint(row=0, col=0, x=10.0, y=50.0),
            ...     GroundControlPoint(row=0, col=8, x=11.0, y=50.0),
            ...     GroundControlPoint(row=8, col=0, x=10.0, y=49.0),
            ...     GroundControlPoint(row=8, col=8, x=11.0, y=49.0),
            ... ]
            >>> ds.set_gcps(pts, 4326)
            >>> ds.raster.GetGCPCount()
            4

            ```
    """
    if self._ds.access == "read_only":
        raise ReadOnlyError(
            "The Dataset is opened read-only. Please read the dataset using "
            "read_only=False to attach GCPs."
        )
    gcp_list = list(gcps)
    if not gcp_list:
        raise ValueError("set_gcps requires at least one GroundControlPoint.")
    wkt = sr_from_user_input(projection).ExportToWkt()
    self._ds.raster.SetGCPs([point.to_gdal() for point in gcp_list], wkt)

set_rpcs(rpc)

Attach rational-polynomial coefficients (an RPC sensor model).

Values are stringified before being written to GDAL's "RPC" metadata domain. The dataset must be writable.

Parameters:

Name	Type	Description	Default
rpc	Mapping[str, str | float]	The RPC coefficient mapping. Must contain every required key: the five *_OFF offsets, the five *_SCALE scales, and the four coefficient vectors (LINE_NUM_COEFF, LINE_DEN_COEFF, SAMP_NUM_COEFF, SAMP_DEN_COEFF).	required

Raises:

Type	Description
ReadOnlyError	The dataset is opened read-only.
ValueError	One or more required RPC keys are missing (the error lists them).

Examples:

• Round-trip a coefficient set through the RPC domain:

  >>> import numpy as np
  >>> from pyramids.dataset import Dataset
  >>> rpc = {k: "0" for k in (
  ...     "LINE_OFF", "SAMP_OFF", "LAT_OFF", "LONG_OFF", "HEIGHT_OFF",
  ...     "LINE_SCALE", "SAMP_SCALE", "LAT_SCALE", "LONG_SCALE",
  ...     "HEIGHT_SCALE", "LINE_NUM_COEFF", "LINE_DEN_COEFF",
  ...     "SAMP_NUM_COEFF", "SAMP_DEN_COEFF",
  ... )}
  >>> ds = Dataset.create_from_array(
  ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
  ... )
  >>> ds.set_rpcs(rpc)
  >>> ds.rpcs["HEIGHT_OFF"]
  '0'
  

Source code in src/pyramids/dataset/engines/georef.py

def set_rpcs(self: Georef, rpc: Mapping[str, str | float]) -> None:
    """Attach rational-polynomial coefficients (an RPC sensor model).

    Values are stringified before being written to GDAL's ``"RPC"`` metadata
    domain. The dataset must be writable.

    Args:
        rpc: The RPC coefficient mapping. Must contain every required key:
            the five ``*_OFF`` offsets, the five ``*_SCALE`` scales, and the
            four coefficient vectors (``LINE_NUM_COEFF``, ``LINE_DEN_COEFF``,
            ``SAMP_NUM_COEFF``, ``SAMP_DEN_COEFF``).

    Raises:
        ReadOnlyError: The dataset is opened read-only.
        ValueError: One or more required RPC keys are missing (the error
            lists them).

    Examples:
        - Round-trip a coefficient set through the RPC domain:
            ```python
            >>> import numpy as np
            >>> from pyramids.dataset import Dataset
            >>> rpc = {k: "0" for k in (
            ...     "LINE_OFF", "SAMP_OFF", "LAT_OFF", "LONG_OFF", "HEIGHT_OFF",
            ...     "LINE_SCALE", "SAMP_SCALE", "LAT_SCALE", "LONG_SCALE",
            ...     "HEIGHT_SCALE", "LINE_NUM_COEFF", "LINE_DEN_COEFF",
            ...     "SAMP_NUM_COEFF", "SAMP_DEN_COEFF",
            ... )}
            >>> ds = Dataset.create_from_array(
            ...     np.ones((4, 4), "float32"), top_left_corner=(0.0, 4.0), cell_size=1.0
            ... )
            >>> ds.set_rpcs(rpc)
            >>> ds.rpcs["HEIGHT_OFF"]
            '0'

            ```
    """
    if self._ds.access == "read_only":
        raise ReadOnlyError(
            "The Dataset is opened read-only. Please read the dataset using "
            "read_only=False to attach RPC metadata."
        )
    missing = _REQUIRED_RPC_KEYS - set(rpc)
    if missing:
        raise ValueError(
            f"RPC metadata is missing required keys: {sorted(missing)}"
        )
    stringified = {key: str(value) for key, value in rpc.items()}
    self._ds.raster.SetMetadata(stringified, "RPC")

pyramids.dataset._gcp.GroundControlPoint dataclass

An immutable ground-control point: a pixel mapped to a map coordinate.

A GCP states "pixel (col, row) is at map coordinate (x, y[, z])". The pixel axes follow GDAL: col is the pixel (x) index and row is the line (y) index, both measured from the raster's top-left corner.

Parameters:

Name	Type	Description	Default
row	float	Pixel-space line (y / row) index of the control point.	required
col	float	Pixel-space pixel (x / column) index of the control point.	required
x	float	Map-space X (easting / longitude) at that pixel.	required
y	float	Map-space Y (northing / latitude) at that pixel.	required
z	float	Map-space Z (elevation). Defaults to 0.0.	0.0
id	str | None	Optional point identifier. Defaults to None.	None
info	str | None	Optional free-text label. Defaults to None.	None

Examples:

• Build a point and read its fields:

  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> gcp = GroundControlPoint(row=0.0, col=0.0, x=10.0, y=50.0)
  >>> (gcp.col, gcp.row, gcp.x, gcp.y, gcp.z)
  (0.0, 0.0, 10.0, 50.0, 0.0)
  

• Round-trip through the GDAL representation:

  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> gcp = GroundControlPoint(row=4.0, col=2.0, x=11.5, y=46.2, id="tl")
  >>> back = GroundControlPoint.from_gdal(gcp.to_gdal())
  >>> (back.col, back.row, back.x, back.y, back.id)
  (2.0, 4.0, 11.5, 46.2, 'tl')
  

Source code in src/pyramids/dataset/_gcp.py

@dataclass(frozen=True)
class GroundControlPoint:
    """An immutable ground-control point: a pixel mapped to a map coordinate.

    A GCP states "pixel ``(col, row)`` is at map coordinate ``(x, y[, z])``". The
    pixel axes follow GDAL: ``col`` is the pixel (x) index and ``row`` is the line
    (y) index, both measured from the raster's top-left corner.

    Args:
        row: Pixel-space line (y / row) index of the control point.
        col: Pixel-space pixel (x / column) index of the control point.
        x: Map-space X (easting / longitude) at that pixel.
        y: Map-space Y (northing / latitude) at that pixel.
        z: Map-space Z (elevation). Defaults to ``0.0``.
        id: Optional point identifier. Defaults to ``None``.
        info: Optional free-text label. Defaults to ``None``.

    Examples:
        - Build a point and read its fields:
            ```python
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> gcp = GroundControlPoint(row=0.0, col=0.0, x=10.0, y=50.0)
            >>> (gcp.col, gcp.row, gcp.x, gcp.y, gcp.z)
            (0.0, 0.0, 10.0, 50.0, 0.0)

            ```
        - Round-trip through the GDAL representation:
            ```python
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> gcp = GroundControlPoint(row=4.0, col=2.0, x=11.5, y=46.2, id="tl")
            >>> back = GroundControlPoint.from_gdal(gcp.to_gdal())
            >>> (back.col, back.row, back.x, back.y, back.id)
            (2.0, 4.0, 11.5, 46.2, 'tl')

            ```
    """

    row: float
    col: float
    x: float
    y: float
    z: float = 0.0
    id: str | None = None
    info: str | None = None

    def to_gdal(self) -> gdal.GCP:
        """Convert to an :class:`osgeo.gdal.GCP`.

        Returns:
            gdal.GCP: a GDAL GCP whose ``pixel``/``line`` are this point's
            ``col``/``row`` and whose ``(x, y, z)`` are the map coordinate.

        Examples:
            - The GDAL GCP carries the pixel and map coordinates:
                ```python
                >>> from pyramids.dataset._gcp import GroundControlPoint
                >>> g = GroundControlPoint(row=3.0, col=7.0, x=1.0, y=2.0).to_gdal()
                >>> (g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY)
                (7.0, 3.0, 1.0, 2.0)

                ```
        """
        return gdal.GCP(
            self.x, self.y, self.z, self.col, self.row, self.info or "", self.id or ""
        )

    @classmethod
    def from_gdal(cls, gcp: gdal.GCP) -> GroundControlPoint:
        """Build a :class:`GroundControlPoint` from an :class:`osgeo.gdal.GCP`.

        Args:
            gcp: A GDAL GCP (fields ``GCPPixel``/``GCPLine``/``GCPX``/``GCPY``/
                ``GCPZ``/``Id``/``Info``). Empty ``Id``/``Info`` map to ``None``.

        Returns:
            GroundControlPoint: the equivalent value object.

        Examples:
            - Convert a GDAL GCP into the value object:
                ```python
                >>> from osgeo import gdal
                >>> from pyramids.dataset._gcp import GroundControlPoint
                >>> gp = GroundControlPoint.from_gdal(gdal.GCP(1.0, 2.0, 0.0, 7.0, 3.0))
                >>> (gp.col, gp.row, gp.x, gp.y)
                (7.0, 3.0, 1.0, 2.0)

                ```
        """
        return cls(
            row=gcp.GCPLine,
            col=gcp.GCPPixel,
            x=gcp.GCPX,
            y=gcp.GCPY,
            z=gcp.GCPZ,
            id=gcp.Id or None,
            info=gcp.Info or None,
        )

from_gdal(gcp) classmethod

Build a :class:GroundControlPoint from an :class:osgeo.gdal.GCP.

Parameters:

Name	Type	Description	Default
gcp	GCP	A GDAL GCP (fields GCPPixel/GCPLine/GCPX/GCPY/ GCPZ/Id/Info). Empty Id/Info map to None.	required

Returns:

Name	Type	Description
GroundControlPoint	GroundControlPoint	the equivalent value object.

Examples:

• Convert a GDAL GCP into the value object:

  >>> from osgeo import gdal
  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> gp = GroundControlPoint.from_gdal(gdal.GCP(1.0, 2.0, 0.0, 7.0, 3.0))
  >>> (gp.col, gp.row, gp.x, gp.y)
  (7.0, 3.0, 1.0, 2.0)
  

Source code in src/pyramids/dataset/_gcp.py

@classmethod
def from_gdal(cls, gcp: gdal.GCP) -> GroundControlPoint:
    """Build a :class:`GroundControlPoint` from an :class:`osgeo.gdal.GCP`.

    Args:
        gcp: A GDAL GCP (fields ``GCPPixel``/``GCPLine``/``GCPX``/``GCPY``/
            ``GCPZ``/``Id``/``Info``). Empty ``Id``/``Info`` map to ``None``.

    Returns:
        GroundControlPoint: the equivalent value object.

    Examples:
        - Convert a GDAL GCP into the value object:
            ```python
            >>> from osgeo import gdal
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> gp = GroundControlPoint.from_gdal(gdal.GCP(1.0, 2.0, 0.0, 7.0, 3.0))
            >>> (gp.col, gp.row, gp.x, gp.y)
            (7.0, 3.0, 1.0, 2.0)

            ```
    """
    return cls(
        row=gcp.GCPLine,
        col=gcp.GCPPixel,
        x=gcp.GCPX,
        y=gcp.GCPY,
        z=gcp.GCPZ,
        id=gcp.Id or None,
        info=gcp.Info or None,
    )

to_gdal()

Convert to an :class:osgeo.gdal.GCP.

Returns:

Type	Description
GCP	gdal.GCP: a GDAL GCP whose pixel/line are this point's
GCP	col/row and whose (x, y, z) are the map coordinate.

Examples:

• The GDAL GCP carries the pixel and map coordinates:

  >>> from pyramids.dataset._gcp import GroundControlPoint
  >>> g = GroundControlPoint(row=3.0, col=7.0, x=1.0, y=2.0).to_gdal()
  >>> (g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY)
  (7.0, 3.0, 1.0, 2.0)
  

Source code in src/pyramids/dataset/_gcp.py

def to_gdal(self) -> gdal.GCP:
    """Convert to an :class:`osgeo.gdal.GCP`.

    Returns:
        gdal.GCP: a GDAL GCP whose ``pixel``/``line`` are this point's
        ``col``/``row`` and whose ``(x, y, z)`` are the map coordinate.

    Examples:
        - The GDAL GCP carries the pixel and map coordinates:
            ```python
            >>> from pyramids.dataset._gcp import GroundControlPoint
            >>> g = GroundControlPoint(row=3.0, col=7.0, x=1.0, y=2.0).to_gdal()
            >>> (g.GCPPixel, g.GCPLine, g.GCPX, g.GCPY)
            (7.0, 3.0, 1.0, 2.0)

            ```
    """
    return gdal.GCP(
        self.x, self.y, self.z, self.col, self.row, self.info or "", self.id or ""
    )