raster

Raster data structure.

Bounds

Bases: NamedTuple

Bounding box coordinates for a raster.

Attributes:

Name	Type	Description
xmin	float	The minimum x-coordinate.
ymin	float	The minimum y-coordinate.
xmax	float	The maximum x-coordinate.
ymax	float	The maximum y-coordinate.

Source code in src/rastr/raster.py

class Bounds(NamedTuple):
    """Bounding box coordinates for a raster.

    Attributes:
        xmin: The minimum x-coordinate.
        ymin: The minimum y-coordinate.
        xmax: The maximum x-coordinate.
        ymax: The maximum y-coordinate.
    """

    xmin: float
    ymin: float
    xmax: float
    ymax: float

Raster

Bases: BaseModel

2-dimensional raster and metadata.

Source code in src/rastr/raster.py

class Raster(BaseModel):
    """2-dimensional raster and metadata."""

    arr: InstanceOf[np.ndarray]
    raster_meta: RasterMeta

    @field_validator("arr")
    @classmethod
    def check_2d_array(cls, v: NDArray) -> NDArray:
        """Validator to ensure the cell array is 2D."""
        if v.ndim != 2:
            msg = "Cell array must be 2D"
            raise RasterCellArrayShapeError(msg)
        return v

    @property
    def meta(self) -> RasterMeta:
        """Alias for raster_meta."""
        return self.raster_meta

    @meta.setter
    def meta(self, value: RasterMeta) -> None:
        self.raster_meta = value

    @property
    def shape(self) -> tuple[int, ...]:
        """Shape of the raster array."""
        return self.arr.shape

    @property
    def crs(self) -> CRS:
        """Convenience property to access the CRS via meta."""
        return self.meta.crs

    @crs.setter
    def crs(self, value: CRS) -> None:
        """Set the CRS via meta."""
        self.meta.crs = value

    @property
    def transform(self) -> Affine:
        """Convenience property to access the transform via meta."""
        return self.meta.transform

    @transform.setter
    def transform(self, value: Affine) -> None:
        """Set the transform via meta."""
        self.meta.transform = value

    @property
    def origin(self) -> tuple[float, float]:
        """The grid origin (x, y) — the corner of the first pixel.

        This is the translation component (c, f) of the affine transform.
        For north-up rasters this corresponds to (xmin, ymax); for south-up
        rasters it corresponds to (xmin, ymin).
        """
        return (self.transform.c, self.transform.f)

    @origin.setter
    def origin(self, value: tuple[float, float]) -> None:
        """Set the grid origin via the transform."""
        x, y = value
        t = self.transform
        self.transform = Affine(t.a, t.b, x, t.d, t.e, y)

    @property
    def cell_size(self) -> tuple[float, float]:
        """Convenience property to access the cell size via meta."""
        return self.meta.cell_size

    @property
    def has_square_cells(self) -> bool:
        """Check if the raster has square cells."""
        return self.meta.has_square_cells

    @property
    def square_cell_size(self) -> float:
        """Convenience property to access the square cell size via meta."""
        return self.meta.square_cell_size

    def __init__(
        self,
        *,
        arr: ArrayLike,
        meta: RasterMeta | None = None,
        raster_meta: RasterMeta | None = None,
    ) -> None:
        arr = np.asarray(arr)

        # Set the meta
        if meta is not None and raster_meta is not None:
            msg = (
                "Only one of 'meta' or 'raster_meta' should be provided, they are "
                "aliases."
            )
            raise ValueError(msg)
        elif meta is not None and raster_meta is None:
            raster_meta = meta
        elif meta is None and raster_meta is not None:
            pass
        else:
            # Don't need to mention `'meta'` to simplify the messaging.
            msg = "The attribute 'raster_meta' is required."
            raise ValueError(msg)

        super().__init__(arr=arr, raster_meta=raster_meta)

    def __eq__(self, other: object) -> bool:
        """Check equality of two Raster objects."""
        if not isinstance(other, Raster):
            return NotImplemented
        return (
            np.array_equal(self.arr, other.arr, equal_nan=True)
            and self.raster_meta == other.raster_meta
        )

    def is_like(self, other: Raster) -> bool:
        """Check if two Raster objects have the same metadata and shape.

        Args:
            other: Another Raster to compare with.

        Returns:
            True if both rasters have the same meta and shape attributes.
        """
        return self.meta == other.meta and self.shape == other.shape

    __hash__ = BaseModel.__hash__

    def __add__(self, other: float | Self | ArrayLike) -> Self:
        cls = self.__class__
        if isinstance(other, float | int):
            new_arr = self.arr + other
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        elif isinstance(other, Raster):
            if self.raster_meta != other.raster_meta:
                msg = (
                    "Rasters must have the same metadata (e.g. CRS, cell size, etc.) "
                    "to be added"
                )
                raise ValueError(msg)
            if self.arr.shape != other.arr.shape:
                msg = (
                    "Rasters must have the same shape to be added:\n"
                    f"{self.arr.shape} != {other.arr.shape}"
                )
                raise ValueError(msg)
            new_arr = self.arr + other.arr
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        else:
            return NotImplemented

    def __radd__(self, other: float | ArrayLike) -> Self:
        return self + other

    def __mul__(self, other: float | Self | ArrayLike) -> Self:
        cls = self.__class__
        if isinstance(other, float | int):
            new_arr = self.arr * other
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        elif isinstance(other, Raster):
            if self.raster_meta != other.raster_meta:
                msg = (
                    "Rasters must have the same metadata (e.g. CRS, cell size, etc.) "
                    "to be multiplied"
                )
                raise ValueError(msg)
            if self.arr.shape != other.arr.shape:
                msg = "Rasters must have the same shape to be multiplied"
                raise ValueError(msg)
            new_arr = self.arr * other.arr
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        else:
            return NotImplemented

    def __rmul__(self, other: float | ArrayLike) -> Self:
        return self * other

    def __truediv__(self, other: float | Self | ArrayLike) -> Self:
        cls = self.__class__
        if isinstance(other, float | int):
            new_arr = self.arr / other
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        elif isinstance(other, Raster):
            if self.raster_meta != other.raster_meta:
                msg = (
                    "Rasters must have the same metadata (e.g. CRS, cell size, etc.) "
                    "to be divided"
                )
                raise ValueError(msg)
            if self.arr.shape != other.arr.shape:
                msg = "Rasters must have the same shape to be divided"
                raise ValueError(msg)
            new_arr = self.arr / other.arr
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        else:
            return NotImplemented

    def __rtruediv__(self, other: float) -> Self:
        if not isinstance(other, float | int):
            return NotImplemented
        cls = self.__class__
        new_arr = other / self.arr
        return cls(arr=new_arr, raster_meta=self.raster_meta)

    def __sub__(self, other: float | Self | ArrayLike) -> Self:
        cls = self.__class__
        if isinstance(other, float | int):
            new_arr = self.arr - other
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        elif isinstance(other, Raster):
            if self.raster_meta != other.raster_meta:
                msg = (
                    "Rasters must have the same metadata (e.g. CRS, cell size, etc.) "
                    "to be subtracted"
                )
                raise ValueError(msg)
            if self.arr.shape != other.arr.shape:
                msg = "Rasters must have the same shape to be subtracted"
                raise ValueError(msg)
            new_arr = self.arr - other.arr
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        else:
            return NotImplemented

    def __rsub__(self, other: float | ArrayLike) -> Self:
        cls = self.__class__
        if isinstance(other, float | int):
            new_arr = other - self.arr
            return cls(arr=new_arr, raster_meta=self.raster_meta)
        else:
            return NotImplemented

    def __neg__(self) -> Self:
        cls = self.__class__
        return cls(arr=-self.arr, raster_meta=self.raster_meta)

    def __array_ufunc__(
        self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any
    ) -> Self:
        """Support NumPy ufuncs between ndarrays and Rasters.

        Carries out the ufunc on the underlying array when the ndarray operand
        has the same shape as this Raster.  Returns ``NotImplemented`` for
        unsupported ufuncs, unsupported methods, or shape mismatches.
        """

        # There are over 60 ufuncs and not all output 2D arrays. For example, reduction
        # ops like np.sum and np.min return scalars while some ufuncs can supposedly
        # return tuples (np.modf). Some would also create types with are currently
        # unsupported by rastr, like np.greater creating bool. Therefore, we choose
        # to whitelist supported ufuncs, and aim to gradually expand this list. This
        # approach allows us to fail fast with a clear error message when an unsupported
        # ufunc is used, rather than silently passing and potentially producing
        # unexpected results.
        if method != "__call__" or ufunc not in _RASTER_SUPPORTED_UFUNCS:
            return NotImplemented
        new_inputs = []
        for inp in inputs:
            if isinstance(inp, np.ndarray):
                if inp.shape != self.shape:
                    msg = (
                        f"Cannot apply ufunc '{ufunc.__name__}' between ndarray of"
                        f" shape {inp.shape} and Raster of shape {self.shape}:"
                        " shapes must be equal."
                    )
                    raise ValueError(msg)
                new_inputs.append(inp)
            elif isinstance(inp, Raster):
                new_inputs.append(inp.arr)
            else:
                new_inputs.append(inp)
        cls = self.__class__
        result_arr = ufunc(*new_inputs, **kwargs)
        return cls(arr=result_arr, raster_meta=self.raster_meta)

    def abs(self) -> Self:
        """Compute the absolute value of the raster.

        Returns a new raster with the absolute value of each cell. The original raster
        is not modified.

        Returns:
            A new Raster instance with the absolute values.
        """
        cls = self.__class__
        return cls(arr=np.abs(self.arr), raster_meta=self.raster_meta)

    def log(self) -> Self:
        """Compute the natural logarithm of the raster.

        Returns a new raster with the natural logarithm of each cell. The original
        raster is not modified.

        Returns:
            A new Raster instance with the natural logarithm values.
        """
        cls = self.__class__
        return cls(arr=np.log(self.arr), raster_meta=self.raster_meta)

    def exp(self) -> Self:
        """Compute the exponential of the raster.

        Returns a new raster with the exponential (e^x) of each cell. The original
        raster is not modified.

        Returns:
            A new Raster instance with the exponential values.
        """
        cls = self.__class__
        return cls(arr=np.exp(self.arr), raster_meta=self.raster_meta)

    def clamp(
        self,
        a_min: float | None = None,
        a_max: float | None = None,
    ) -> Self:
        """Clip (clamp) values to a specified range.

        Returns a new raster with values clipped to [a_min, a_max].

        Args:
            a_min: Minimum value. Values below this will be set to a_min.
                   If None, no minimum clipping is applied.
            a_max: Maximum value. Values above this will be set to a_max.
                   If None, no maximum clipping is applied.

        Returns:
            A new Raster instance with clipped values.
        """
        cls = self.__class__
        return cls(
            arr=np.clip(self.arr, a_min, a_max),
            raster_meta=self.raster_meta,
        )

    def astype(self, dtype: DTypeLike) -> Self:
        """Cast the raster array to a specified dtype.

        Returns a new raster with the array cast to the given dtype. The original
        raster is not modified.

        Args:
            dtype: Target data type (e.g. ``"float32"``, ``np.int16``).

        Returns:
            A new Raster instance with the array cast to the specified dtype.
        """
        cls = self.__class__
        return cls(arr=self.arr.astype(dtype), raster_meta=self.raster_meta)

    def set_crs(self, crs: CRS | str, *, allow_override: bool = False) -> Self:
        """Set the CRS of the raster without reprojecting.

        To reproject the raster data to a different CRS, use `to_crs()` instead.

        Args:
            crs: The coordinate reference system to assign. Can be a CRS object or a
                 string that can be parsed by pyproj (e.g., 'EPSG:4326').
            allow_override: If False (default), raises an error if the raster already
                           has a CRS that differs from the one being set. If True,
                           allows overriding any existing CRS.

        Returns:
            A new Raster instance with the updated CRS and unchanged array data.

        Raises:
            ValueError: If the raster already has a CRS that differs from the one
                       being set and allow_override is False.

        """
        crs_obj = CRS.from_user_input(crs)

        # Check if raster already has a different CRS
        if (
            self.raster_meta.crs is not None
            and self.raster_meta.crs != crs_obj
            and not allow_override
        ):
            msg = (
                f"The raster already has a CRS ({self.raster_meta.crs}) which is "
                f"different from the one provided ({crs_obj}). Use "
                f"allow_override=True to override."
            )
            raise ValueError(msg)

        new_meta = RasterMeta(
            crs=crs_obj,
            transform=self.raster_meta.transform,
        )
        return self.__class__(arr=self.arr, raster_meta=new_meta)

    def set_origin(self, *, x: float | None = None, y: float | None = None) -> Self:
        """Set the transform origin without modifying pixel data.

        The origin is the corner of the first pixel in the raster grid —
        the translation component (c, f) of the affine transform. For
        north-up rasters this corresponds to (xmin, ymax); for south-up
        rasters it corresponds to (xmin, ymin).

        Unspecified axes retain their current value. If neither axis is
        provided, returns an unchanged copy.

        Args:
            x: New x-coordinate of the grid origin. If None, keeps current.
            y: New y-coordinate of the grid origin. If None, keeps current.

        Returns:
            A new Raster with the updated transform origin and unchanged array data.

        Example:
            Normalize a raster with non-standard longitude (e.g., -200° to -170°)
            to standard EPSG:4326 range::

                raster = raster.set_origin(x=raster.origin[0] + 360)
        """
        t = self.raster_meta.transform
        new_x = x if x is not None else t.c
        new_y = y if y is not None else t.f
        new_transform = Affine(t.a, t.b, new_x, t.d, t.e, new_y)
        new_meta = RasterMeta(
            crs=self.raster_meta.crs,
            transform=new_transform,
        )
        return self.__class__(arr=self.arr, raster_meta=new_meta)

    @property
    def cell_centre_coords(self) -> NDArray[np.float64]:
        """Get the coordinates of the cell centres in the raster."""
        return self.raster_meta.get_cell_centre_coords(self.arr.shape)

    @property
    def cell_x_coords(self) -> NDArray[np.float64]:
        """Get the x coordinates of the cell centres in the raster."""
        return self.raster_meta.get_cell_x_coords(self.arr.shape[1])

    @property
    def cell_y_coords(self) -> NDArray[np.float64]:
        """Get the y coordinates of the cell centres in the raster."""
        return self.raster_meta.get_cell_y_coords(self.arr.shape[0])

    @contextmanager
    def to_rasterio_dataset(
        self,
    ) -> Generator[DatasetReader | BufferedDatasetWriter | DatasetWriter]:
        """Create a rasterio in-memory dataset from the Raster object.

        Example:
            >>> raster = Raster.example()
            >>> with raster.to_rasterio_dataset() as dataset:
            >>>     ...
        """
        memfile = MemoryFile()

        height, width = self.arr.shape

        try:
            with memfile.open(
                driver="GTiff",
                height=height,
                width=width,
                count=1,  # Assuming a single band; adjust as necessary
                dtype=self.arr.dtype,
                crs=self.raster_meta.crs.to_wkt(),
                transform=self.raster_meta.transform,
            ) as dataset:
                dataset.write(self.arr, 1)

            # Yield the dataset for reading
            with memfile.open() as dataset:
                yield dataset
        finally:
            memfile.close()

    @overload
    def sample(
        self,
        xy: Collection[tuple[float, float]] | Collection[Point] | ArrayLike,
        *,
        na_action: Literal["raise", "ignore"] = "raise",
    ) -> NDArray: ...
    @overload
    def sample(
        self,
        xy: tuple[float, float] | Point,
        *,
        na_action: Literal["raise", "ignore"] = "raise",
    ) -> float: ...
    def sample(
        self,
        xy: Collection[tuple[float, float]]
        | Collection[Point]
        | ArrayLike
        | tuple[float, float]
        | Point
        | gpd.GeoDataFrame,
        *,
        na_action: Literal["raise", "ignore"] = "raise",
    ) -> NDArray | float:
        """Sample raster values at GeoSeries locations and return sampled values.

        Args:
            xy: A list of (x, y) coordinates, shapely Point objects, or a
                GeoDataFrame containing Point geometries to sample the raster at.
            na_action: Action to take when a NaN value is encountered in the input xy.
                       Options are "raise" (raise an error) or "ignore" (replace with
                       NaN).

        Returns:
            A list of sampled raster values for each geometry in the GeoSeries.
        """
        # If this function is too slow, consider the optimizations detailed here:
        # https://rdrn.me/optimising-sampling/

        import geopandas as gpd

        if isinstance(xy, gpd.GeoDataFrame):
            xy = _gdf_to_xy(xy, self.crs)
            singleton = False

        # Convert shapely Points to coordinate tuples if needed
        elif isinstance(xy, Point):
            xy = [(xy.x, xy.y)]
            singleton = True
        elif (
            isinstance(xy, Collection)
            and len(xy) > 0
            and isinstance(next(iter(xy)), Point)
        ):
            xy = [(point.x, point.y) for point in xy]  # pyright: ignore[reportAttributeAccessIssue]
            singleton = False
        elif (
            isinstance(xy, tuple)
            and len(xy) == 2
            and isinstance(next(iter(xy)), (float, int))
        ):
            xy = [xy]  # pyright: ignore[reportAssignmentType]
            singleton = True
        else:
            singleton = False

        xy = np.asarray(xy, dtype=float)

        if len(xy) == 0:
            # Short-circuit
            return np.array([], dtype=float)

        # Create in-memory rasterio dataset from the incumbent Raster object
        with self.to_rasterio_dataset() as dataset:
            if dataset.count != 1:
                msg = "Only single band rasters are supported."
                raise NotImplementedError(msg)

            xy_arr = np.array(xy)

            # Determine the indexes of any x,y coordinates where either is NaN.
            # We will drop these indexes for the purposes of calling .sample, but
            # then we will add NaN values back in at the end, inserting NaN into the
            # results array.
            xy_is_nan = np.isnan(xy_arr).any(axis=1)
            xy_nan_idxs = list(np.atleast_1d(np.squeeze(np.nonzero(xy_is_nan))))
            xy_arr = xy_arr[~xy_is_nan]

            if na_action == "raise" and len(xy_nan_idxs) > 0:
                nan_error_msg = "NaN value found in input coordinates"
                raise ValueError(nan_error_msg)

            # Sample the raster in-memory dataset (e.g. PGA values) at the coordinates
            samples = list(
                rasterio.sample.sample_gen(
                    dataset,
                    xy_arr,
                    indexes=1,  # Single band raster, N.B. rasterio is 1-indexed
                    masked=True,
                )
            )

            # Convert the sampled values to a NumPy array and set masked values to NaN
            raster_values = np.array(
                [s.data[0] if not numpy.ma.getmask(s) else np.nan for s in samples]
            ).astype(float)

            if len(xy_nan_idxs) > 0:
                # Insert NaN values back into the results array
                # This is tricky because all the indexes get offset once we remove
                # elements.
                offset_xy_nan_idxs = xy_nan_idxs - np.arange(len(xy_nan_idxs))
                raster_values = np.insert(
                    raster_values,
                    offset_xy_nan_idxs,
                    np.nan,
                    axis=0,
                )

        if singleton:
            (raster_value,) = raster_values
            return raster_value

        return raster_values

    @property
    def bounds(self) -> Bounds:
        """Bounding box of the raster as a named tuple with xmin, ymin, xmax, ymax.

        The bounds represent the outer edges of the raster cells.
        """
        x1, y1, x2, y2 = rasterio.transform.array_bounds(
            height=self.arr.shape[0],
            width=self.arr.shape[1],
            transform=self.raster_meta.transform,
        )
        xmin, xmax = sorted([x1, x2])
        ymin, ymax = sorted([y1, y2])
        return Bounds(
            xmin=float(xmin), ymin=float(ymin), xmax=float(xmax), ymax=float(ymax)
        )

    @property
    def bbox(self) -> Polygon:
        """Bounding box of the raster as a shapely polygon."""
        xmin, ymin, xmax, ymax = self.bounds
        return Polygon(
            [
                (xmin, ymin),
                (xmin, ymax),
                (xmax, ymax),
                (xmax, ymin),
                (xmin, ymin),
            ]
        )

    def explore(  # noqa: PLR0913 c.f. geopandas.explore which also has many input args
        self,
        *,
        m: Map | None = None,
        opacity: float = 1.0,
        colormap: str
        | Callable[[float], tuple[float, float, float, float]] = "viridis",
        cbar_label: str | None = None,
        vmin: float | None = None,
        vmax: float | None = None,
    ) -> Map:
        """Display the raster on a folium map."""
        if not FOLIUM_INSTALLED:
            msg = "The 'folium' package is required for 'explore()'."
            raise ImportError(msg)

        import folium.raster_layers

        if m is None:
            m = folium.Map()

        if vmin is not None and vmax is not None and vmax <= vmin:
            msg = "'vmin' must be less than 'vmax'."
            raise ValueError(msg)

        cmap_func = _get_colormap_function(colormap)

        # Transform bounds to WGS84 using pyproj directly
        wgs84_crs = CRS.from_epsg(4326)
        transformer = Transformer.from_crs(
            self.raster_meta.crs, wgs84_crs, always_xy=True
        )

        # Get the corner points of the bounding box
        raster_xmin, raster_ymin, raster_xmax, raster_ymax = self.bounds
        corner_points = [
            (raster_xmin, raster_ymin),
            (raster_xmin, raster_ymax),
            (raster_xmax, raster_ymax),
            (raster_xmax, raster_ymin),
        ]

        # Transform all corner points to WGS84
        transformed_points = [transformer.transform(x, y) for x, y in corner_points]

        # Find the bounding box of the transformed points
        transformed_xs, transformed_ys = zip(*transformed_points, strict=True)
        xmin, xmax = min(transformed_xs), max(transformed_xs)
        ymin, ymax = min(transformed_ys), max(transformed_ys)

        # Normalize the array to [0, 1] for colormap mapping
        _vmin, _vmax = _get_vmin_vmax(self, vmin=vmin, vmax=vmax)
        arr = self.normalize(vmin=_vmin, vmax=_vmax).arr

        # Finally, need to determine whether to flip the image based on negative Affine
        # coefficients
        flip_x = self.raster_meta.transform.a < 0
        flip_y = self.raster_meta.transform.e > 0
        if flip_x:
            arr = np.flip(arr, axis=1)
        if flip_y:
            arr = np.flip(arr, axis=0)

        bounds = [[ymin, xmin], [ymax, xmax]]
        img = folium.raster_layers.ImageOverlay(
            image=arr,
            bounds=bounds,
            opacity=opacity,
            colormap=cmap_func,
            mercator_project=True,
        )

        img.add_to(m)

        # Add a colorbar legend
        if BRANCA_INSTALLED:
            cbar = _map_colorbar(colormap=cmap_func, vmin=_vmin, vmax=_vmax)
            if cbar_label:
                cbar.caption = cbar_label
            cbar.add_to(m)

        m.fit_bounds(bounds)

        return m

    def normalize(
        self, *, vmin: float | None = None, vmax: float | None = None
    ) -> Self:
        """Normalize the raster values to the range [0, 1].

        If custom vmin and vmax values are provided, values below vmin will be set to 0,
        and values above vmax will be set to 1.

        Args:
            vmin: Minimum value for normalization. Values below this will be set to 0.
                  If None, the minimum value in the array is used.
            vmax: Maximum value for normalization. Values above this will be set to 1.
                  If None, the maximum value in the array is used.
        """
        _vmin, _vmax = _get_vmin_vmax(self, vmin=vmin, vmax=vmax)

        arr = self.arr.copy()
        if _vmax > _vmin:
            arr = (arr - _vmin) / (_vmax - _vmin)
            arr = np.clip(arr, 0, 1)
        else:
            arr = np.zeros_like(arr)
        return self.__class__(arr=arr, raster_meta=self.raster_meta)

    def to_clipboard(self) -> None:
        """Copy the raster cell array to the clipboard."""
        import pandas as pd

        pd.DataFrame(self.arr).to_clipboard(index=False, header=False)

    def plot(
        self,
        *,
        ax: Axes | None = None,
        cbar_label: str | None = None,
        basemap: bool = False,
        cmap: str = "viridis",
        suppressed: Collection[float] | float = tuple(),
        **kwargs: Any,
    ) -> Axes:
        """Plot the raster on a matplotlib axis.

        Args:
            ax: A matplotlib axes object to plot on. If None, a new figure will be
                created.
            cbar_label: Label for the colorbar. If None, no label is added.
            basemap: Whether to add a basemap. Currently not implemented.
            cmap: Colormap to use for the plot.
            suppressed: Values to suppress from the plot (i.e. not display). This can be
                        useful for zeroes especially.
            **kwargs: Additional keyword arguments to pass to `rasterio.plot.show()`.
                      This includes parameters like `alpha` for transparency, and `vmin`
                      and `vmax` for controlling the color scale limits.
        """
        if not MATPLOTLIB_INSTALLED:
            msg = "The 'matplotlib' package is required for 'plot()'."
            raise ImportError(msg)

        from matplotlib import pyplot as plt
        from mpl_toolkits.axes_grid1 import make_axes_locatable

        suppressed = np.array(suppressed)

        if ax is None:
            _, _ax = plt.subplots()
            _ax: Axes
            ax = _ax

        if basemap:
            msg = "Basemap plotting is not yet implemented."
            raise NotImplementedError(msg)

        model = self.model_copy()
        model.arr = model.arr.copy()

        # Get extent of the unsuppressed values in array index coordinates
        suppressed_mask = np.isin(model.arr, suppressed)
        (x_unsuppressed,) = np.nonzero((~suppressed_mask).any(axis=0))
        (y_unsuppressed,) = np.nonzero((~suppressed_mask).any(axis=1))

        if len(x_unsuppressed) == 0 or len(y_unsuppressed) == 0:
            msg = "Raster contains no unsuppressed values; cannot plot."
            raise ValueError(msg)

        # N.B. these are array index coordinates, so np.min and np.max are safe since
        # they cannot encounter NaN values.
        min_x_unsuppressed = np.min(x_unsuppressed)
        max_x_unsuppressed = np.max(x_unsuppressed)
        min_y_unsuppressed = np.min(y_unsuppressed)
        max_y_unsuppressed = np.max(y_unsuppressed)

        # Transform to raster CRS
        x1, y1 = self.raster_meta.transform * (  # type: ignore[reportAssignmentType] overloaded tuple size in affine
            min_x_unsuppressed,
            min_y_unsuppressed,
        )
        x2, y2 = self.raster_meta.transform * (  # type: ignore[reportAssignmentType]
            max_x_unsuppressed + 1,
            max_y_unsuppressed + 1,
        )
        xmin, xmax = sorted([x1, x2])
        ymin, ymax = sorted([y1, y2])

        model.arr[suppressed_mask] = np.nan

        img, *_ = model.rio_show(ax=ax, cmap=cmap, with_bounds=True, **kwargs)

        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)

        ax.set_aspect("equal", "box")
        ax.set_yticklabels([])
        ax.set_xticklabels([])
        ax.tick_params(left=False, bottom=False, top=False, right=False)

        divider = make_axes_locatable(ax)
        cax = divider.append_axes("right", size="5%", pad=0.05)
        fig = ax.get_figure()
        if fig is not None:
            fig.colorbar(img, label=cbar_label, cax=cax)
        return ax

    def rio_show(self, **kwargs: Any) -> list[AxesImage]:
        """Plot the raster using rasterio's built-in plotting function.

        This is useful for lower-level access to rasterio's plotting capabilities.
        Generally, the `plot()` method is preferred for most use cases.

        Args:
            **kwargs: Keyword arguments to pass to `rasterio.plot.show()`. This includes
                      parameters like `alpha` for transparency, and `with_bounds` to
                      control whether to plot in spatial coordinates or array index
                      coordinates.
        """
        with self.to_rasterio_dataset() as dataset:
            return rasterio.plot.show(dataset, **kwargs).get_images()

    def as_geodataframe(self, name: str = "value") -> gpd.GeoDataFrame:
        """Create a GeoDataFrame representation of the raster."""
        import geopandas as gpd

        polygons = create_fishnet(bounds=self.bounds, res=self.raster_meta.cell_size)
        point_tuples = [polygon.centroid.coords[0] for polygon in polygons]
        raster_gdf = gpd.GeoDataFrame(
            {
                "geometry": polygons,
                name: self.sample(point_tuples, na_action="ignore"),
            },
            crs=self.raster_meta.crs,
        )

        return raster_gdf

    def gdf(self, name: str = "value") -> gpd.GeoDataFrame:
        """Create a GeoDataFrame representation of the raster.

        Alias for `as_geodataframe()`.
        """
        return self.as_geodataframe(name=name)

    def to_file(self, path: Path | str, **kwargs: Any) -> None:
        """Write the raster to a file.

        Args:
            path: Path to output file.
            **kwargs: Additional keyword arguments to pass to `rasterio.open()`. If
                      `nodata` is provided, NaN values in the raster will be replaced
                      with the nodata value.
        """
        from rastr.io_ import write_raster  # noqa: PLC0415

        return write_raster(self, path=path, **kwargs)

    def __str__(self) -> str:
        cls = self.__class__
        mean = self.mean()
        return f"{cls.__name__}(shape={self.arr.shape}, {mean=})"

    def __repr__(self) -> str:
        return str(self)

    @classmethod
    def example(cls) -> Self:
        """Create an example Raster."""
        # Peaks dataset style example
        n = 256
        x = np.linspace(-3, 3, n)
        y = np.linspace(-3, 3, n)
        x, y = np.meshgrid(x, y)
        z = np.exp(-(x**2) - y**2) * np.sin(3 * np.sqrt(x**2 + y**2))
        arr = z.astype(np.float32)

        raster_meta = RasterMeta.example()
        return cls(arr=arr, raster_meta=raster_meta)

    @classmethod
    def full_like(cls, other: Raster, *, fill_value: float) -> Self:
        """Create a raster with the same metadata as another but filled with a constant.

        Args:
            other: The raster to copy metadata from.
            fill_value: The constant value to fill all cells with.

        Returns:
            A new raster with the same shape and metadata as `other`, but with all cells
            set to `fill_value`.
        """
        arr = np.full(other.shape, fill_value, dtype=np.float32)
        return cls(arr=arr, raster_meta=other.raster_meta)

    @classmethod
    def read_file(cls, filename: Path | str, crs: CRS | str | None = None) -> Self:
        """Read raster data from a file and return an in-memory Raster object.

        Args:
            filename: Path to the raster file.
            crs: Optional coordinate reference system to override the file's CRS.
        """
        # Import here to avoid circular import (rastr.io imports Raster)
        from rastr.io_ import read_raster_inmem  # noqa: PLC0415

        return read_raster_inmem(filename, crs=crs, cls=cls)

    @classmethod
    def read_mosaic_dir(
        cls, mosaic_dir: Path | str, *, glob: str = "*.tif", crs: CRS | None = None
    ) -> Self:
        """Read a raster mosaic from a directory and return an in-memory Raster object.

        Args:
            mosaic_dir: Path to the directory containing raster files.
            glob: Glob pattern to match raster files. Default is "*.tif".
            crs: Optional coordinate reference system to override the files' CRS.

        Returns:
            A Raster object containing the mosaicked data from all matching files.
        """
        # Import here to avoid circular import (rastr.io imports Raster)
        from rastr.io_ import read_raster_mosaic_inmem  # noqa: PLC0415

        raster_obj = read_raster_mosaic_inmem(mosaic_dir, glob=glob, crs=crs)
        return cls(arr=raster_obj.arr, raster_meta=raster_obj.raster_meta)

    @overload
    def apply(
        self,
        func: Callable[Concatenate[np.ndarray, P], np.ndarray],
        *,
        raw: Literal[True],
        **kwargs: Any,
    ) -> Self: ...
    @overload
    def apply(
        self,
        func: Callable[Concatenate[float, P], float]
        | Callable[Concatenate[np.ndarray, P], np.ndarray],
        *,
        raw: Literal[False] = False,
        **kwargs: Any,
    ) -> Self: ...
    def apply(self, func: Callable, *, raw: bool = False, **kwargs: Any) -> Self:
        """Apply a function element-wise to the raster array.

        Creates a new raster instance with the same metadata (CRS, transform, etc.)
        but with the data array transformed by the provided function. The original
        raster is not modified.

        Args:
            func: The function to apply to the raster array. If `raw` is True, this
                  function should accept and return a NumPy array. If `raw` is False,
                  this function should accept and return a single float value.
            raw: If True, the function is applied directly to the entire array at
                 once. If False, the function is applied element-wise to each cell
                 in the array using `np.vectorize()`. Default is False.
            **kwargs: Additional keyword arguments to pass to the function.
        """
        new_raster = self.model_copy()
        if raw:
            new_arr = func(self.arr, **kwargs)
        else:
            new_arr = np.vectorize(func)(self.arr, **kwargs)
        new_raster.arr = np.asarray(new_arr)
        return new_raster

    def max(self) -> float:
        """Get the maximum value in the raster, ignoring NaN values.

        Returns:
            The maximum value in the raster. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanmax(self.arr))

    def min(self) -> float:
        """Get the minimum value in the raster, ignoring NaN values.

        Returns:
            The minimum value in the raster. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanmin(self.arr))

    def mean(self) -> float:
        """Get the mean value in the raster, ignoring NaN values.

        Returns:
            The mean value in the raster. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanmean(self.arr))

    def std(self) -> float:
        """Get the standard deviation of values in the raster, ignoring NaN values.

        Returns:
            The standard deviation of the raster. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanstd(self.arr))

    def quantile(self, q: float) -> float:
        """Get the specified quantile value in the raster, ignoring NaN values.

        Args:
            q: Quantile to compute, must be between 0 and 1 inclusive.

        Returns:
            The quantile value. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanquantile(self.arr, q))

    def median(self) -> float:
        """Get the median value in the raster, ignoring NaN values.

        This is equivalent to quantile(0.5).

        Returns:
            The median value in the raster. Returns NaN if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nanmedian(self.arr))

    def sum(self) -> float:
        """Get the sum of all values in the raster, ignoring NaN values.

        Returns:
            The sum of all values in the raster. Returns zero if all values are NaN.
        """
        with suppress_slice_warning():
            return float(np.nansum(self.arr))

    def unique(self) -> NDArray:
        """Get the unique cell values in the raster, including NaN.

        Returns:
            Array of unique values in the raster.
        """
        return np.unique(self.arr.flatten())

    def fillna(self, value: float) -> Self:
        """Fill NaN values in the raster with a specified value.

        See also `extrapolate()` for filling NaN values using extrapolation from data.
        """
        filled_arr = np.nan_to_num(self.arr, nan=value)
        new_raster = self.model_copy()
        new_raster.arr = filled_arr
        return new_raster

    def replace(
        self, to_replace: float | dict[float, float], value: float | None = None
    ) -> Self:
        """Replace values in the raster with other values.

        Creates a new raster with the specified values replaced. This is useful for
        operations like replacing zeros with NaNs, or vice versa.

        The method supports two interfaces:
        1. Single replacement: `raster.replace(to_replace=0, value=np.nan)`
        2. Multiple replacements using a dictionary:
           `raster.replace({0: np.nan, -999: np.nan})`

        Args:
            to_replace: Value to be replaced, or a dictionary mapping values to
                        their replacements.
            value: Replacement value. Required when to_replace is a float, must be
                   None when to_replace is a dict.

        Examples:
            >>> # Replace a single value
            >>> raster.replace(to_replace=0, value=np.nan)
            >>> # Replace multiple values
            >>> raster.replace({0: np.nan, -999: np.nan})
        """
        # Determine the replacement map
        if isinstance(to_replace, dict):
            if value is not None:
                msg = "value must be None when to_replace is a dict"
                raise ValueError(msg)
            map_ = to_replace
        else:
            if value is None:
                msg = "value must be specified when to_replace is a float"
                raise ValueError(msg)
            map_ = {to_replace: value}

        # Start with a copy of the array
        replaced_arr = self.arr.copy()

        # Check if we need to convert to float (if assigning NaN to non-float array)
        needs_float = any(
            np.isnan(new_val) for new_val in map_.values()
        ) and not np.issubdtype(replaced_arr.dtype, np.floating)
        if needs_float:
            replaced_arr = replaced_arr.astype(float)

        # Apply each replacement based on the original array values
        # to prevent chained replacements
        for old_val, new_val in map_.items():
            # Handle NaN specially since NaN != NaN
            if np.isnan(old_val):
                mask = np.isnan(self.arr)
            else:
                mask = self.arr == old_val

            replaced_arr[mask] = new_val

        new_raster = self.model_copy()
        new_raster.arr = replaced_arr
        return new_raster

    def copy(self) -> Self:  # type: ignore[override]
        """Create a copy of the raster.

        This method wraps `model_copy()` for convenience.

        Returns:
            A new Raster instance.
        """
        return self.model_copy(deep=True)

    def get_xy(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
        """Get the x and y coordinates of the raster cell centres in meshgrid format.

        Returns the coordinates of the cell centres as two separate 2D arrays in
        meshgrid format, where each array has the same shape as the raster data array.

        Returns:
            A tuple of (x, y) coordinate arrays where:
            - x: 2D array of x-coordinates of cell centres
            - y: 2D array of y-coordinates of cell centres
            Both arrays have the same shape as the raster data array.
        """
        coords = self.raster_meta.get_cell_centre_coords(self.arr.shape)
        return coords[:, :, 0], coords[:, :, 1]

    def contour(
        self, levels: Collection[float] | NDArray, *, smoothing: bool = True
    ) -> gpd.GeoDataFrame:
        """Create contour lines from the raster data, optionally with smoothing.

        The contour lines are returned as a GeoDataFrame with the contours dissolved
        by level, resulting in one row per contour level. Each row contains a
        (Multi)LineString geometry representing all contour lines for that level,
        and the contour level value in a column named 'level'.

        Consider calling `blur()` before this method to smooth the raster data before
        contouring, to denoise the contours.

        Args:
            levels: A collection or array of contour levels to generate. The contour
                    lines will be generated for each level in this sequence.
            smoothing: Defaults to true, which corresponds to applying a smoothing
                       algorithm to the contour lines. At the moment, this is the
                       Catmull-Rom spline algorithm. If set to False, the raw
                       contours will be returned without any smoothing.
        """
        import geopandas as gpd
        import skimage.measure

        all_levels = []
        all_geoms = []
        for level in levels:
            # If this is the maximum or minimum level, perturb it ever-so-slightly to
            # ensure we get contours at the edges of the raster
            perturbed_level = level
            if level == self.max():
                perturbed_level -= CONTOUR_PERTURB_EPS
            elif level == self.min():
                perturbed_level += CONTOUR_PERTURB_EPS

            contours = skimage.measure.find_contours(
                self.arr,
                level=perturbed_level,
            )

            # Construct shapely LineString objects
            # Convert to CRS from array index coordinates to raster CRS
            geoms = [
                LineString(
                    np.array(
                        rasterio.transform.xy(self.raster_meta.transform, *contour.T)
                    ).T
                )
                for contour in contours
                # Contour lines need at least three distinct points to avoid
                # degenerate geometries
                if np.unique(contour, axis=0).shape[0] > 2
            ]

            # Apply smoothing if requested
            if smoothing:
                geoms = [catmull_rom_smooth(geom) for geom in geoms]

            all_geoms.extend(geoms)
            all_levels.extend([level] * len(geoms))

        contour_gdf = gpd.GeoDataFrame(
            data={
                "level": all_levels,
            },
            geometry=all_geoms,
            crs=self.raster_meta.crs,
        )

        # Dissolve contours by level to merge all contour lines of the same level
        contour_gdf = contour_gdf.dissolve(by="level", as_index=False, sort=True)
        contour_gdf = contour_gdf.set_geometry("geometry")
        contour_gdf["geometry"] = contour_gdf["geometry"].apply(cast_multilinestring)  # pyright: ignore[reportArgumentType]
        return contour_gdf

    def sobel(self) -> Self:
        """Compute the Sobel gradient magnitude of the raster.

        This is effectively a discrete differentiation operator, computing an
        approximation of the magnitude of the gradient of the image intensity function.

        Borders are treated using half-sample symmetric sampling, i.e. repeating the
        border values. Be aware that this can lead to edge artifacts and under-estimate
        the gradient along the border pixels.

        Returns:
            New raster containing the gradient magnitude in units of raster cell units
            per unit distance (e.g. per meter).
        """
        if not self.has_square_cells:
            msg = "Sobel filter currently only supports square rasters."
            raise NotImplementedError(msg)

        from skimage import filters

        new_raster = self.model_copy()
        # Scale by cell size to convert to per unit distance
        new_raster.arr = filters.sobel(self.arr) / self.square_cell_size
        return new_raster

    def blur(self, sigma: float, *, preserve_nan: bool = True) -> Self:
        """Apply a Gaussian blur to the raster data.

        Args:
            sigma: Standard deviation for Gaussian kernel, in units of geographic
                   coordinate distance (e.g. meters). A larger sigma results in a more
                   blurred image.
            preserve_nan: If True, applies NaN-safe blurring by extrapolating NaN values
                          before blurring and restoring them afterwards. This prevents
                          NaNs from spreading into valid data during the blur operation.
        """
        from scipy.ndimage import gaussian_filter

        cell_width, cell_height = self.raster_meta.cell_size
        cell_sigma = (sigma / cell_height, sigma / cell_width)

        if preserve_nan:
            # Save the original NaN mask
            nan_mask = np.isnan(self.arr)

            # If there are no NaNs, just apply regular blur
            if not np.any(nan_mask):
                blurred_array = gaussian_filter(self.arr, sigma=cell_sigma)
            else:
                # Extrapolate to fill NaN values temporarily
                extrapolated_arr = fillna_nearest_neighbours(arr=self.arr)

                # Apply blur to the extrapolated array
                blurred_array = gaussian_filter(extrapolated_arr, sigma=cell_sigma)

                # Restore original NaN values
                blurred_array = np.where(nan_mask, np.nan, blurred_array)
        else:
            blurred_array = gaussian_filter(self.arr, sigma=cell_sigma)

        new_raster = self.model_copy()
        new_raster.arr = blurred_array
        return new_raster

    def dilate(self, radius: float) -> Self:
        """Apply a morphological dilation to the raster using a disk footprint.

        Morphological dilation sets the value of a pixel to the maximum over all pixel
        values within a local neighborhood centered about it.

        Dilation enlarges bright regions and shrinks dark regions. This is useful e.g.
        to find a region nearby a steep edge after applying a Sobel filter.

        The radius parameter controls the dilation extent. This is rounded up to the
        nearest integer multiple of the cell size, since dilation is a discrete
        operation. To reduce inaccuracies due to this rounding, consider resampling
        the raster to a smaller cell size before applying dilation.

        NaN values in the original raster are preserved in their original locations.
        They are temporarily filled during dilation to avoid spreading into valid data,
        then restored after the operation completes. The output raster maintains the
        same shape as the input.

        Args:
            radius: Radius of the disk footprint used in dilation, in units of
                    geographic coordinate distance (e.g. meters).

        Returns:
            New raster with dilated values. NaN locations are preserved from the
            original raster.
        """
        from skimage import morphology

        if not self.has_square_cells:
            msg = "Dilate currently only supports rasters with square cells."
            raise NotImplementedError(msg)

        cell_size = self.square_cell_size

        # Round up to nearest cell
        cell_radius = int(np.ceil(radius / cell_size))

        # Calculate actual radius based on rounded cell count
        radius_m = cell_radius * cell_size

        # Store original NaN mask and shape
        original_nan_mask = np.isnan(self.arr)
        original_shape = self.arr.shape

        # Handle all-NaN case
        if np.all(original_nan_mask):
            return self.model_copy()

        # Pad the raster with non-consequential values to avoid edge effects
        fill_val = self.min() - 1.0
        new_raster = self.model_copy()
        new_raster = new_raster.pad(width=radius_m, value=fill_val)

        # Replace NaNs with fill_val to avoid issues during dilation
        new_raster.arr[np.isnan(new_raster.arr)] = fill_val
        new_raster.arr = morphology.dilation(
            new_raster.arr, morphology.disk(cell_radius)
        )

        # Crop back to original size using array slicing
        new_raster.arr = new_raster.arr[
            cell_radius : cell_radius + original_shape[0],
            cell_radius : cell_radius + original_shape[1],
        ]
        # Restore original metadata
        new_raster = self.__class__(arr=new_raster.arr, meta=self.meta)

        # Restore original NaN values
        new_raster.arr[original_nan_mask] = np.nan

        return new_raster

    def extrapolate(self, method: Literal["nearest"] = "nearest") -> Self:
        """Extrapolate the raster data to fill NaN values.

        See also `fillna()` for filling NaN values with a specific value.

        If the raster is all-NaN, this method will return a copy of the raster without
        changing the NaN values.

        Args:
            method: The method to use for extrapolation. Currently only 'nearest' is
                    supported, which fills NaN values with the nearest non-NaN value.
        """
        if method not in ("nearest",):
            msg = f"Unsupported extrapolation method: {method}"
            raise NotImplementedError(msg)

        raster = self.model_copy()
        raster.arr = fillna_nearest_neighbours(arr=self.arr)

        return raster

    def pad(self, width: float, *, value: float = np.nan) -> Self:
        """Extend the raster by adding a constant fill value around the edges.

        By default, the padding value is NaN, but this can be changed via the
        `value` parameter.

        This grows the raster by adding padding around all edges. New cells are
        filled with the constant `value`.

        If the width is not an exact multiple of the cell size, the padding may be
        slightly larger than the specified width, i.e. the value is rounded up to
        the nearest whole number of cells. For non-square cells, padding is computed
        independently in x and y directions.

        Args:
            width: The width of the padding, in the same units as the raster CRS
                   (e.g. meters). This defines how far from the edge the padding
                   extends.
            value: The constant value to use for padding. Default is NaN.
        """
        cell_width, cell_height = self.raster_meta.cell_size

        # Calculate number of cells to pad in each direction
        pad_cols = int(np.ceil(width / cell_width))
        pad_rows = int(np.ceil(width / cell_height))

        # Get current bounds
        xmin, ymin, xmax, ymax = self.bounds

        # Calculate new bounds with padding
        new_xmin = xmin - (pad_cols * cell_width)
        new_ymin = ymin - (pad_rows * cell_height)
        new_xmax = xmax + (pad_cols * cell_width)
        new_ymax = ymax + (pad_rows * cell_height)

        # Create padded array
        new_height = self.arr.shape[0] + 2 * pad_rows
        new_width = self.arr.shape[1] + 2 * pad_cols

        # Create new array filled with the padding value
        padded_arr = np.full((new_height, new_width), value, dtype=self.arr.dtype)

        # Copy original array into the center of the padded array
        padded_arr[
            pad_rows : pad_rows + self.arr.shape[0],
            pad_cols : pad_cols + self.arr.shape[1],
        ] = self.arr

        # Create new transform for the padded raster
        new_transform = rasterio.transform.from_bounds(
            west=new_xmin,
            south=new_ymin,
            east=new_xmax,
            north=new_ymax,
            width=new_width,
            height=new_height,
        )

        # Create new raster metadata
        new_meta = RasterMeta(
            crs=self.raster_meta.crs,
            transform=new_transform,
        )

        return self.__class__(arr=padded_arr, raster_meta=new_meta)

    def crop(
        self,
        bounds: tuple[float, float, float, float] | Bounds | ArrayLike,
        *,
        strategy: Literal["underflow", "overflow"] = "underflow",
    ) -> Self:
        """Crop the raster to the specified bounds as (minx, miny, maxx, maxy).

        Args:
            bounds: A tuple of (minx, miny, maxx, maxy) defining the bounds to crop to.
            strategy:   The cropping strategy to use. 'underflow' will crop the raster
                        to be fully within the bounds, ignoring any cells that are
                        partially outside the bounds. 'overflow' will instead include
                        cells that intersect the bounds, ensuring the bounds area
                        remains covered with cells.

        Returns:
            A new Raster instance cropped to the specified bounds.
        """

        bounds = np.asarray(bounds)
        if len(bounds) != 4:
            msg = (
                f"bounds must be a sequence of length 4 (minx, miny, maxx, maxy); "
                f"got length {len(bounds)}"
            )
            raise ValueError(msg)

        minx, miny, maxx, maxy = bounds
        arr = self.arr

        # Get half cell sizes for cropping
        cell_width, cell_height = self.raster_meta.cell_size
        half_cell_width = cell_width / 2
        half_cell_height = cell_height / 2

        # Get the cell centre coordinates as 1D arrays
        x_coords = self.cell_x_coords
        y_coords = self.cell_y_coords

        # Get the indices to crop the array
        if strategy == "underflow":
            x_idx = (x_coords >= minx + half_cell_width) & (
                x_coords <= maxx - half_cell_width
            )
            y_idx = (y_coords >= miny + half_cell_height) & (
                y_coords <= maxy - half_cell_height
            )
        elif strategy == "overflow":
            x_idx = (x_coords > minx - half_cell_width) & (
                x_coords < maxx + half_cell_width
            )
            y_idx = (y_coords > miny - half_cell_height) & (
                y_coords < maxy + half_cell_height
            )
        else:
            msg = f"Unsupported cropping strategy: {strategy}"
            raise NotImplementedError(msg)

        # Crop the array
        cropped_arr = arr[np.ix_(y_idx, x_idx)]

        # Check the shape of the cropped array
        if cropped_arr.size == 0:
            msg = "Cropped array is empty; no cells within the specified bounds."
            raise ValueError(msg)

        # Recalculate the transform for the cropped raster
        x_coords = x_coords[x_idx]
        y_coords = y_coords[y_idx]
        transform = rasterio.transform.from_bounds(
            west=x_coords.min() - half_cell_width,
            south=y_coords.min() - half_cell_height,
            east=x_coords.max() + half_cell_width,
            north=y_coords.max() + half_cell_height,
            width=cropped_arr.shape[1],
            height=cropped_arr.shape[0],
        )

        # Update the raster
        cls = self.__class__
        new_meta = RasterMeta(crs=self.raster_meta.crs, transform=transform)
        return cls(arr=cropped_arr, raster_meta=new_meta)

    def to_bounds(
        self,
        bounds: tuple[float, float, float, float] | Bounds | ArrayLike,
        *,
        strategy: Literal["underflow", "overflow"] = "underflow",
    ) -> Self:
        """Resize the raster to the specified bounds, cropping or padding as needed.

        This method behaves like `crop()` when the bounds are within the current raster,
        but also allows expanding the raster by padding with NaN-valued cells when the
        bounds extend beyond the current raster boundaries.

        Args:
            bounds: A tuple of (minx, miny, maxx, maxy) defining the target bounds.
            strategy:
                The strategy to use when determining cell inclusion at boundaries.
                'underflow' will only include cells fully within the bounds,
                'overflow' will include cells that intersect the bounds.

        Returns:
            A new Raster instance resized to the specified bounds.
        """
        bounds = np.asarray(bounds)
        if len(bounds) != 4:
            msg = (
                f"bounds must be a sequence of length 4 (minx, miny, maxx, maxy); "
                f"got length {len(bounds)}"
            )
            raise ValueError(msg)

        if strategy not in ("underflow", "overflow"):
            msg = f"Unsupported strategy: {strategy}"
            raise NotImplementedError(msg)

        target_minx, target_miny, target_maxx, target_maxy = bounds
        cell_width, cell_height = self.raster_meta.cell_size
        half_cell_width = cell_width / 2
        half_cell_height = cell_height / 2

        sign = 1 if strategy == "underflow" else -1
        x_range = (
            target_minx + sign * half_cell_width,
            target_maxx - sign * half_cell_width,
        )
        y_range = (
            target_miny + sign * half_cell_height,
            target_maxy - sign * half_cell_height,
        )

        if x_range[1] < x_range[0] or y_range[1] < y_range[0]:
            msg = "No cells within the specified bounds."
            raise ValueError(msg)

        # Get current cell coordinates
        current_x_coords = self.cell_x_coords
        current_y_coords = self.cell_y_coords

        # Determine the grid alignment based on the first cell center
        first_x = current_x_coords[0]
        first_y = current_y_coords[0]

        # Calculate indices for target grid (aligned to original grid)
        # x-direction is always ascending
        x_start = np.ceil((x_range[0] - first_x) / cell_width)
        x_end = np.floor((x_range[1] - first_x) / cell_width)

        # y-direction may be ascending or descending depending on transform
        # Determine direction from current_y_coords
        y_ascending = (
            len(current_y_coords) > 1 and current_y_coords[1] > current_y_coords[0]
        )

        if y_ascending:
            y_start = np.ceil((y_range[0] - first_y) / cell_height)
            y_end = np.floor((y_range[1] - first_y) / cell_height)
            y_indices = np.arange(y_start, y_end + 1, dtype=int)
            target_y_coords = first_y + y_indices * cell_height
        else:
            y_start = np.ceil((first_y - y_range[1]) / cell_height)
            y_end = np.floor((first_y - y_range[0]) / cell_height)
            y_indices = np.arange(y_start, y_end + 1, dtype=int)
            target_y_coords = first_y - y_indices * cell_height

        # Generate target coordinates
        x_indices = np.arange(x_start, x_end + 1, dtype=int)

        if len(x_indices) == 0 or len(y_indices) == 0:
            msg = "No cells within the specified bounds."
            raise ValueError(msg)

        target_x_coords = first_x + x_indices * cell_width

        output_shape = (len(target_y_coords), len(target_x_coords))
        output_arr = np.full(output_shape, np.nan, dtype=float)

        # For each target coordinate, find matching indices in current raster
        # Use broadcasting to create boolean masks
        x_match_mask = (
            np.abs(target_x_coords[:, np.newaxis] - current_x_coords[np.newaxis, :])
            < COORD_MATCH_TOLERANCE
        )
        y_match_mask = (
            np.abs(target_y_coords[:, np.newaxis] - current_y_coords[np.newaxis, :])
            < COORD_MATCH_TOLERANCE
        )

        y_target_indices, y_current_indices = np.where(y_match_mask)
        x_target_indices, x_current_indices = np.where(x_match_mask)

        for ty_idx, cy_idx in zip(y_target_indices, y_current_indices, strict=False):
            for tx_idx, cx_idx in zip(
                x_target_indices, x_current_indices, strict=False
            ):
                output_arr[ty_idx, tx_idx] = self.arr[cy_idx, cx_idx]

        transform = rasterio.transform.from_bounds(
            west=target_x_coords.min() - half_cell_width,
            south=target_y_coords.min() - half_cell_height,
            east=target_x_coords.max() + half_cell_width,
            north=target_y_coords.max() + half_cell_height,
            width=output_arr.shape[1],
            height=output_arr.shape[0],
        )

        # Create the new raster
        cls = self.__class__
        new_meta = RasterMeta(crs=self.raster_meta.crs, transform=transform)
        return cls(arr=output_arr, raster_meta=new_meta)

    def taper_border(self, width: float, *, limit: float = 0.0) -> Self:
        """Taper values to a limiting value around the border of the raster.

        By default, the borders are tapered to zero, but this can be changed via the
        `limit` parameter.

        This keeps the raster size the same, overwriting values in the border area.
        To instead grow the raster, consider using `pad()` followed by `taper_border()`.

        The tapering is linear from the cell centres around the border of the raster,
        so the value at the edge of the raster will be equal to `limit`.

        Args:
            width: The width of the taper, in the same units as the raster CRS
                   (e.g. meters). This defines how far from the edge the tapering
                   starts.
            limit: The limiting value to taper to at the edges. Default is zero.
        """

        cell_width, cell_height = self.raster_meta.cell_size
        width_in_cols = width / cell_width
        width_in_rows = width / cell_height

        # Calculate the distance from the edge in cell units
        arr_height, arr_width = self.arr.shape
        y_indices, x_indices = np.indices((int(arr_height), int(arr_width)))
        dist_from_left = x_indices
        dist_from_right = arr_width - 1 - x_indices
        dist_from_top = y_indices
        dist_from_bottom = arr_height - 1 - y_indices
        dist_from_edge_cols = np.minimum(dist_from_left, dist_from_right)
        dist_from_edge_rows = np.minimum(dist_from_top, dist_from_bottom)

        # Distance from cell centres to nearest border edge in CRS units
        dist_from_edge = np.minimum(
            dist_from_edge_cols * cell_width,
            dist_from_edge_rows * cell_height,
        )

        # Maintain edge-band behavior by rounding up separately per axis
        within_x_band = dist_from_edge_cols < np.ceil(width_in_cols)
        within_y_band = dist_from_edge_rows < np.ceil(width_in_rows)
        mask = within_x_band | within_y_band
        masked_dist_arr = np.where(mask, dist_from_edge, np.nan)
        masked_arr = np.where(mask, self.arr, np.nan)

        # Calculate the tapering factor based on the distance from the edge
        taper_factor = np.clip(masked_dist_arr / width, 0.0, 1.0)
        tapered_values = limit + (masked_arr - limit) * taper_factor

        # Create the new raster array
        new_arr = self.arr.copy()
        new_arr[mask] = tapered_values[mask]
        new_raster = self.model_copy()
        new_raster.arr = new_arr

        return new_raster

    def clip(
        self,
        polygon: BaseGeometry,
        *,
        strategy: Literal["centres"] = "centres",
    ) -> Self:
        """Clip the raster to the specified polygon, replacing cells outside with NaN.

        The clipping strategy determines how to handle cells that are partially
        within the polygon. Currently, only the 'centres' strategy is supported, which
        retains cells whose centres fall within the polygon.

        Args:
            polygon: A shapely Polygon or MultiPolygon defining the area to clip to.
            strategy: The clipping strategy to use. Currently only 'centres' is
                      supported, which retains cells whose centres fall within the
                      polygon.

        Returns:
            A new Raster with cells outside the polygon set to NaN.

        Raises:
            TypeError: If the provided geometry is not a Polygon or MultiPolygon.
        """
        if not isinstance(polygon, Polygon | MultiPolygon):
            msg = (
                f"Only Polygon and MultiPolygon geometries are supported for clipping, "
                f"got {type(polygon).__name__}"
            )
            raise TypeError(msg)

        if strategy != "centres":
            msg = f"Unsupported clipping strategy: {strategy}"
            raise NotImplementedError(msg)

        mask_raster = self._polygon_indicator(polygon)

        raster = self.model_copy()
        raster.arr = np.where(mask_raster.arr, raster.arr, np.nan)

        return raster

    def _trim_value(self, *, value_mask: NDArray[np.bool_], value_name: str) -> Self:
        """Crop the raster by trimming away slices matching the mask at the edges.

        Args:
            value_mask: Boolean mask where True indicates values to trim
            value_name: Name of the value type for error messages (e.g., 'NaN', 'zero')
        """
        arr = self.arr

        # Check if the entire array matches the mask
        if np.all(value_mask):
            msg = f"Cannot crop raster: all values are {value_name}"
            raise ValueError(msg)

        # Find rows and columns that are not all matching the mask
        row_mask = np.all(value_mask, axis=1)
        col_mask = np.all(value_mask, axis=0)

        # Find the bounding indices
        (row_indices,) = np.where(~row_mask)
        (col_indices,) = np.where(~col_mask)

        min_row, max_row = row_indices[0], row_indices[-1]
        min_col, max_col = col_indices[0], col_indices[-1]

        # Crop the array
        cropped_arr = arr[min_row : max_row + 1, min_col : max_col + 1]

        # Shift the transform by the number of pixels cropped (min_col, min_row)
        new_transform = (
            self.raster_meta.transform
            * rasterio.transform.Affine.translation(min_col, min_row)
        )

        # Create new metadata
        new_meta = RasterMeta(
            crs=self.raster_meta.crs,
            transform=new_transform,
        )

        return self.__class__(arr=cropped_arr, raster_meta=new_meta)

    def trim_nan(self) -> Self:
        """Crop the raster by trimming away all-NaN slices at the edges.

        This effectively trims the raster to the smallest bounding box that contains all
        of the non-NaN values. Note that this does not guarantee no NaN values at all
        around the edges, only that there won't be entire edges which are all-NaN.

        Consider using `.extrapolate()` for further cleanup of NaN values.
        """
        return self._trim_value(value_mask=np.isnan(self.arr), value_name="NaN")

    def trim_zeros(self) -> Self:
        """Crop the raster by trimming away all-zero slices at the edges.

        This effectively trims the raster to the smallest bounding box that contains all
        of the non-zero values. Note that this does not guarantee no zero values at all
        around the edges, only that there won't be entire edges which are all-zero.
        """
        return self._trim_value(value_mask=(self.arr == 0), value_name="zero")

    def resample(
        self,
        cell_size: tuple[float, float] | float,
        *,
        method: Literal["bilinear"] = "bilinear",
    ) -> Self:
        """Resample the raster data to a new resolution.

        If the new cell size is not an exact multiple of the current cell size, the
        overall raster bounds may increase slightly. The affine transform will keep
        the same shift, i.e. the top-left corner of the raster will remain in the same'
        coordinate location. A corollary is that the overall centre of the raster bounds
        will not necessary be the same as the original raster.

        Args:
            cell_size: The desired cell size for the resampled raster. This can be a
            single float for square cells, or a tuple of (cell_width, cell_height) for
            rectangular cells.
            method: The resampling method to use. Only 'bilinear' is supported.
        """
        if method not in ("bilinear",):
            msg = f"Unsupported resampling method: {method}"
            raise NotImplementedError(msg)

        target_cell_width, target_cell_height = _ensure_pair(cell_size)
        source_cell_width, source_cell_height = self.raster_meta.cell_size

        x_factor = source_cell_width / target_cell_width
        y_factor = source_cell_height / target_cell_height

        cls = self.__class__
        # Use the rasterio dataset with proper context management
        with self.to_rasterio_dataset() as dataset:
            # N.B. the new height and width may increase slightly.
            new_height = int(np.ceil(dataset.height * y_factor))
            new_width = int(np.ceil(dataset.width * x_factor))

            # Resample via rasterio
            (new_arr,) = dataset.read(  # Assume exactly one band
                out_shape=(dataset.count, new_height, new_width),
                resampling=Resampling.bilinear,
            )

            # Create new RasterMeta with the exact requested cell size
            new_raster_meta = RasterMeta(
                transform=Affine(
                    target_cell_width,
                    dataset.transform.b,
                    dataset.transform.c,
                    dataset.transform.d,
                    -target_cell_height,
                    dataset.transform.f,
                ),
                crs=self.raster_meta.crs,
            )

            return cls(arr=new_arr, raster_meta=new_raster_meta)

    def replace_polygon(
        self,
        polygon: BaseGeometry | dict[BaseGeometry, float],
        value: float | None = None,
    ) -> Self:
        """Replace values within the specified polygon(s) with other values.

        Creates a new raster with the specified values replaced. This is useful for
        operations like masking or setting regions to NaN.

        The method supports two interfaces:
        1. Single replacement: `raster.replace_polygon(polygon1, value=np.nan)`
        2. Multiple replacements using a dictionary:
           `raster.replace_polygon({polygon1: 0, polygon2: 1})`

        Args:
            polygon: Geometry to replace, or dict mapping geometries to values.
            value: Replacement value. Required if polygon is a geometry, None if polygon
                   is a dict.

        Examples:
            >>> # Replace a single polygon
            >>> raster.replace_polygon(polygon1, value=np.nan)
            >>> # Replace multiple polygons
            >>> raster.replace_polygon({polygon1: 0, polygon2: 1})
        """
        # Normalize input to dict format
        if isinstance(polygon, dict):
            if value is not None:
                msg = "value must be None when polygon is a dict"
                raise ValueError(msg)
            replacements = polygon
        else:
            if value is None:
                msg = "value must be specified when polygon is a geometry"
                raise ValueError(msg)
            replacements = {polygon: value}

        # Validate all geometries upfront
        for geom in replacements.keys():
            if not isinstance(geom, Polygon | MultiPolygon):
                msg = (
                    f"Only Polygon and MultiPolygon geometries are supported, "
                    f"got {type(geom).__name__}"
                )
                raise TypeError(msg)

        # Create copy and convert to float if needed
        raster = self.model_copy()
        needs_float = any(
            isinstance(v, float) or (v is not None and np.isnan(v))
            for v in replacements.values()
        )
        if needs_float and not np.issubdtype(raster.arr.dtype, np.floating):
            raster.arr = raster.arr.astype(float)

        # Apply all replacements
        for geom, val in replacements.items():
            mask_raster = self._polygon_indicator(geom)
            raster.arr = np.where(mask_raster.arr, val, raster.arr)

        return raster

    def _polygon_indicator(self, geom: BaseGeometry) -> Self:
        """Create a binary mask with 1 inside the polygon, 0 outside.

        The mask has the same shape and transform as the raster.

        Args:
            geom:
                A shapely geometry (typically Polygon or MultiPolygon)
                to rasterize.

        Returns:
            Raster:
                A new Raster where cells inside the geometry are 1,
                and outside are 0.
        """

        arr = rasterio.features.rasterize(
            [(geom, 1)],
            fill=0,
            out_shape=self.shape,
            transform=self.meta.transform,
            dtype=np.uint8,
        )

        # Create a new raster with the binary mask array
        raster = self.model_copy()
        raster.arr = arr

        return raster

bbox property

Bounding box of the raster as a shapely polygon.

bounds property

Bounding box of the raster as a named tuple with xmin, ymin, xmax, ymax.

The bounds represent the outer edges of the raster cells.

cell_centre_coords property

Get the coordinates of the cell centres in the raster.

cell_size property

Convenience property to access the cell size via meta.

cell_x_coords property

Get the x coordinates of the cell centres in the raster.

cell_y_coords property

Get the y coordinates of the cell centres in the raster.

crs property writable

Convenience property to access the CRS via meta.

has_square_cells property

Check if the raster has square cells.

meta property writable

Alias for raster_meta.

origin property writable

The grid origin (x, y) — the corner of the first pixel.

This is the translation component (c, f) of the affine transform. For north-up rasters this corresponds to (xmin, ymax); for south-up rasters it corresponds to (xmin, ymin).

shape property

Shape of the raster array.

square_cell_size property

Convenience property to access the square cell size via meta.

transform property writable

Convenience property to access the transform via meta.

__array_ufunc__(ufunc, method, *inputs, **kwargs)

Support NumPy ufuncs between ndarrays and Rasters.

Carries out the ufunc on the underlying array when the ndarray operand has the same shape as this Raster. Returns NotImplemented for unsupported ufuncs, unsupported methods, or shape mismatches.

Source code in src/rastr/raster.py

def __array_ufunc__(
    self, ufunc: np.ufunc, method: str, *inputs: Any, **kwargs: Any
) -> Self:
    """Support NumPy ufuncs between ndarrays and Rasters.

    Carries out the ufunc on the underlying array when the ndarray operand
    has the same shape as this Raster.  Returns ``NotImplemented`` for
    unsupported ufuncs, unsupported methods, or shape mismatches.
    """

    # There are over 60 ufuncs and not all output 2D arrays. For example, reduction
    # ops like np.sum and np.min return scalars while some ufuncs can supposedly
    # return tuples (np.modf). Some would also create types with are currently
    # unsupported by rastr, like np.greater creating bool. Therefore, we choose
    # to whitelist supported ufuncs, and aim to gradually expand this list. This
    # approach allows us to fail fast with a clear error message when an unsupported
    # ufunc is used, rather than silently passing and potentially producing
    # unexpected results.
    if method != "__call__" or ufunc not in _RASTER_SUPPORTED_UFUNCS:
        return NotImplemented
    new_inputs = []
    for inp in inputs:
        if isinstance(inp, np.ndarray):
            if inp.shape != self.shape:
                msg = (
                    f"Cannot apply ufunc '{ufunc.__name__}' between ndarray of"
                    f" shape {inp.shape} and Raster of shape {self.shape}:"
                    " shapes must be equal."
                )
                raise ValueError(msg)
            new_inputs.append(inp)
        elif isinstance(inp, Raster):
            new_inputs.append(inp.arr)
        else:
            new_inputs.append(inp)
    cls = self.__class__
    result_arr = ufunc(*new_inputs, **kwargs)
    return cls(arr=result_arr, raster_meta=self.raster_meta)

__eq__(other)

Check equality of two Raster objects.

Source code in src/rastr/raster.py

def __eq__(self, other: object) -> bool:
    """Check equality of two Raster objects."""
    if not isinstance(other, Raster):
        return NotImplemented
    return (
        np.array_equal(self.arr, other.arr, equal_nan=True)
        and self.raster_meta == other.raster_meta
    )

abs()

Compute the absolute value of the raster.

Returns a new raster with the absolute value of each cell. The original raster is not modified.

Returns:

Type	Description
Self	A new Raster instance with the absolute values.

Source code in src/rastr/raster.py

def abs(self) -> Self:
    """Compute the absolute value of the raster.

    Returns a new raster with the absolute value of each cell. The original raster
    is not modified.

    Returns:
        A new Raster instance with the absolute values.
    """
    cls = self.__class__
    return cls(arr=np.abs(self.arr), raster_meta=self.raster_meta)

apply(func, *, raw=False, **kwargs)

apply(func: Callable[Concatenate[np.ndarray, P], np.ndarray], *, raw: Literal[True], **kwargs: Any) -> Self

apply(func: Callable[Concatenate[float, P], float] | Callable[Concatenate[np.ndarray, P], np.ndarray], *, raw: Literal[False] = False, **kwargs: Any) -> Self

Apply a function element-wise to the raster array.

Creates a new raster instance with the same metadata (CRS, transform, etc.) but with the data array transformed by the provided function. The original raster is not modified.

Parameters:

Name	Type	Description	Default
func	Callable	The function to apply to the raster array. If raw is True, this function should accept and return a NumPy array. If raw is False, this function should accept and return a single float value.	required
raw	bool	If True, the function is applied directly to the entire array at once. If False, the function is applied element-wise to each cell in the array using np.vectorize(). Default is False.	False
**kwargs	Any	Additional keyword arguments to pass to the function.	{}

Source code in src/rastr/raster.py

def apply(self, func: Callable, *, raw: bool = False, **kwargs: Any) -> Self:
    """Apply a function element-wise to the raster array.

    Creates a new raster instance with the same metadata (CRS, transform, etc.)
    but with the data array transformed by the provided function. The original
    raster is not modified.

    Args:
        func: The function to apply to the raster array. If `raw` is True, this
              function should accept and return a NumPy array. If `raw` is False,
              this function should accept and return a single float value.
        raw: If True, the function is applied directly to the entire array at
             once. If False, the function is applied element-wise to each cell
             in the array using `np.vectorize()`. Default is False.
        **kwargs: Additional keyword arguments to pass to the function.
    """
    new_raster = self.model_copy()
    if raw:
        new_arr = func(self.arr, **kwargs)
    else:
        new_arr = np.vectorize(func)(self.arr, **kwargs)
    new_raster.arr = np.asarray(new_arr)
    return new_raster

as_geodataframe(name='value')

Create a GeoDataFrame representation of the raster.

Source code in src/rastr/raster.py

def as_geodataframe(self, name: str = "value") -> gpd.GeoDataFrame:
    """Create a GeoDataFrame representation of the raster."""
    import geopandas as gpd

    polygons = create_fishnet(bounds=self.bounds, res=self.raster_meta.cell_size)
    point_tuples = [polygon.centroid.coords[0] for polygon in polygons]
    raster_gdf = gpd.GeoDataFrame(
        {
            "geometry": polygons,
            name: self.sample(point_tuples, na_action="ignore"),
        },
        crs=self.raster_meta.crs,
    )

    return raster_gdf

astype(dtype)

Cast the raster array to a specified dtype.

Returns a new raster with the array cast to the given dtype. The original raster is not modified.

Parameters:

Name	Type	Description	Default
dtype	DTypeLike	Target data type (e.g. "float32", np.int16).	required

Returns:

Type	Description
Self	A new Raster instance with the array cast to the specified dtype.

Source code in src/rastr/raster.py

def astype(self, dtype: DTypeLike) -> Self:
    """Cast the raster array to a specified dtype.

    Returns a new raster with the array cast to the given dtype. The original
    raster is not modified.

    Args:
        dtype: Target data type (e.g. ``"float32"``, ``np.int16``).

    Returns:
        A new Raster instance with the array cast to the specified dtype.
    """
    cls = self.__class__
    return cls(arr=self.arr.astype(dtype), raster_meta=self.raster_meta)

blur(sigma, *, preserve_nan=True)

Apply a Gaussian blur to the raster data.

Parameters:

Name	Type	Description	Default
sigma	float	Standard deviation for Gaussian kernel, in units of geographic coordinate distance (e.g. meters). A larger sigma results in a more blurred image.	required
preserve_nan	bool	If True, applies NaN-safe blurring by extrapolating NaN values before blurring and restoring them afterwards. This prevents NaNs from spreading into valid data during the blur operation.	True

Source code in src/rastr/raster.py

def blur(self, sigma: float, *, preserve_nan: bool = True) -> Self:
    """Apply a Gaussian blur to the raster data.

    Args:
        sigma: Standard deviation for Gaussian kernel, in units of geographic
               coordinate distance (e.g. meters). A larger sigma results in a more
               blurred image.
        preserve_nan: If True, applies NaN-safe blurring by extrapolating NaN values
                      before blurring and restoring them afterwards. This prevents
                      NaNs from spreading into valid data during the blur operation.
    """
    from scipy.ndimage import gaussian_filter

    cell_width, cell_height = self.raster_meta.cell_size
    cell_sigma = (sigma / cell_height, sigma / cell_width)

    if preserve_nan:
        # Save the original NaN mask
        nan_mask = np.isnan(self.arr)

        # If there are no NaNs, just apply regular blur
        if not np.any(nan_mask):
            blurred_array = gaussian_filter(self.arr, sigma=cell_sigma)
        else:
            # Extrapolate to fill NaN values temporarily
            extrapolated_arr = fillna_nearest_neighbours(arr=self.arr)

            # Apply blur to the extrapolated array
            blurred_array = gaussian_filter(extrapolated_arr, sigma=cell_sigma)

            # Restore original NaN values
            blurred_array = np.where(nan_mask, np.nan, blurred_array)
    else:
        blurred_array = gaussian_filter(self.arr, sigma=cell_sigma)

    new_raster = self.model_copy()
    new_raster.arr = blurred_array
    return new_raster

check_2d_array(v) classmethod

Validator to ensure the cell array is 2D.

Source code in src/rastr/raster.py

@field_validator("arr")
@classmethod
def check_2d_array(cls, v: NDArray) -> NDArray:
    """Validator to ensure the cell array is 2D."""
    if v.ndim != 2:
        msg = "Cell array must be 2D"
        raise RasterCellArrayShapeError(msg)
    return v

clamp(a_min=None, a_max=None)

Clip (clamp) values to a specified range.

Returns a new raster with values clipped to [a_min, a_max].

Parameters:

Name	Type	Description	Default
a_min	float | None	Minimum value. Values below this will be set to a_min. If None, no minimum clipping is applied.	None
a_max	float | None	Maximum value. Values above this will be set to a_max. If None, no maximum clipping is applied.	None

Returns:

Type	Description
Self	A new Raster instance with clipped values.

Source code in src/rastr/raster.py

def clamp(
    self,
    a_min: float | None = None,
    a_max: float | None = None,
) -> Self:
    """Clip (clamp) values to a specified range.

    Returns a new raster with values clipped to [a_min, a_max].

    Args:
        a_min: Minimum value. Values below this will be set to a_min.
               If None, no minimum clipping is applied.
        a_max: Maximum value. Values above this will be set to a_max.
               If None, no maximum clipping is applied.

    Returns:
        A new Raster instance with clipped values.
    """
    cls = self.__class__
    return cls(
        arr=np.clip(self.arr, a_min, a_max),
        raster_meta=self.raster_meta,
    )

clip(polygon, *, strategy='centres')

Clip the raster to the specified polygon, replacing cells outside with NaN.

The clipping strategy determines how to handle cells that are partially within the polygon. Currently, only the 'centres' strategy is supported, which retains cells whose centres fall within the polygon.

Parameters:

Name	Type	Description	Default
polygon	BaseGeometry	A shapely Polygon or MultiPolygon defining the area to clip to.	required
strategy	Literal['centres']	The clipping strategy to use. Currently only 'centres' is supported, which retains cells whose centres fall within the polygon.	'centres'

Returns:

Type	Description
Self	A new Raster with cells outside the polygon set to NaN.

Raises:

Type	Description
TypeError	If the provided geometry is not a Polygon or MultiPolygon.

Source code in src/rastr/raster.py

def clip(
    self,
    polygon: BaseGeometry,
    *,
    strategy: Literal["centres"] = "centres",
) -> Self:
    """Clip the raster to the specified polygon, replacing cells outside with NaN.

    The clipping strategy determines how to handle cells that are partially
    within the polygon. Currently, only the 'centres' strategy is supported, which
    retains cells whose centres fall within the polygon.

    Args:
        polygon: A shapely Polygon or MultiPolygon defining the area to clip to.
        strategy: The clipping strategy to use. Currently only 'centres' is
                  supported, which retains cells whose centres fall within the
                  polygon.

    Returns:
        A new Raster with cells outside the polygon set to NaN.

    Raises:
        TypeError: If the provided geometry is not a Polygon or MultiPolygon.
    """
    if not isinstance(polygon, Polygon | MultiPolygon):
        msg = (
            f"Only Polygon and MultiPolygon geometries are supported for clipping, "
            f"got {type(polygon).__name__}"
        )
        raise TypeError(msg)

    if strategy != "centres":
        msg = f"Unsupported clipping strategy: {strategy}"
        raise NotImplementedError(msg)

    mask_raster = self._polygon_indicator(polygon)

    raster = self.model_copy()
    raster.arr = np.where(mask_raster.arr, raster.arr, np.nan)

    return raster

contour(levels, *, smoothing=True)

Create contour lines from the raster data, optionally with smoothing.

The contour lines are returned as a GeoDataFrame with the contours dissolved by level, resulting in one row per contour level. Each row contains a (Multi)LineString geometry representing all contour lines for that level, and the contour level value in a column named 'level'.

Consider calling blur() before this method to smooth the raster data before contouring, to denoise the contours.

Parameters:

Name	Type	Description	Default
levels	Collection[float] | NDArray	A collection or array of contour levels to generate. The contour lines will be generated for each level in this sequence.	required
smoothing	bool	Defaults to true, which corresponds to applying a smoothing algorithm to the contour lines. At the moment, this is the Catmull-Rom spline algorithm. If set to False, the raw contours will be returned without any smoothing.	True

Source code in src/rastr/raster.py

def contour(
    self, levels: Collection[float] | NDArray, *, smoothing: bool = True
) -> gpd.GeoDataFrame:
    """Create contour lines from the raster data, optionally with smoothing.

    The contour lines are returned as a GeoDataFrame with the contours dissolved
    by level, resulting in one row per contour level. Each row contains a
    (Multi)LineString geometry representing all contour lines for that level,
    and the contour level value in a column named 'level'.

    Consider calling `blur()` before this method to smooth the raster data before
    contouring, to denoise the contours.

    Args:
        levels: A collection or array of contour levels to generate. The contour
                lines will be generated for each level in this sequence.
        smoothing: Defaults to true, which corresponds to applying a smoothing
                   algorithm to the contour lines. At the moment, this is the
                   Catmull-Rom spline algorithm. If set to False, the raw
                   contours will be returned without any smoothing.
    """
    import geopandas as gpd
    import skimage.measure

    all_levels = []
    all_geoms = []
    for level in levels:
        # If this is the maximum or minimum level, perturb it ever-so-slightly to
        # ensure we get contours at the edges of the raster
        perturbed_level = level
        if level == self.max():
            perturbed_level -= CONTOUR_PERTURB_EPS
        elif level == self.min():
            perturbed_level += CONTOUR_PERTURB_EPS

        contours = skimage.measure.find_contours(
            self.arr,
            level=perturbed_level,
        )

        # Construct shapely LineString objects
        # Convert to CRS from array index coordinates to raster CRS
        geoms = [
            LineString(
                np.array(
                    rasterio.transform.xy(self.raster_meta.transform, *contour.T)
                ).T
            )
            for contour in contours
            # Contour lines need at least three distinct points to avoid
            # degenerate geometries
            if np.unique(contour, axis=0).shape[0] > 2
        ]

        # Apply smoothing if requested
        if smoothing:
            geoms = [catmull_rom_smooth(geom) for geom in geoms]

        all_geoms.extend(geoms)
        all_levels.extend([level] * len(geoms))

    contour_gdf = gpd.GeoDataFrame(
        data={
            "level": all_levels,
        },
        geometry=all_geoms,
        crs=self.raster_meta.crs,
    )

    # Dissolve contours by level to merge all contour lines of the same level
    contour_gdf = contour_gdf.dissolve(by="level", as_index=False, sort=True)
    contour_gdf = contour_gdf.set_geometry("geometry")
    contour_gdf["geometry"] = contour_gdf["geometry"].apply(cast_multilinestring)  # pyright: ignore[reportArgumentType]
    return contour_gdf

copy()

Create a copy of the raster.

This method wraps model_copy() for convenience.

Returns:

Type	Description
Self	A new Raster instance.

Source code in src/rastr/raster.py

def copy(self) -> Self:  # type: ignore[override]
    """Create a copy of the raster.

    This method wraps `model_copy()` for convenience.

    Returns:
        A new Raster instance.
    """
    return self.model_copy(deep=True)

crop(bounds, *, strategy='underflow')

Crop the raster to the specified bounds as (minx, miny, maxx, maxy).

Parameters:

Name	Type	Description	Default
bounds	tuple[float, float, float, float] | Bounds | ArrayLike	A tuple of (minx, miny, maxx, maxy) defining the bounds to crop to.	required
strategy	Literal['underflow', 'overflow']	The cropping strategy to use. 'underflow' will crop the raster to be fully within the bounds, ignoring any cells that are partially outside the bounds. 'overflow' will instead include cells that intersect the bounds, ensuring the bounds area remains covered with cells.	'underflow'

Returns:

Type	Description
Self	A new Raster instance cropped to the specified bounds.

Source code in src/rastr/raster.py

def crop(
    self,
    bounds: tuple[float, float, float, float] | Bounds | ArrayLike,
    *,
    strategy: Literal["underflow", "overflow"] = "underflow",
) -> Self:
    """Crop the raster to the specified bounds as (minx, miny, maxx, maxy).

    Args:
        bounds: A tuple of (minx, miny, maxx, maxy) defining the bounds to crop to.
        strategy:   The cropping strategy to use. 'underflow' will crop the raster
                    to be fully within the bounds, ignoring any cells that are
                    partially outside the bounds. 'overflow' will instead include
                    cells that intersect the bounds, ensuring the bounds area
                    remains covered with cells.

    Returns:
        A new Raster instance cropped to the specified bounds.
    """

    bounds = np.asarray(bounds)
    if len(bounds) != 4:
        msg = (
            f"bounds must be a sequence of length 4 (minx, miny, maxx, maxy); "
            f"got length {len(bounds)}"
        )
        raise ValueError(msg)

    minx, miny, maxx, maxy = bounds
    arr = self.arr

    # Get half cell sizes for cropping
    cell_width, cell_height = self.raster_meta.cell_size
    half_cell_width = cell_width / 2
    half_cell_height = cell_height / 2

    # Get the cell centre coordinates as 1D arrays
    x_coords = self.cell_x_coords
    y_coords = self.cell_y_coords

    # Get the indices to crop the array
    if strategy == "underflow":
        x_idx = (x_coords >= minx + half_cell_width) & (
            x_coords <= maxx - half_cell_width
        )
        y_idx = (y_coords >= miny + half_cell_height) & (
            y_coords <= maxy - half_cell_height
        )
    elif strategy == "overflow":
        x_idx = (x_coords > minx - half_cell_width) & (
            x_coords < maxx + half_cell_width
        )
        y_idx = (y_coords > miny - half_cell_height) & (
            y_coords < maxy + half_cell_height
        )
    else:
        msg = f"Unsupported cropping strategy: {strategy}"
        raise NotImplementedError(msg)

    # Crop the array
    cropped_arr = arr[np.ix_(y_idx, x_idx)]

    # Check the shape of the cropped array
    if cropped_arr.size == 0:
        msg = "Cropped array is empty; no cells within the specified bounds."
        raise ValueError(msg)

    # Recalculate the transform for the cropped raster
    x_coords = x_coords[x_idx]
    y_coords = y_coords[y_idx]
    transform = rasterio.transform.from_bounds(
        west=x_coords.min() - half_cell_width,
        south=y_coords.min() - half_cell_height,
        east=x_coords.max() + half_cell_width,
        north=y_coords.max() + half_cell_height,
        width=cropped_arr.shape[1],
        height=cropped_arr.shape[0],
    )

    # Update the raster
    cls = self.__class__
    new_meta = RasterMeta(crs=self.raster_meta.crs, transform=transform)
    return cls(arr=cropped_arr, raster_meta=new_meta)

dilate(radius)

Apply a morphological dilation to the raster using a disk footprint.

Morphological dilation sets the value of a pixel to the maximum over all pixel values within a local neighborhood centered about it.

Dilation enlarges bright regions and shrinks dark regions. This is useful e.g. to find a region nearby a steep edge after applying a Sobel filter.

The radius parameter controls the dilation extent. This is rounded up to the nearest integer multiple of the cell size, since dilation is a discrete operation. To reduce inaccuracies due to this rounding, consider resampling the raster to a smaller cell size before applying dilation.

NaN values in the original raster are preserved in their original locations. They are temporarily filled during dilation to avoid spreading into valid data, then restored after the operation completes. The output raster maintains the same shape as the input.

Parameters:

Name	Type	Description	Default
radius	float	Radius of the disk footprint used in dilation, in units of geographic coordinate distance (e.g. meters).	required

Returns:

Type	Description
Self	New raster with dilated values. NaN locations are preserved from the
Self	original raster.

Source code in src/rastr/raster.py

def dilate(self, radius: float) -> Self:
    """Apply a morphological dilation to the raster using a disk footprint.

    Morphological dilation sets the value of a pixel to the maximum over all pixel
    values within a local neighborhood centered about it.

    Dilation enlarges bright regions and shrinks dark regions. This is useful e.g.
    to find a region nearby a steep edge after applying a Sobel filter.

    The radius parameter controls the dilation extent. This is rounded up to the
    nearest integer multiple of the cell size, since dilation is a discrete
    operation. To reduce inaccuracies due to this rounding, consider resampling
    the raster to a smaller cell size before applying dilation.

    NaN values in the original raster are preserved in their original locations.
    They are temporarily filled during dilation to avoid spreading into valid data,
    then restored after the operation completes. The output raster maintains the
    same shape as the input.

    Args:
        radius: Radius of the disk footprint used in dilation, in units of
                geographic coordinate distance (e.g. meters).

    Returns:
        New raster with dilated values. NaN locations are preserved from the
        original raster.
    """
    from skimage import morphology

    if not self.has_square_cells:
        msg = "Dilate currently only supports rasters with square cells."
        raise NotImplementedError(msg)

    cell_size = self.square_cell_size

    # Round up to nearest cell
    cell_radius = int(np.ceil(radius / cell_size))

    # Calculate actual radius based on rounded cell count
    radius_m = cell_radius * cell_size

    # Store original NaN mask and shape
    original_nan_mask = np.isnan(self.arr)
    original_shape = self.arr.shape

    # Handle all-NaN case
    if np.all(original_nan_mask):
        return self.model_copy()

    # Pad the raster with non-consequential values to avoid edge effects
    fill_val = self.min() - 1.0
    new_raster = self.model_copy()
    new_raster = new_raster.pad(width=radius_m, value=fill_val)

    # Replace NaNs with fill_val to avoid issues during dilation
    new_raster.arr[np.isnan(new_raster.arr)] = fill_val
    new_raster.arr = morphology.dilation(
        new_raster.arr, morphology.disk(cell_radius)
    )

    # Crop back to original size using array slicing
    new_raster.arr = new_raster.arr[
        cell_radius : cell_radius + original_shape[0],
        cell_radius : cell_radius + original_shape[1],
    ]
    # Restore original metadata
    new_raster = self.__class__(arr=new_raster.arr, meta=self.meta)

    # Restore original NaN values
    new_raster.arr[original_nan_mask] = np.nan

    return new_raster

example() classmethod

Create an example Raster.

Source code in src/rastr/raster.py

@classmethod
def example(cls) -> Self:
    """Create an example Raster."""
    # Peaks dataset style example
    n = 256
    x = np.linspace(-3, 3, n)
    y = np.linspace(-3, 3, n)
    x, y = np.meshgrid(x, y)
    z = np.exp(-(x**2) - y**2) * np.sin(3 * np.sqrt(x**2 + y**2))
    arr = z.astype(np.float32)

    raster_meta = RasterMeta.example()
    return cls(arr=arr, raster_meta=raster_meta)

exp()

Compute the exponential of the raster.

Returns a new raster with the exponential (e^x) of each cell. The original raster is not modified.

Returns:

Type	Description
Self	A new Raster instance with the exponential values.

Source code in src/rastr/raster.py

def exp(self) -> Self:
    """Compute the exponential of the raster.

    Returns a new raster with the exponential (e^x) of each cell. The original
    raster is not modified.

    Returns:
        A new Raster instance with the exponential values.
    """
    cls = self.__class__
    return cls(arr=np.exp(self.arr), raster_meta=self.raster_meta)

explore(*, m=None, opacity=1.0, colormap='viridis', cbar_label=None, vmin=None, vmax=None)

Display the raster on a folium map.

Source code in src/rastr/raster.py

def explore(  # noqa: PLR0913 c.f. geopandas.explore which also has many input args
    self,
    *,
    m: Map | None = None,
    opacity: float = 1.0,
    colormap: str
    | Callable[[float], tuple[float, float, float, float]] = "viridis",
    cbar_label: str | None = None,
    vmin: float | None = None,
    vmax: float | None = None,
) -> Map:
    """Display the raster on a folium map."""
    if not FOLIUM_INSTALLED:
        msg = "The 'folium' package is required for 'explore()'."
        raise ImportError(msg)

    import folium.raster_layers

    if m is None:
        m = folium.Map()

    if vmin is not None and vmax is not None and vmax <= vmin:
        msg = "'vmin' must be less than 'vmax'."
        raise ValueError(msg)

    cmap_func = _get_colormap_function(colormap)

    # Transform bounds to WGS84 using pyproj directly
    wgs84_crs = CRS.from_epsg(4326)
    transformer = Transformer.from_crs(
        self.raster_meta.crs, wgs84_crs, always_xy=True
    )

    # Get the corner points of the bounding box
    raster_xmin, raster_ymin, raster_xmax, raster_ymax = self.bounds
    corner_points = [
        (raster_xmin, raster_ymin),
        (raster_xmin, raster_ymax),
        (raster_xmax, raster_ymax),
        (raster_xmax, raster_ymin),
    ]

    # Transform all corner points to WGS84
    transformed_points = [transformer.transform(x, y) for x, y in corner_points]

    # Find the bounding box of the transformed points
    transformed_xs, transformed_ys = zip(*transformed_points, strict=True)
    xmin, xmax = min(transformed_xs), max(transformed_xs)
    ymin, ymax = min(transformed_ys), max(transformed_ys)

    # Normalize the array to [0, 1] for colormap mapping
    _vmin, _vmax = _get_vmin_vmax(self, vmin=vmin, vmax=vmax)
    arr = self.normalize(vmin=_vmin, vmax=_vmax).arr

    # Finally, need to determine whether to flip the image based on negative Affine
    # coefficients
    flip_x = self.raster_meta.transform.a < 0
    flip_y = self.raster_meta.transform.e > 0
    if flip_x:
        arr = np.flip(arr, axis=1)
    if flip_y:
        arr = np.flip(arr, axis=0)

    bounds = [[ymin, xmin], [ymax, xmax]]
    img = folium.raster_layers.ImageOverlay(
        image=arr,
        bounds=bounds,
        opacity=opacity,
        colormap=cmap_func,
        mercator_project=True,
    )

    img.add_to(m)

    # Add a colorbar legend
    if BRANCA_INSTALLED:
        cbar = _map_colorbar(colormap=cmap_func, vmin=_vmin, vmax=_vmax)
        if cbar_label:
            cbar.caption = cbar_label
        cbar.add_to(m)

    m.fit_bounds(bounds)

    return m

extrapolate(method='nearest')

Extrapolate the raster data to fill NaN values.

See also fillna() for filling NaN values with a specific value.

If the raster is all-NaN, this method will return a copy of the raster without changing the NaN values.

Parameters:

Name	Type	Description	Default
method	Literal['nearest']	The method to use for extrapolation. Currently only 'nearest' is supported, which fills NaN values with the nearest non-NaN value.	'nearest'

Source code in src/rastr/raster.py

def extrapolate(self, method: Literal["nearest"] = "nearest") -> Self:
    """Extrapolate the raster data to fill NaN values.

    See also `fillna()` for filling NaN values with a specific value.

    If the raster is all-NaN, this method will return a copy of the raster without
    changing the NaN values.

    Args:
        method: The method to use for extrapolation. Currently only 'nearest' is
                supported, which fills NaN values with the nearest non-NaN value.
    """
    if method not in ("nearest",):
        msg = f"Unsupported extrapolation method: {method}"
        raise NotImplementedError(msg)

    raster = self.model_copy()
    raster.arr = fillna_nearest_neighbours(arr=self.arr)

    return raster

fillna(value)

Fill NaN values in the raster with a specified value.

See also extrapolate() for filling NaN values using extrapolation from data.

Source code in src/rastr/raster.py

def fillna(self, value: float) -> Self:
    """Fill NaN values in the raster with a specified value.

    See also `extrapolate()` for filling NaN values using extrapolation from data.
    """
    filled_arr = np.nan_to_num(self.arr, nan=value)
    new_raster = self.model_copy()
    new_raster.arr = filled_arr
    return new_raster

full_like(other, *, fill_value) classmethod

Create a raster with the same metadata as another but filled with a constant.

Parameters:

Name	Type	Description	Default
other	Raster	The raster to copy metadata from.	required
fill_value	float	The constant value to fill all cells with.	required

Returns:

Type	Description
Self	A new raster with the same shape and metadata as other, but with all cells
Self	set to fill_value.

Source code in src/rastr/raster.py

@classmethod
def full_like(cls, other: Raster, *, fill_value: float) -> Self:
    """Create a raster with the same metadata as another but filled with a constant.

    Args:
        other: The raster to copy metadata from.
        fill_value: The constant value to fill all cells with.

    Returns:
        A new raster with the same shape and metadata as `other`, but with all cells
        set to `fill_value`.
    """
    arr = np.full(other.shape, fill_value, dtype=np.float32)
    return cls(arr=arr, raster_meta=other.raster_meta)

gdf(name='value')

Create a GeoDataFrame representation of the raster.

Alias for as_geodataframe().

Source code in src/rastr/raster.py

def gdf(self, name: str = "value") -> gpd.GeoDataFrame:
    """Create a GeoDataFrame representation of the raster.

    Alias for `as_geodataframe()`.
    """
    return self.as_geodataframe(name=name)

get_xy()

Get the x and y coordinates of the raster cell centres in meshgrid format.

Returns the coordinates of the cell centres as two separate 2D arrays in meshgrid format, where each array has the same shape as the raster data array.

Returns:

Type	Description
NDArray[float64]	A tuple of (x, y) coordinate arrays where:
NDArray[float64]	x: 2D array of x-coordinates of cell centres
tuple[NDArray[float64], NDArray[float64]]	y: 2D array of y-coordinates of cell centres
tuple[NDArray[float64], NDArray[float64]]	Both arrays have the same shape as the raster data array.

Source code in src/rastr/raster.py

def get_xy(self) -> tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Get the x and y coordinates of the raster cell centres in meshgrid format.

    Returns the coordinates of the cell centres as two separate 2D arrays in
    meshgrid format, where each array has the same shape as the raster data array.

    Returns:
        A tuple of (x, y) coordinate arrays where:
        - x: 2D array of x-coordinates of cell centres
        - y: 2D array of y-coordinates of cell centres
        Both arrays have the same shape as the raster data array.
    """
    coords = self.raster_meta.get_cell_centre_coords(self.arr.shape)
    return coords[:, :, 0], coords[:, :, 1]

is_like(other)

Check if two Raster objects have the same metadata and shape.

Parameters:

Name	Type	Description	Default
other	Raster	Another Raster to compare with.	required

Returns:

Type	Description
bool	True if both rasters have the same meta and shape attributes.

Source code in src/rastr/raster.py

def is_like(self, other: Raster) -> bool:
    """Check if two Raster objects have the same metadata and shape.

    Args:
        other: Another Raster to compare with.

    Returns:
        True if both rasters have the same meta and shape attributes.
    """
    return self.meta == other.meta and self.shape == other.shape

log()

Compute the natural logarithm of the raster.

Returns a new raster with the natural logarithm of each cell. The original raster is not modified.

Returns:

Type	Description
Self	A new Raster instance with the natural logarithm values.

Source code in src/rastr/raster.py

def log(self) -> Self:
    """Compute the natural logarithm of the raster.

    Returns a new raster with the natural logarithm of each cell. The original
    raster is not modified.

    Returns:
        A new Raster instance with the natural logarithm values.
    """
    cls = self.__class__
    return cls(arr=np.log(self.arr), raster_meta=self.raster_meta)

max()

Get the maximum value in the raster, ignoring NaN values.

Returns:

Type	Description
float	The maximum value in the raster. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def max(self) -> float:
    """Get the maximum value in the raster, ignoring NaN values.

    Returns:
        The maximum value in the raster. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanmax(self.arr))

mean()

Get the mean value in the raster, ignoring NaN values.

Returns:

Type	Description
float	The mean value in the raster. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def mean(self) -> float:
    """Get the mean value in the raster, ignoring NaN values.

    Returns:
        The mean value in the raster. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanmean(self.arr))

median()

Get the median value in the raster, ignoring NaN values.

This is equivalent to quantile(0.5).

Returns:

Type	Description
float	The median value in the raster. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def median(self) -> float:
    """Get the median value in the raster, ignoring NaN values.

    This is equivalent to quantile(0.5).

    Returns:
        The median value in the raster. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanmedian(self.arr))

min()

Get the minimum value in the raster, ignoring NaN values.

Returns:

Type	Description
float	The minimum value in the raster. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def min(self) -> float:
    """Get the minimum value in the raster, ignoring NaN values.

    Returns:
        The minimum value in the raster. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanmin(self.arr))

normalize(*, vmin=None, vmax=None)

Normalize the raster values to the range [0, 1].

If custom vmin and vmax values are provided, values below vmin will be set to 0, and values above vmax will be set to 1.

Parameters:

Name	Type	Description	Default
vmin	float | None	Minimum value for normalization. Values below this will be set to 0. If None, the minimum value in the array is used.	None
vmax	float | None	Maximum value for normalization. Values above this will be set to 1. If None, the maximum value in the array is used.	None

Source code in src/rastr/raster.py

def normalize(
    self, *, vmin: float | None = None, vmax: float | None = None
) -> Self:
    """Normalize the raster values to the range [0, 1].

    If custom vmin and vmax values are provided, values below vmin will be set to 0,
    and values above vmax will be set to 1.

    Args:
        vmin: Minimum value for normalization. Values below this will be set to 0.
              If None, the minimum value in the array is used.
        vmax: Maximum value for normalization. Values above this will be set to 1.
              If None, the maximum value in the array is used.
    """
    _vmin, _vmax = _get_vmin_vmax(self, vmin=vmin, vmax=vmax)

    arr = self.arr.copy()
    if _vmax > _vmin:
        arr = (arr - _vmin) / (_vmax - _vmin)
        arr = np.clip(arr, 0, 1)
    else:
        arr = np.zeros_like(arr)
    return self.__class__(arr=arr, raster_meta=self.raster_meta)

pad(width, *, value=np.nan)

Extend the raster by adding a constant fill value around the edges.

By default, the padding value is NaN, but this can be changed via the value parameter.

This grows the raster by adding padding around all edges. New cells are filled with the constant value.

If the width is not an exact multiple of the cell size, the padding may be slightly larger than the specified width, i.e. the value is rounded up to the nearest whole number of cells. For non-square cells, padding is computed independently in x and y directions.

Parameters:

Name	Type	Description	Default
width	float	The width of the padding, in the same units as the raster CRS (e.g. meters). This defines how far from the edge the padding extends.	required
value	float	The constant value to use for padding. Default is NaN.	nan

Source code in src/rastr/raster.py

def pad(self, width: float, *, value: float = np.nan) -> Self:
    """Extend the raster by adding a constant fill value around the edges.

    By default, the padding value is NaN, but this can be changed via the
    `value` parameter.

    This grows the raster by adding padding around all edges. New cells are
    filled with the constant `value`.

    If the width is not an exact multiple of the cell size, the padding may be
    slightly larger than the specified width, i.e. the value is rounded up to
    the nearest whole number of cells. For non-square cells, padding is computed
    independently in x and y directions.

    Args:
        width: The width of the padding, in the same units as the raster CRS
               (e.g. meters). This defines how far from the edge the padding
               extends.
        value: The constant value to use for padding. Default is NaN.
    """
    cell_width, cell_height = self.raster_meta.cell_size

    # Calculate number of cells to pad in each direction
    pad_cols = int(np.ceil(width / cell_width))
    pad_rows = int(np.ceil(width / cell_height))

    # Get current bounds
    xmin, ymin, xmax, ymax = self.bounds

    # Calculate new bounds with padding
    new_xmin = xmin - (pad_cols * cell_width)
    new_ymin = ymin - (pad_rows * cell_height)
    new_xmax = xmax + (pad_cols * cell_width)
    new_ymax = ymax + (pad_rows * cell_height)

    # Create padded array
    new_height = self.arr.shape[0] + 2 * pad_rows
    new_width = self.arr.shape[1] + 2 * pad_cols

    # Create new array filled with the padding value
    padded_arr = np.full((new_height, new_width), value, dtype=self.arr.dtype)

    # Copy original array into the center of the padded array
    padded_arr[
        pad_rows : pad_rows + self.arr.shape[0],
        pad_cols : pad_cols + self.arr.shape[1],
    ] = self.arr

    # Create new transform for the padded raster
    new_transform = rasterio.transform.from_bounds(
        west=new_xmin,
        south=new_ymin,
        east=new_xmax,
        north=new_ymax,
        width=new_width,
        height=new_height,
    )

    # Create new raster metadata
    new_meta = RasterMeta(
        crs=self.raster_meta.crs,
        transform=new_transform,
    )

    return self.__class__(arr=padded_arr, raster_meta=new_meta)

plot(*, ax=None, cbar_label=None, basemap=False, cmap='viridis', suppressed=tuple(), **kwargs)

Plot the raster on a matplotlib axis.

Parameters:

Name	Type	Description	Default
ax	Axes | None	A matplotlib axes object to plot on. If None, a new figure will be created.	None
cbar_label	str | None	Label for the colorbar. If None, no label is added.	None
basemap	bool	Whether to add a basemap. Currently not implemented.	False
cmap	str	Colormap to use for the plot.	'viridis'
suppressed	Collection[float] | float	Values to suppress from the plot (i.e. not display). This can be useful for zeroes especially.	tuple()
**kwargs	Any	Additional keyword arguments to pass to rasterio.plot.show(). This includes parameters like alpha for transparency, and vmin and vmax for controlling the color scale limits.	{}

Source code in src/rastr/raster.py

def plot(
    self,
    *,
    ax: Axes | None = None,
    cbar_label: str | None = None,
    basemap: bool = False,
    cmap: str = "viridis",
    suppressed: Collection[float] | float = tuple(),
    **kwargs: Any,
) -> Axes:
    """Plot the raster on a matplotlib axis.

    Args:
        ax: A matplotlib axes object to plot on. If None, a new figure will be
            created.
        cbar_label: Label for the colorbar. If None, no label is added.
        basemap: Whether to add a basemap. Currently not implemented.
        cmap: Colormap to use for the plot.
        suppressed: Values to suppress from the plot (i.e. not display). This can be
                    useful for zeroes especially.
        **kwargs: Additional keyword arguments to pass to `rasterio.plot.show()`.
                  This includes parameters like `alpha` for transparency, and `vmin`
                  and `vmax` for controlling the color scale limits.
    """
    if not MATPLOTLIB_INSTALLED:
        msg = "The 'matplotlib' package is required for 'plot()'."
        raise ImportError(msg)

    from matplotlib import pyplot as plt
    from mpl_toolkits.axes_grid1 import make_axes_locatable

    suppressed = np.array(suppressed)

    if ax is None:
        _, _ax = plt.subplots()
        _ax: Axes
        ax = _ax

    if basemap:
        msg = "Basemap plotting is not yet implemented."
        raise NotImplementedError(msg)

    model = self.model_copy()
    model.arr = model.arr.copy()

    # Get extent of the unsuppressed values in array index coordinates
    suppressed_mask = np.isin(model.arr, suppressed)
    (x_unsuppressed,) = np.nonzero((~suppressed_mask).any(axis=0))
    (y_unsuppressed,) = np.nonzero((~suppressed_mask).any(axis=1))

    if len(x_unsuppressed) == 0 or len(y_unsuppressed) == 0:
        msg = "Raster contains no unsuppressed values; cannot plot."
        raise ValueError(msg)

    # N.B. these are array index coordinates, so np.min and np.max are safe since
    # they cannot encounter NaN values.
    min_x_unsuppressed = np.min(x_unsuppressed)
    max_x_unsuppressed = np.max(x_unsuppressed)
    min_y_unsuppressed = np.min(y_unsuppressed)
    max_y_unsuppressed = np.max(y_unsuppressed)

    # Transform to raster CRS
    x1, y1 = self.raster_meta.transform * (  # type: ignore[reportAssignmentType] overloaded tuple size in affine
        min_x_unsuppressed,
        min_y_unsuppressed,
    )
    x2, y2 = self.raster_meta.transform * (  # type: ignore[reportAssignmentType]
        max_x_unsuppressed + 1,
        max_y_unsuppressed + 1,
    )
    xmin, xmax = sorted([x1, x2])
    ymin, ymax = sorted([y1, y2])

    model.arr[suppressed_mask] = np.nan

    img, *_ = model.rio_show(ax=ax, cmap=cmap, with_bounds=True, **kwargs)

    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)

    ax.set_aspect("equal", "box")
    ax.set_yticklabels([])
    ax.set_xticklabels([])
    ax.tick_params(left=False, bottom=False, top=False, right=False)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    fig = ax.get_figure()
    if fig is not None:
        fig.colorbar(img, label=cbar_label, cax=cax)
    return ax

quantile(q)

Get the specified quantile value in the raster, ignoring NaN values.

Parameters:

Name	Type	Description	Default
q	float	Quantile to compute, must be between 0 and 1 inclusive.	required

Returns:

Type	Description
float	The quantile value. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def quantile(self, q: float) -> float:
    """Get the specified quantile value in the raster, ignoring NaN values.

    Args:
        q: Quantile to compute, must be between 0 and 1 inclusive.

    Returns:
        The quantile value. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanquantile(self.arr, q))

read_file(filename, crs=None) classmethod

Read raster data from a file and return an in-memory Raster object.

Parameters:

Name	Type	Description	Default
filename	Path | str	Path to the raster file.	required
crs	CRS | str | None	Optional coordinate reference system to override the file's CRS.	None

Source code in src/rastr/raster.py

@classmethod
def read_file(cls, filename: Path | str, crs: CRS | str | None = None) -> Self:
    """Read raster data from a file and return an in-memory Raster object.

    Args:
        filename: Path to the raster file.
        crs: Optional coordinate reference system to override the file's CRS.
    """
    # Import here to avoid circular import (rastr.io imports Raster)
    from rastr.io_ import read_raster_inmem  # noqa: PLC0415

    return read_raster_inmem(filename, crs=crs, cls=cls)

read_mosaic_dir(mosaic_dir, *, glob='*.tif', crs=None) classmethod

Read a raster mosaic from a directory and return an in-memory Raster object.

Parameters:

Name	Type	Description	Default
mosaic_dir	Path | str	Path to the directory containing raster files.	required
glob	str	Glob pattern to match raster files. Default is "*.tif".	'*.tif'
crs	CRS | None	Optional coordinate reference system to override the files' CRS.	None

Returns:

Type	Description
Self	A Raster object containing the mosaicked data from all matching files.

Source code in src/rastr/raster.py

@classmethod
def read_mosaic_dir(
    cls, mosaic_dir: Path | str, *, glob: str = "*.tif", crs: CRS | None = None
) -> Self:
    """Read a raster mosaic from a directory and return an in-memory Raster object.

    Args:
        mosaic_dir: Path to the directory containing raster files.
        glob: Glob pattern to match raster files. Default is "*.tif".
        crs: Optional coordinate reference system to override the files' CRS.

    Returns:
        A Raster object containing the mosaicked data from all matching files.
    """
    # Import here to avoid circular import (rastr.io imports Raster)
    from rastr.io_ import read_raster_mosaic_inmem  # noqa: PLC0415

    raster_obj = read_raster_mosaic_inmem(mosaic_dir, glob=glob, crs=crs)
    return cls(arr=raster_obj.arr, raster_meta=raster_obj.raster_meta)

replace(to_replace, value=None)

Replace values in the raster with other values.

Creates a new raster with the specified values replaced. This is useful for operations like replacing zeros with NaNs, or vice versa.

The method supports two interfaces: 1. Single replacement: raster.replace(to_replace=0, value=np.nan) 2. Multiple replacements using a dictionary: raster.replace({0: np.nan, -999: np.nan})

Parameters:

Name	Type	Description	Default
to_replace	float | dict[float, float]	Value to be replaced, or a dictionary mapping values to their replacements.	required
value	float | None	Replacement value. Required when to_replace is a float, must be None when to_replace is a dict.	None

Examples:

>>> # Replace a single value
>>> raster.replace(to_replace=0, value=np.nan)
>>> # Replace multiple values
>>> raster.replace({0: np.nan, -999: np.nan})

Source code in src/rastr/raster.py

def replace(
    self, to_replace: float | dict[float, float], value: float | None = None
) -> Self:
    """Replace values in the raster with other values.

    Creates a new raster with the specified values replaced. This is useful for
    operations like replacing zeros with NaNs, or vice versa.

    The method supports two interfaces:
    1. Single replacement: `raster.replace(to_replace=0, value=np.nan)`
    2. Multiple replacements using a dictionary:
       `raster.replace({0: np.nan, -999: np.nan})`

    Args:
        to_replace: Value to be replaced, or a dictionary mapping values to
                    their replacements.
        value: Replacement value. Required when to_replace is a float, must be
               None when to_replace is a dict.

    Examples:
        >>> # Replace a single value
        >>> raster.replace(to_replace=0, value=np.nan)
        >>> # Replace multiple values
        >>> raster.replace({0: np.nan, -999: np.nan})
    """
    # Determine the replacement map
    if isinstance(to_replace, dict):
        if value is not None:
            msg = "value must be None when to_replace is a dict"
            raise ValueError(msg)
        map_ = to_replace
    else:
        if value is None:
            msg = "value must be specified when to_replace is a float"
            raise ValueError(msg)
        map_ = {to_replace: value}

    # Start with a copy of the array
    replaced_arr = self.arr.copy()

    # Check if we need to convert to float (if assigning NaN to non-float array)
    needs_float = any(
        np.isnan(new_val) for new_val in map_.values()
    ) and not np.issubdtype(replaced_arr.dtype, np.floating)
    if needs_float:
        replaced_arr = replaced_arr.astype(float)

    # Apply each replacement based on the original array values
    # to prevent chained replacements
    for old_val, new_val in map_.items():
        # Handle NaN specially since NaN != NaN
        if np.isnan(old_val):
            mask = np.isnan(self.arr)
        else:
            mask = self.arr == old_val

        replaced_arr[mask] = new_val

    new_raster = self.model_copy()
    new_raster.arr = replaced_arr
    return new_raster

replace_polygon(polygon, value=None)

Replace values within the specified polygon(s) with other values.

Creates a new raster with the specified values replaced. This is useful for operations like masking or setting regions to NaN.

The method supports two interfaces: 1. Single replacement: raster.replace_polygon(polygon1, value=np.nan) 2. Multiple replacements using a dictionary: raster.replace_polygon({polygon1: 0, polygon2: 1})

Parameters:

Name	Type	Description	Default
polygon	BaseGeometry | dict[BaseGeometry, float]	Geometry to replace, or dict mapping geometries to values.	required
value	float | None	Replacement value. Required if polygon is a geometry, None if polygon is a dict.	None

Examples:

>>> # Replace a single polygon
>>> raster.replace_polygon(polygon1, value=np.nan)
>>> # Replace multiple polygons
>>> raster.replace_polygon({polygon1: 0, polygon2: 1})

Source code in src/rastr/raster.py

def replace_polygon(
    self,
    polygon: BaseGeometry | dict[BaseGeometry, float],
    value: float | None = None,
) -> Self:
    """Replace values within the specified polygon(s) with other values.

    Creates a new raster with the specified values replaced. This is useful for
    operations like masking or setting regions to NaN.

    The method supports two interfaces:
    1. Single replacement: `raster.replace_polygon(polygon1, value=np.nan)`
    2. Multiple replacements using a dictionary:
       `raster.replace_polygon({polygon1: 0, polygon2: 1})`

    Args:
        polygon: Geometry to replace, or dict mapping geometries to values.
        value: Replacement value. Required if polygon is a geometry, None if polygon
               is a dict.

    Examples:
        >>> # Replace a single polygon
        >>> raster.replace_polygon(polygon1, value=np.nan)
        >>> # Replace multiple polygons
        >>> raster.replace_polygon({polygon1: 0, polygon2: 1})
    """
    # Normalize input to dict format
    if isinstance(polygon, dict):
        if value is not None:
            msg = "value must be None when polygon is a dict"
            raise ValueError(msg)
        replacements = polygon
    else:
        if value is None:
            msg = "value must be specified when polygon is a geometry"
            raise ValueError(msg)
        replacements = {polygon: value}

    # Validate all geometries upfront
    for geom in replacements.keys():
        if not isinstance(geom, Polygon | MultiPolygon):
            msg = (
                f"Only Polygon and MultiPolygon geometries are supported, "
                f"got {type(geom).__name__}"
            )
            raise TypeError(msg)

    # Create copy and convert to float if needed
    raster = self.model_copy()
    needs_float = any(
        isinstance(v, float) or (v is not None and np.isnan(v))
        for v in replacements.values()
    )
    if needs_float and not np.issubdtype(raster.arr.dtype, np.floating):
        raster.arr = raster.arr.astype(float)

    # Apply all replacements
    for geom, val in replacements.items():
        mask_raster = self._polygon_indicator(geom)
        raster.arr = np.where(mask_raster.arr, val, raster.arr)

    return raster

resample(cell_size, *, method='bilinear')

Resample the raster data to a new resolution.

If the new cell size is not an exact multiple of the current cell size, the overall raster bounds may increase slightly. The affine transform will keep the same shift, i.e. the top-left corner of the raster will remain in the same' coordinate location. A corollary is that the overall centre of the raster bounds will not necessary be the same as the original raster.

Parameters:

Name	Type	Description	Default
cell_size	tuple[float, float] | float	The desired cell size for the resampled raster. This can be a	required
method	Literal['bilinear']	The resampling method to use. Only 'bilinear' is supported.	'bilinear'

Source code in src/rastr/raster.py

def resample(
    self,
    cell_size: tuple[float, float] | float,
    *,
    method: Literal["bilinear"] = "bilinear",
) -> Self:
    """Resample the raster data to a new resolution.

    If the new cell size is not an exact multiple of the current cell size, the
    overall raster bounds may increase slightly. The affine transform will keep
    the same shift, i.e. the top-left corner of the raster will remain in the same'
    coordinate location. A corollary is that the overall centre of the raster bounds
    will not necessary be the same as the original raster.

    Args:
        cell_size: The desired cell size for the resampled raster. This can be a
        single float for square cells, or a tuple of (cell_width, cell_height) for
        rectangular cells.
        method: The resampling method to use. Only 'bilinear' is supported.
    """
    if method not in ("bilinear",):
        msg = f"Unsupported resampling method: {method}"
        raise NotImplementedError(msg)

    target_cell_width, target_cell_height = _ensure_pair(cell_size)
    source_cell_width, source_cell_height = self.raster_meta.cell_size

    x_factor = source_cell_width / target_cell_width
    y_factor = source_cell_height / target_cell_height

    cls = self.__class__
    # Use the rasterio dataset with proper context management
    with self.to_rasterio_dataset() as dataset:
        # N.B. the new height and width may increase slightly.
        new_height = int(np.ceil(dataset.height * y_factor))
        new_width = int(np.ceil(dataset.width * x_factor))

        # Resample via rasterio
        (new_arr,) = dataset.read(  # Assume exactly one band
            out_shape=(dataset.count, new_height, new_width),
            resampling=Resampling.bilinear,
        )

        # Create new RasterMeta with the exact requested cell size
        new_raster_meta = RasterMeta(
            transform=Affine(
                target_cell_width,
                dataset.transform.b,
                dataset.transform.c,
                dataset.transform.d,
                -target_cell_height,
                dataset.transform.f,
            ),
            crs=self.raster_meta.crs,
        )

        return cls(arr=new_arr, raster_meta=new_raster_meta)

rio_show(**kwargs)

Plot the raster using rasterio's built-in plotting function.

This is useful for lower-level access to rasterio's plotting capabilities. Generally, the plot() method is preferred for most use cases.

Parameters:

Name	Type	Description	Default
**kwargs	Any	Keyword arguments to pass to rasterio.plot.show(). This includes parameters like alpha for transparency, and with_bounds to control whether to plot in spatial coordinates or array index coordinates.	{}

Source code in src/rastr/raster.py

def rio_show(self, **kwargs: Any) -> list[AxesImage]:
    """Plot the raster using rasterio's built-in plotting function.

    This is useful for lower-level access to rasterio's plotting capabilities.
    Generally, the `plot()` method is preferred for most use cases.

    Args:
        **kwargs: Keyword arguments to pass to `rasterio.plot.show()`. This includes
                  parameters like `alpha` for transparency, and `with_bounds` to
                  control whether to plot in spatial coordinates or array index
                  coordinates.
    """
    with self.to_rasterio_dataset() as dataset:
        return rasterio.plot.show(dataset, **kwargs).get_images()

sample(xy, *, na_action='raise')

sample(xy: Collection[tuple[float, float]] | Collection[Point] | ArrayLike, *, na_action: Literal['raise', 'ignore'] = 'raise') -> NDArray

sample(xy: tuple[float, float] | Point, *, na_action: Literal['raise', 'ignore'] = 'raise') -> float

Sample raster values at GeoSeries locations and return sampled values.

Parameters:

Name	Type	Description	Default
xy	Collection[tuple[float, float]] | Collection[Point] | ArrayLike | tuple[float, float] | Point | GeoDataFrame	A list of (x, y) coordinates, shapely Point objects, or a GeoDataFrame containing Point geometries to sample the raster at.	required
na_action	Literal['raise', 'ignore']	Action to take when a NaN value is encountered in the input xy. Options are "raise" (raise an error) or "ignore" (replace with NaN).	'raise'

Returns:

Type	Description
NDArray | float	A list of sampled raster values for each geometry in the GeoSeries.

Source code in src/rastr/raster.py

def sample(
    self,
    xy: Collection[tuple[float, float]]
    | Collection[Point]
    | ArrayLike
    | tuple[float, float]
    | Point
    | gpd.GeoDataFrame,
    *,
    na_action: Literal["raise", "ignore"] = "raise",
) -> NDArray | float:
    """Sample raster values at GeoSeries locations and return sampled values.

    Args:
        xy: A list of (x, y) coordinates, shapely Point objects, or a
            GeoDataFrame containing Point geometries to sample the raster at.
        na_action: Action to take when a NaN value is encountered in the input xy.
                   Options are "raise" (raise an error) or "ignore" (replace with
                   NaN).

    Returns:
        A list of sampled raster values for each geometry in the GeoSeries.
    """
    # If this function is too slow, consider the optimizations detailed here:
    # https://rdrn.me/optimising-sampling/

    import geopandas as gpd

    if isinstance(xy, gpd.GeoDataFrame):
        xy = _gdf_to_xy(xy, self.crs)
        singleton = False

    # Convert shapely Points to coordinate tuples if needed
    elif isinstance(xy, Point):
        xy = [(xy.x, xy.y)]
        singleton = True
    elif (
        isinstance(xy, Collection)
        and len(xy) > 0
        and isinstance(next(iter(xy)), Point)
    ):
        xy = [(point.x, point.y) for point in xy]  # pyright: ignore[reportAttributeAccessIssue]
        singleton = False
    elif (
        isinstance(xy, tuple)
        and len(xy) == 2
        and isinstance(next(iter(xy)), (float, int))
    ):
        xy = [xy]  # pyright: ignore[reportAssignmentType]
        singleton = True
    else:
        singleton = False

    xy = np.asarray(xy, dtype=float)

    if len(xy) == 0:
        # Short-circuit
        return np.array([], dtype=float)

    # Create in-memory rasterio dataset from the incumbent Raster object
    with self.to_rasterio_dataset() as dataset:
        if dataset.count != 1:
            msg = "Only single band rasters are supported."
            raise NotImplementedError(msg)

        xy_arr = np.array(xy)

        # Determine the indexes of any x,y coordinates where either is NaN.
        # We will drop these indexes for the purposes of calling .sample, but
        # then we will add NaN values back in at the end, inserting NaN into the
        # results array.
        xy_is_nan = np.isnan(xy_arr).any(axis=1)
        xy_nan_idxs = list(np.atleast_1d(np.squeeze(np.nonzero(xy_is_nan))))
        xy_arr = xy_arr[~xy_is_nan]

        if na_action == "raise" and len(xy_nan_idxs) > 0:
            nan_error_msg = "NaN value found in input coordinates"
            raise ValueError(nan_error_msg)

        # Sample the raster in-memory dataset (e.g. PGA values) at the coordinates
        samples = list(
            rasterio.sample.sample_gen(
                dataset,
                xy_arr,
                indexes=1,  # Single band raster, N.B. rasterio is 1-indexed
                masked=True,
            )
        )

        # Convert the sampled values to a NumPy array and set masked values to NaN
        raster_values = np.array(
            [s.data[0] if not numpy.ma.getmask(s) else np.nan for s in samples]
        ).astype(float)

        if len(xy_nan_idxs) > 0:
            # Insert NaN values back into the results array
            # This is tricky because all the indexes get offset once we remove
            # elements.
            offset_xy_nan_idxs = xy_nan_idxs - np.arange(len(xy_nan_idxs))
            raster_values = np.insert(
                raster_values,
                offset_xy_nan_idxs,
                np.nan,
                axis=0,
            )

    if singleton:
        (raster_value,) = raster_values
        return raster_value

    return raster_values

set_crs(crs, *, allow_override=False)

Set the CRS of the raster without reprojecting.

To reproject the raster data to a different CRS, use to_crs() instead.

Parameters:

Name	Type	Description	Default
crs	CRS | str	The coordinate reference system to assign. Can be a CRS object or a string that can be parsed by pyproj (e.g., 'EPSG:4326').	required
allow_override	bool	If False (default), raises an error if the raster already has a CRS that differs from the one being set. If True, allows overriding any existing CRS.	False

Returns:

Type	Description
Self	A new Raster instance with the updated CRS and unchanged array data.

Raises:

Type	Description
ValueError	If the raster already has a CRS that differs from the one being set and allow_override is False.

Source code in src/rastr/raster.py

def set_crs(self, crs: CRS | str, *, allow_override: bool = False) -> Self:
    """Set the CRS of the raster without reprojecting.

    To reproject the raster data to a different CRS, use `to_crs()` instead.

    Args:
        crs: The coordinate reference system to assign. Can be a CRS object or a
             string that can be parsed by pyproj (e.g., 'EPSG:4326').
        allow_override: If False (default), raises an error if the raster already
                       has a CRS that differs from the one being set. If True,
                       allows overriding any existing CRS.

    Returns:
        A new Raster instance with the updated CRS and unchanged array data.

    Raises:
        ValueError: If the raster already has a CRS that differs from the one
                   being set and allow_override is False.

    """
    crs_obj = CRS.from_user_input(crs)

    # Check if raster already has a different CRS
    if (
        self.raster_meta.crs is not None
        and self.raster_meta.crs != crs_obj
        and not allow_override
    ):
        msg = (
            f"The raster already has a CRS ({self.raster_meta.crs}) which is "
            f"different from the one provided ({crs_obj}). Use "
            f"allow_override=True to override."
        )
        raise ValueError(msg)

    new_meta = RasterMeta(
        crs=crs_obj,
        transform=self.raster_meta.transform,
    )
    return self.__class__(arr=self.arr, raster_meta=new_meta)

set_origin(*, x=None, y=None)

Set the transform origin without modifying pixel data.

The origin is the corner of the first pixel in the raster grid — the translation component (c, f) of the affine transform. For north-up rasters this corresponds to (xmin, ymax); for south-up rasters it corresponds to (xmin, ymin).

Unspecified axes retain their current value. If neither axis is provided, returns an unchanged copy.

Parameters:

Name	Type	Description	Default
x	float | None	New x-coordinate of the grid origin. If None, keeps current.	None
y	float | None	New y-coordinate of the grid origin. If None, keeps current.	None

Returns:

Type	Description
Self	A new Raster with the updated transform origin and unchanged array data.

Example

Normalize a raster with non-standard longitude (e.g., -200° to -170°) to standard EPSG:4326 range::

raster = raster.set_origin(x=raster.origin[0] + 360)

Source code in src/rastr/raster.py

def set_origin(self, *, x: float | None = None, y: float | None = None) -> Self:
    """Set the transform origin without modifying pixel data.

    The origin is the corner of the first pixel in the raster grid —
    the translation component (c, f) of the affine transform. For
    north-up rasters this corresponds to (xmin, ymax); for south-up
    rasters it corresponds to (xmin, ymin).

    Unspecified axes retain their current value. If neither axis is
    provided, returns an unchanged copy.

    Args:
        x: New x-coordinate of the grid origin. If None, keeps current.
        y: New y-coordinate of the grid origin. If None, keeps current.

    Returns:
        A new Raster with the updated transform origin and unchanged array data.

    Example:
        Normalize a raster with non-standard longitude (e.g., -200° to -170°)
        to standard EPSG:4326 range::

            raster = raster.set_origin(x=raster.origin[0] + 360)
    """
    t = self.raster_meta.transform
    new_x = x if x is not None else t.c
    new_y = y if y is not None else t.f
    new_transform = Affine(t.a, t.b, new_x, t.d, t.e, new_y)
    new_meta = RasterMeta(
        crs=self.raster_meta.crs,
        transform=new_transform,
    )
    return self.__class__(arr=self.arr, raster_meta=new_meta)

sobel()

Compute the Sobel gradient magnitude of the raster.

This is effectively a discrete differentiation operator, computing an approximation of the magnitude of the gradient of the image intensity function.

Borders are treated using half-sample symmetric sampling, i.e. repeating the border values. Be aware that this can lead to edge artifacts and under-estimate the gradient along the border pixels.

Returns:

Type	Description
Self	New raster containing the gradient magnitude in units of raster cell units
Self	per unit distance (e.g. per meter).

Source code in src/rastr/raster.py

def sobel(self) -> Self:
    """Compute the Sobel gradient magnitude of the raster.

    This is effectively a discrete differentiation operator, computing an
    approximation of the magnitude of the gradient of the image intensity function.

    Borders are treated using half-sample symmetric sampling, i.e. repeating the
    border values. Be aware that this can lead to edge artifacts and under-estimate
    the gradient along the border pixels.

    Returns:
        New raster containing the gradient magnitude in units of raster cell units
        per unit distance (e.g. per meter).
    """
    if not self.has_square_cells:
        msg = "Sobel filter currently only supports square rasters."
        raise NotImplementedError(msg)

    from skimage import filters

    new_raster = self.model_copy()
    # Scale by cell size to convert to per unit distance
    new_raster.arr = filters.sobel(self.arr) / self.square_cell_size
    return new_raster

std()

Get the standard deviation of values in the raster, ignoring NaN values.

Returns:

Type	Description
float	The standard deviation of the raster. Returns NaN if all values are NaN.

Source code in src/rastr/raster.py

def std(self) -> float:
    """Get the standard deviation of values in the raster, ignoring NaN values.

    Returns:
        The standard deviation of the raster. Returns NaN if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nanstd(self.arr))

sum()

Get the sum of all values in the raster, ignoring NaN values.

Returns:

Type	Description
float	The sum of all values in the raster. Returns zero if all values are NaN.

Source code in src/rastr/raster.py

def sum(self) -> float:
    """Get the sum of all values in the raster, ignoring NaN values.

    Returns:
        The sum of all values in the raster. Returns zero if all values are NaN.
    """
    with suppress_slice_warning():
        return float(np.nansum(self.arr))

taper_border(width, *, limit=0.0)

Taper values to a limiting value around the border of the raster.

By default, the borders are tapered to zero, but this can be changed via the limit parameter.

This keeps the raster size the same, overwriting values in the border area. To instead grow the raster, consider using pad() followed by taper_border().

The tapering is linear from the cell centres around the border of the raster, so the value at the edge of the raster will be equal to limit.

Parameters:

Name	Type	Description	Default
width	float	The width of the taper, in the same units as the raster CRS (e.g. meters). This defines how far from the edge the tapering starts.	required
limit	float	The limiting value to taper to at the edges. Default is zero.	0.0

Source code in src/rastr/raster.py

def taper_border(self, width: float, *, limit: float = 0.0) -> Self:
    """Taper values to a limiting value around the border of the raster.

    By default, the borders are tapered to zero, but this can be changed via the
    `limit` parameter.

    This keeps the raster size the same, overwriting values in the border area.
    To instead grow the raster, consider using `pad()` followed by `taper_border()`.

    The tapering is linear from the cell centres around the border of the raster,
    so the value at the edge of the raster will be equal to `limit`.

    Args:
        width: The width of the taper, in the same units as the raster CRS
               (e.g. meters). This defines how far from the edge the tapering
               starts.
        limit: The limiting value to taper to at the edges. Default is zero.
    """

    cell_width, cell_height = self.raster_meta.cell_size
    width_in_cols = width / cell_width
    width_in_rows = width / cell_height

    # Calculate the distance from the edge in cell units
    arr_height, arr_width = self.arr.shape
    y_indices, x_indices = np.indices((int(arr_height), int(arr_width)))
    dist_from_left = x_indices
    dist_from_right = arr_width - 1 - x_indices
    dist_from_top = y_indices
    dist_from_bottom = arr_height - 1 - y_indices
    dist_from_edge_cols = np.minimum(dist_from_left, dist_from_right)
    dist_from_edge_rows = np.minimum(dist_from_top, dist_from_bottom)

    # Distance from cell centres to nearest border edge in CRS units
    dist_from_edge = np.minimum(
        dist_from_edge_cols * cell_width,
        dist_from_edge_rows * cell_height,
    )

    # Maintain edge-band behavior by rounding up separately per axis
    within_x_band = dist_from_edge_cols < np.ceil(width_in_cols)
    within_y_band = dist_from_edge_rows < np.ceil(width_in_rows)
    mask = within_x_band | within_y_band
    masked_dist_arr = np.where(mask, dist_from_edge, np.nan)
    masked_arr = np.where(mask, self.arr, np.nan)

    # Calculate the tapering factor based on the distance from the edge
    taper_factor = np.clip(masked_dist_arr / width, 0.0, 1.0)
    tapered_values = limit + (masked_arr - limit) * taper_factor

    # Create the new raster array
    new_arr = self.arr.copy()
    new_arr[mask] = tapered_values[mask]
    new_raster = self.model_copy()
    new_raster.arr = new_arr

    return new_raster

to_bounds(bounds, *, strategy='underflow')

Resize the raster to the specified bounds, cropping or padding as needed.

This method behaves like crop() when the bounds are within the current raster, but also allows expanding the raster by padding with NaN-valued cells when the bounds extend beyond the current raster boundaries.

Parameters:

Name	Type	Description	Default
bounds	tuple[float, float, float, float] | Bounds | ArrayLike	A tuple of (minx, miny, maxx, maxy) defining the target bounds.	required
strategy	Literal['underflow', 'overflow']	The strategy to use when determining cell inclusion at boundaries. 'underflow' will only include cells fully within the bounds, 'overflow' will include cells that intersect the bounds.	'underflow'

Returns:

Type	Description
Self	A new Raster instance resized to the specified bounds.

Source code in src/rastr/raster.py

def to_bounds(
    self,
    bounds: tuple[float, float, float, float] | Bounds | ArrayLike,
    *,
    strategy: Literal["underflow", "overflow"] = "underflow",
) -> Self:
    """Resize the raster to the specified bounds, cropping or padding as needed.

    This method behaves like `crop()` when the bounds are within the current raster,
    but also allows expanding the raster by padding with NaN-valued cells when the
    bounds extend beyond the current raster boundaries.

    Args:
        bounds: A tuple of (minx, miny, maxx, maxy) defining the target bounds.
        strategy:
            The strategy to use when determining cell inclusion at boundaries.
            'underflow' will only include cells fully within the bounds,
            'overflow' will include cells that intersect the bounds.

    Returns:
        A new Raster instance resized to the specified bounds.
    """
    bounds = np.asarray(bounds)
    if len(bounds) != 4:
        msg = (
            f"bounds must be a sequence of length 4 (minx, miny, maxx, maxy); "
            f"got length {len(bounds)}"
        )
        raise ValueError(msg)

    if strategy not in ("underflow", "overflow"):
        msg = f"Unsupported strategy: {strategy}"
        raise NotImplementedError(msg)

    target_minx, target_miny, target_maxx, target_maxy = bounds
    cell_width, cell_height = self.raster_meta.cell_size
    half_cell_width = cell_width / 2
    half_cell_height = cell_height / 2

    sign = 1 if strategy == "underflow" else -1
    x_range = (
        target_minx + sign * half_cell_width,
        target_maxx - sign * half_cell_width,
    )
    y_range = (
        target_miny + sign * half_cell_height,
        target_maxy - sign * half_cell_height,
    )

    if x_range[1] < x_range[0] or y_range[1] < y_range[0]:
        msg = "No cells within the specified bounds."
        raise ValueError(msg)

    # Get current cell coordinates
    current_x_coords = self.cell_x_coords
    current_y_coords = self.cell_y_coords

    # Determine the grid alignment based on the first cell center
    first_x = current_x_coords[0]
    first_y = current_y_coords[0]

    # Calculate indices for target grid (aligned to original grid)
    # x-direction is always ascending
    x_start = np.ceil((x_range[0] - first_x) / cell_width)
    x_end = np.floor((x_range[1] - first_x) / cell_width)

    # y-direction may be ascending or descending depending on transform
    # Determine direction from current_y_coords
    y_ascending = (
        len(current_y_coords) > 1 and current_y_coords[1] > current_y_coords[0]
    )

    if y_ascending:
        y_start = np.ceil((y_range[0] - first_y) / cell_height)
        y_end = np.floor((y_range[1] - first_y) / cell_height)
        y_indices = np.arange(y_start, y_end + 1, dtype=int)
        target_y_coords = first_y + y_indices * cell_height
    else:
        y_start = np.ceil((first_y - y_range[1]) / cell_height)
        y_end = np.floor((first_y - y_range[0]) / cell_height)
        y_indices = np.arange(y_start, y_end + 1, dtype=int)
        target_y_coords = first_y - y_indices * cell_height

    # Generate target coordinates
    x_indices = np.arange(x_start, x_end + 1, dtype=int)

    if len(x_indices) == 0 or len(y_indices) == 0:
        msg = "No cells within the specified bounds."
        raise ValueError(msg)

    target_x_coords = first_x + x_indices * cell_width

    output_shape = (len(target_y_coords), len(target_x_coords))
    output_arr = np.full(output_shape, np.nan, dtype=float)

    # For each target coordinate, find matching indices in current raster
    # Use broadcasting to create boolean masks
    x_match_mask = (
        np.abs(target_x_coords[:, np.newaxis] - current_x_coords[np.newaxis, :])
        < COORD_MATCH_TOLERANCE
    )
    y_match_mask = (
        np.abs(target_y_coords[:, np.newaxis] - current_y_coords[np.newaxis, :])
        < COORD_MATCH_TOLERANCE
    )

    y_target_indices, y_current_indices = np.where(y_match_mask)
    x_target_indices, x_current_indices = np.where(x_match_mask)

    for ty_idx, cy_idx in zip(y_target_indices, y_current_indices, strict=False):
        for tx_idx, cx_idx in zip(
            x_target_indices, x_current_indices, strict=False
        ):
            output_arr[ty_idx, tx_idx] = self.arr[cy_idx, cx_idx]

    transform = rasterio.transform.from_bounds(
        west=target_x_coords.min() - half_cell_width,
        south=target_y_coords.min() - half_cell_height,
        east=target_x_coords.max() + half_cell_width,
        north=target_y_coords.max() + half_cell_height,
        width=output_arr.shape[1],
        height=output_arr.shape[0],
    )

    # Create the new raster
    cls = self.__class__
    new_meta = RasterMeta(crs=self.raster_meta.crs, transform=transform)
    return cls(arr=output_arr, raster_meta=new_meta)

to_clipboard()

Copy the raster cell array to the clipboard.

Source code in src/rastr/raster.py

def to_clipboard(self) -> None:
    """Copy the raster cell array to the clipboard."""
    import pandas as pd

    pd.DataFrame(self.arr).to_clipboard(index=False, header=False)

to_file(path, **kwargs)

Write the raster to a file.

Parameters:

Name	Type	Description	Default
path	Path | str	Path to output file.	required
**kwargs	Any	Additional keyword arguments to pass to rasterio.open(). If nodata is provided, NaN values in the raster will be replaced with the nodata value.	{}

Source code in src/rastr/raster.py

def to_file(self, path: Path | str, **kwargs: Any) -> None:
    """Write the raster to a file.

    Args:
        path: Path to output file.
        **kwargs: Additional keyword arguments to pass to `rasterio.open()`. If
                  `nodata` is provided, NaN values in the raster will be replaced
                  with the nodata value.
    """
    from rastr.io_ import write_raster  # noqa: PLC0415

    return write_raster(self, path=path, **kwargs)

to_rasterio_dataset()

Create a rasterio in-memory dataset from the Raster object.

Example

raster = Raster.example() with raster.to_rasterio_dataset() as dataset: ...

Source code in src/rastr/raster.py

@contextmanager
def to_rasterio_dataset(
    self,
) -> Generator[DatasetReader | BufferedDatasetWriter | DatasetWriter]:
    """Create a rasterio in-memory dataset from the Raster object.

    Example:
        >>> raster = Raster.example()
        >>> with raster.to_rasterio_dataset() as dataset:
        >>>     ...
    """
    memfile = MemoryFile()

    height, width = self.arr.shape

    try:
        with memfile.open(
            driver="GTiff",
            height=height,
            width=width,
            count=1,  # Assuming a single band; adjust as necessary
            dtype=self.arr.dtype,
            crs=self.raster_meta.crs.to_wkt(),
            transform=self.raster_meta.transform,
        ) as dataset:
            dataset.write(self.arr, 1)

        # Yield the dataset for reading
        with memfile.open() as dataset:
            yield dataset
    finally:
        memfile.close()

trim_nan()

Crop the raster by trimming away all-NaN slices at the edges.

This effectively trims the raster to the smallest bounding box that contains all of the non-NaN values. Note that this does not guarantee no NaN values at all around the edges, only that there won't be entire edges which are all-NaN.

Consider using .extrapolate() for further cleanup of NaN values.

Source code in src/rastr/raster.py

def trim_nan(self) -> Self:
    """Crop the raster by trimming away all-NaN slices at the edges.

    This effectively trims the raster to the smallest bounding box that contains all
    of the non-NaN values. Note that this does not guarantee no NaN values at all
    around the edges, only that there won't be entire edges which are all-NaN.

    Consider using `.extrapolate()` for further cleanup of NaN values.
    """
    return self._trim_value(value_mask=np.isnan(self.arr), value_name="NaN")

trim_zeros()

Crop the raster by trimming away all-zero slices at the edges.

This effectively trims the raster to the smallest bounding box that contains all of the non-zero values. Note that this does not guarantee no zero values at all around the edges, only that there won't be entire edges which are all-zero.

Source code in src/rastr/raster.py

def trim_zeros(self) -> Self:
    """Crop the raster by trimming away all-zero slices at the edges.

    This effectively trims the raster to the smallest bounding box that contains all
    of the non-zero values. Note that this does not guarantee no zero values at all
    around the edges, only that there won't be entire edges which are all-zero.
    """
    return self._trim_value(value_mask=(self.arr == 0), value_name="zero")

unique()

Get the unique cell values in the raster, including NaN.

Returns:

Type	Description
NDArray	Array of unique values in the raster.

Source code in src/rastr/raster.py

def unique(self) -> NDArray:
    """Get the unique cell values in the raster, including NaN.

    Returns:
        Array of unique values in the raster.
    """
    return np.unique(self.arr.flatten())

RasterCellArrayShapeError

Bases: ValueError

Custom error for invalid raster cell array shapes.

Source code in src/rastr/raster.py

class RasterCellArrayShapeError(ValueError):
    """Custom error for invalid raster cell array shapes."""

suppress_slice_warning()

Context manager to suppress all-NaN slice warnings from NumPy operations.

An all-NaN slice is a row or a column of an array where every value is NaN. This is common for rasters around the edges, and is almost always a false positive.

Note that even .trim_nan() won't even guarantee there aren't NaN slices, since that method only applies around the edges of a raster, whereas there might be NaN slices in between non-NaN data in the middle.

Source code in src/rastr/raster.py

@contextmanager
def suppress_slice_warning() -> Generator[None, None, None]:
    """Context manager to suppress all-NaN slice warnings from NumPy operations.

    An all-NaN slice is a row or a column of an array where every value is NaN.
    This is common for rasters around the edges, and is almost always a false
    positive.

    Note that even .trim_nan() won't even guarantee there aren't NaN slices,
    since that method only applies around the edges of a raster, whereas there
    might be NaN slices in between non-NaN data in the middle.
    """
    with warnings.catch_warnings():
        warnings.filterwarnings(
            "ignore",
            message="All-NaN slice encountered",
            category=RuntimeWarning,
        )
        yield