Geo — Raster

Raster substrate classes and utilities backed by NumPy arrays.

dissmodel.geo.raster.backend

dissmodel/geo/raster/backend.py

Vectorized engine for cellular automata on raster grids (NumPy 2D arrays).

Responsibility

Provide generic spatial operations (shift, dilate, focal_sum, snapshot) with no domain knowledge — no land-use classes, no CRS, no I/O, no project-specific constants.

Domain models (FloodRasterModel, MangroveRasterModel, …) import RasterBackend and operate on named arrays stored in self.arrays.

Temporal variables

Arrays may be static (y, x) or temporal (time, y, x). Static variables are stored without a time axis and behave exactly as before. Temporal variables are stored with an explicit time coordinate array in self.time_coords.

# static — backward compatible
b.set("slope", slope_arr)
b.get("slope")                   # → (y, x)

# temporal
b.set("dist_roads", roads_arr, time=np.array([2000, 2014, 2020]))
b.get("dist_roads")              # → (time, y, x) full series
b.get("dist_roads", time=2014)   # → (y, x) slice at 2014

CA models always call get(name, time=step) or get(name) for static vars — they never see the time dimension directly.

Minimal example

from dissmodel.geo.raster.backend import RasterBackend, DIRS_MOORE

b = RasterBackend(shape=(100, 100))
b.set("state", np.zeros((100, 100), dtype=np.int8))

state    = b.get("state").copy()          # equivalent to cell.past[attr]
contact  = b.neighbor_contact(state == 1)
for dr, dc in DIRS_MOORE:
    neighbour = RasterBackend.shift2d(state, dr, dc)
    ...
b.arrays["state"] = new_state

DIRS_MOORE = [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)] module-attribute

DIRS_VON_NEUMANN = [(-1, 0), (0, -1), (0, 1), (1, 0)] module-attribute

RasterBackend

Storage and vectorized operations for 2D raster grids.

Replaces TerraME's forEachCell / forEachNeighbor with pure NumPy operations. The backend is shared across multiple models running in the same Environment — each model reads and writes named arrays every step.

Arrays

Stored in self.arrays as np.ndarray of shape (rows, cols) for static variables, or (time, rows, cols) for temporal variables. No names are reserved — domain models define their own ("uso", "alt", "solo", "state", "dist_roads", …).

Time coordinates for temporal variables are stored in self.time_coords as a parallel dict mapping variable name → 1D np.ndarray of time values (int or str), sorted ascending and matching the array's first dimension.

Time lookup rule

get(name, time=t) selects the slice with np.searchsorted and a clamped index — a ceiling lookup: an exact match returns that slice; a t between two coordinates returns the next (later) slice; a t outside the axis clamps to the first/last slice and never raises. Re-setting a temporal variable with a 2D array and time=None demotes it to static (its time axis is removed).

Parameters:

Name	Type	Description	Default
shape	tuple[int, int]	Grid shape as (rows, cols).	required
nodata_value	float | int | None	Sentinel value used to mark cells outside the study extent. When provided, nodata_mask derives the extent mask automatically. Default: None.	None

Examples:

>>> b = RasterBackend(shape=(10, 10))
>>> b.set("state", np.zeros((10, 10), dtype=np.int8))
>>> b.get("state").shape
(10, 10)

>>> b = RasterBackend(shape=(10, 10), nodata_value=-1)
>>> b.nodata_mask   # True = valid cell, False = outside extent

>>> # temporal variable
>>> roads_3d = np.zeros((3, 10, 10))
>>> b.set("dist_roads", roads_3d, time=np.array([2000, 2014, 2020]))
>>> b.get("dist_roads", time=2014).shape
(10, 10)
>>> b.get("dist_roads").shape
(3, 10, 10)

Source code in dissmodel/geo/raster/backend.py

class RasterBackend:
    """
    Storage and vectorized operations for 2D raster grids.

    Replaces TerraME's ``forEachCell`` / ``forEachNeighbor`` with pure NumPy
    operations. The backend is shared across multiple models running in the
    same ``Environment`` — each model reads and writes named arrays every step.

    Arrays
    ------
    Stored in ``self.arrays`` as ``np.ndarray`` of shape ``(rows, cols)``
    for static variables, or ``(time, rows, cols)`` for temporal variables.
    No names are reserved — domain models define their own
    (``"uso"``, ``"alt"``, ``"solo"``, ``"state"``, ``"dist_roads"``, …).

    Time coordinates for temporal variables are stored in ``self.time_coords``
    as a parallel dict mapping variable name → 1D ``np.ndarray`` of time values
    (int or str), sorted ascending and matching the array's first dimension.

    Time lookup rule
    ----------------
    ``get(name, time=t)`` selects the slice with ``np.searchsorted`` and a
    clamped index — a *ceiling* lookup: an exact match returns that slice;
    a ``t`` between two coordinates returns the next (later) slice; a ``t``
    outside the axis clamps to the first/last slice and never raises.
    Re-setting a temporal variable with a 2D array and ``time=None`` demotes
    it to static (its time axis is removed).

    Parameters
    ----------
    shape : tuple[int, int]
        Grid shape as ``(rows, cols)``.
    nodata_value : float | int | None
        Sentinel value used to mark cells outside the study extent.
        When provided, ``nodata_mask`` derives the extent mask automatically.
        Default: ``None``.

    Examples
    --------
    >>> b = RasterBackend(shape=(10, 10))
    >>> b.set("state", np.zeros((10, 10), dtype=np.int8))
    >>> b.get("state").shape
    (10, 10)

    >>> b = RasterBackend(shape=(10, 10), nodata_value=-1)
    >>> b.nodata_mask   # True = valid cell, False = outside extent

    >>> # temporal variable
    >>> roads_3d = np.zeros((3, 10, 10))
    >>> b.set("dist_roads", roads_3d, time=np.array([2000, 2014, 2020]))
    >>> b.get("dist_roads", time=2014).shape
    (10, 10)
    >>> b.get("dist_roads").shape
    (3, 10, 10)
    """

    def __init__(
        self,
        shape: tuple[int, int],
        nodata_value: float | int | None = None,
        transform: Any = None,
        crs: Any = None,
    ) -> None:
        self.shape        = shape
        self.arrays: dict[str, np.ndarray] = {}
        self.time_coords: dict[str, np.ndarray] = {}  # name → 1D time axis
        self.nodata_value = nodata_value

        self.transform    = transform
        self.crs          = crs

    # ── temporal helpers ──────────────────────────────────────────────────────

    def is_temporal(self, name: str) -> bool:
        """Return ``True`` if ``name`` has an associated time axis."""
        return name in self.time_coords

    def time_axis(self, name: str) -> np.ndarray | None:
        """Return the time coordinate array for ``name``, or ``None``."""
        return self.time_coords.get(name)

    # ── extent mask ───────────────────────────────────────────────────────────

    @property
    def nodata_mask(self) -> np.ndarray | None:
        """
        Boolean mask: ``True`` = valid cell, ``False`` = outside extent / nodata.

        Derived in priority order:
        1. ``arrays["mask"]``  — explicit mask band.
        2. ``nodata_value``    — applied over the first available static array.
        3. ``None``            — no information.
        """
        if "mask" in self.arrays:
            return self.arrays["mask"] != 0

        if self.nodata_value is not None and self.arrays:
            # use first static (2D) array for the mask
            for arr in self.arrays.values():
                if arr.ndim == 2:
                    return arr != self.nodata_value

        return None

    # ── read / write ──────────────────────────────────────────────────────────

    def set(
        self,
        name: str,
        array: np.ndarray,
        time: np.ndarray | list | None = None,
    ) -> None:
        """
        Store ``array`` under ``name``.

        Parameters
        ----------
        name : str
            Variable name.
        array : np.ndarray
            Shape ``(y, x)`` for static variables, ``(time, y, x)`` for temporal.
        time : array-like | None
            1D sequence of time coordinate values (int or str) matching the
            first dimension of ``array``. If provided, the variable is marked
            as temporal and ``get(name, time=t)`` will return a 2D slice.
            Must be ``None`` for static (2D) arrays.

        Raises
        ------
        ValueError
            If ``time`` length does not match the first dimension of ``array``.
        """
        arr = np.asarray(array).copy()

        if time is not None:
            t = np.asarray(time)
            if arr.ndim != 3:
                raise ValueError(
                    f"Expected 3D array (time, y, x) when time is given, "
                    f"got shape {arr.shape}"
                )
            if len(t) != arr.shape[0]:
                raise ValueError(
                    f"time length ({len(t)}) must match array first dim ({arr.shape[0]})"
                )
            self.time_coords[name] = t
        else:
            # remove any stale time axis when overwriting with static array
            self.time_coords.pop(name, None)

        self.arrays[name] = arr

    def get(
        self,
        name: str,
        time: int | str | None = None,
    ) -> np.ndarray:
        """
        Return array for ``name``.

        Behaviour
        ---------
        - Static variable (no time axis): always returns ``(y, x)``.
          The ``time`` argument is silently ignored, so CA models can call
          ``get(name, time=step)`` uniformly without checking variable type.
        - Temporal variable + ``time=None``: returns full ``(time, y, x)`` series.
        - Temporal variable + ``time=t``: returns the ``(y, x)`` slice selected
          by a ceiling lookup — exact match returns that slice, a ``t`` between
          coordinates returns the next (later) slice, and out-of-range values
          clamp to the first/last slice without raising.

        Parameters
        ----------
        name : str
        time : int | str | None
            Time value to select. Uses ``np.searchsorted`` with a clamped
            index (ceiling rule).

        Raises
        ------
        KeyError
            If ``name`` is not in ``self.arrays``.
        """
        arr = self.arrays[name]

        if time is None or name not in self.time_coords:
            return arr

        idx = int(np.searchsorted(self.time_coords[name], time))
        # clamp to valid range
        idx = max(0, min(idx, arr.shape[0] - 1))
        return arr[idx]

    def snapshot(self) -> dict[str, np.ndarray]:
        """
        Return a deep copy of all arrays — equivalent to TerraME's ``.past`` mechanism.

        For temporal variables the full ``(time, y, x)`` array is copied.
        Use ``get(name, time=t)`` on the backend directly to snapshot a single
        time slice.

        Returns
        -------
        dict[str, np.ndarray]
        """
        return {k: v.copy() for k, v in self.arrays.items()}

    def rename_band(self, old: str, new: str) -> None:
        """
        Rename an array in-place. No-op if ``old`` does not exist.
        Time coordinates are renamed alongside the array.
        """
        if old in self.arrays:
            self.arrays[new] = self.arrays.pop(old)
            if old in self.time_coords:
                self.time_coords[new] = self.time_coords.pop(old)

    # ── xarray interoperability ───────────────────────────────────────────────

    def to_xarray(self, time: int | None = None):
        """
        Convert the backend to an ``xr.Dataset``.

        Static variables become ``DataArray(y, x)``.
        Temporal variables become ``DataArray(time, y, x)`` with explicit
        time coordinates from ``self.time_coords``.

        If ``time`` is given (simulation step), a scalar ``time`` coordinate
        is added to static variables — useful when assembling multi-step outputs.

        Parameters
        ----------
        time : int | None
            Optional simulation step to attach as a scalar coordinate
            (applies to static variables only).

        Returns
        -------
        xr.Dataset
        """
        try:
            import xarray as xr
        except ImportError:
            raise ImportError(
                "xarray is required for RasterBackend.to_xarray(). "
                "Install it with: pip install xarray"
            )

        rows, cols = self.shape

        if self.transform is not None:
            try:
                xs = np.array([self.transform.c + (c + 0.5) * self.transform.a
                                for c in range(cols)])
                ys = np.array([self.transform.f + (r + 0.5) * self.transform.e
                                for r in range(rows)])
            except AttributeError:
                xs = np.arange(cols, dtype=float)
                ys = np.arange(rows, dtype=float)
        else:
            xs = np.arange(cols, dtype=float)
            ys = np.arange(rows, dtype=float)

        base_coords: dict = {"y": ys, "x": xs}

        # CRS as spatial_ref coordinate (CF / rioxarray convention)
        if self.crs is not None:
            try:
                from pyproj import CRS as ProjCRS
                import xarray as xr
                crs_obj = ProjCRS.from_user_input(self.crs)
                base_coords["spatial_ref"] = xr.DataArray(
                    0,
                    attrs={
                        "crs_wkt":      crs_obj.to_wkt(),
                        "grid_mapping": "spatial_ref",
                    },
                )
            except Exception:
                pass

        data_vars = {}
        for name, arr in self.arrays.items():
            attrs: dict = {}
            if self.nodata_value is not None:
                attrs["_FillValue"] = self.nodata_value
            if self.crs is not None and "spatial_ref" in base_coords:
                attrs["grid_mapping"] = "spatial_ref"

            if name in self.time_coords:
                # temporal variable — emit (time, y, x)
                data_vars[name] = xr.DataArray(
                    arr.copy(),
                    dims=["time", "y", "x"],
                    coords={
                        "time": self.time_coords[name],
                        "y": ys,
                        "x": xs,
                    },
                    attrs=attrs,
                )
            else:
                # static variable — emit (y, x)
                # holds both coordinate arrays and the optional scalar time
                coords: dict[str, Any] = {"y": ys, "x": xs}
                if time is not None:
                    coords["time"] = time
                data_vars[name] = xr.DataArray(
                    arr.copy(),
                    dims=["y", "x"],
                    coords=coords,
                    attrs=attrs,
                )

        ds = xr.Dataset(data_vars, coords=base_coords)
        ds.attrs["Conventions"] = "CF-1.8"
        return ds

    @classmethod
    def from_xarray(cls, ds, nodata_value: float | int | None = None) -> "RasterBackend":
        """
        Build a ``RasterBackend`` from an ``xr.Dataset`` or ``xr.DataArray``.

        Variables with dimensions ``(y, x)`` are imported as static arrays.
        Variables with dimensions ``(time, y, x)`` are imported as temporal
        arrays with their time coordinates stored in ``self.time_coords``.

        Parameters
        ----------
        ds : xr.Dataset | xr.DataArray
        nodata_value : float | int | None

        Returns
        -------
        RasterBackend

        Raises
        ------
        ValueError
            If ``ds`` contains no spatial variables.
        """
        try:
            import xarray as xr
        except ImportError:
            raise ImportError(
                "xarray is required for RasterBackend.from_xarray(). "
                "Install it with: pip install xarray"
            )

        if isinstance(ds, xr.DataArray):
            name = ds.name or "data"
            ds = ds.to_dataset(name=name)

        spatial_vars = {
            name: var
            for name, var in ds.data_vars.items()
            if set(var.dims) >= {"y", "x"} and var.ndim in (2, 3)
        }

        if not spatial_vars:
            raise ValueError(
                "No spatial (y, x) or (time, y, x) variables found in Dataset."
            )

        # infer shape from first 2D variable, or spatial slice of first 3D
        shape_2d = None
        for var in spatial_vars.values():
            if var.ndim == 2:
                yi = var.dims.index("y")
                xi = var.dims.index("x")
                shape_2d = (var.shape[yi], var.shape[xi])
                break
        if shape_2d is None:
            var = next(iter(spatial_vars.values()))
            yi = var.dims.index("y")
            xi = var.dims.index("x")
            shape_2d = (var.shape[yi], var.shape[xi])

        rows, cols = shape_2d

        # recover transform
        transform = None
        try:
            import rasterio.transform
            ys = ds.coords["y"].values
            xs = ds.coords["x"].values
            if len(ys) >= 2 and len(xs) >= 2:
                res_y = float(ys[1] - ys[0])
                res_x = float(xs[1] - xs[0])
                origin_x = float(xs[0]) - res_x / 2
                origin_y = float(ys[0]) - res_y / 2
                transform = rasterio.transform.from_origin(
                    origin_x, origin_y - res_y * (rows - 1), res_x, abs(res_y)
                ) if res_y < 0 else rasterio.transform.Affine(
                    res_x, 0, origin_x, 0, res_y, origin_y
                )
        except Exception:
            pass

        # recover CRS
        crs = None
        if "spatial_ref" in ds.coords:
            try:
                from pyproj import CRS as ProjCRS
                wkt = ds.coords["spatial_ref"].attrs.get("crs_wkt", "")
                if wkt:
                    crs = ProjCRS.from_wkt(wkt)
            except Exception:
                pass

        backend = cls(
            shape=(rows, cols),
            nodata_value=nodata_value,
            transform=transform,
            crs=crs,
        )

        for name, var in spatial_vars.items():
            if var.ndim == 3 and "time" in var.dims:
                # temporal variable
                arr = var.transpose("time", "y", "x").values
                t   = ds.coords["time"].values
                backend.arrays[name]      = arr.copy()
                backend.time_coords[name] = t
            else:
                # static variable
                arr = var.transpose("y", "x").values
                backend.arrays[name] = arr.copy()

        return backend

    # ── spatial operations ────────────────────────────────────────────────────

    @staticmethod
    def shift2d(arr: np.ndarray, dr: int, dc: int) -> np.ndarray:
        """
        Shift ``arr`` by ``(dr, dc)`` rows/columns without wrap-around.
        Edges are filled with zero.

        Parameters
        ----------
        arr : np.ndarray  shape (y, x)
        dr : int
        dc : int

        Returns
        -------
        np.ndarray
        """
        rows, cols = arr.shape
        out = np.zeros_like(arr)
        rs  = slice(max(0, -dr), min(rows, rows - dr))
        rd  = slice(max(0,  dr), min(rows, rows + dr))
        cs_ = slice(max(0, -dc), min(cols, cols - dc))
        cd  = slice(max(0,  dc), min(cols, cols + dc))
        out[rd, cd] = arr[rs, cs_]
        return out

    def band_names(self) -> list[str]:
        """Return the names of all arrays currently stored in the backend."""
        return list(self.arrays.keys())

    def temporal_band_names(self) -> list[str]:
        """Return the names of temporal (time, y, x) arrays."""
        return list(self.time_coords.keys())

    def static_band_names(self) -> list[str]:
        """Return the names of static (y, x) arrays."""
        return [k for k in self.arrays if k not in self.time_coords]

    @staticmethod
    def neighbor_contact(
        condition: np.ndarray,
        neighborhood: list[tuple[int, int]] | None = None,
    ) -> np.ndarray:
        """
        Return a boolean mask where each cell has at least one neighbour
        satisfying ``condition``.

        Parameters
        ----------
        condition : np.ndarray  shape (y, x)
        neighborhood : list[tuple[int, int]] | None

        Returns
        -------
        np.ndarray  bool
        """
        if neighborhood is None:
            return binary_dilation(condition.astype(bool), structure=np.ones((3, 3)))
        result = np.zeros_like(condition, dtype=bool)
        for dr, dc in neighborhood:
            result |= RasterBackend.shift2d(condition.astype(np.int8), dr, dc) > 0
        return result

    def focal_sum(
        self,
        name: str,
        neighborhood: list[tuple[int, int]] = DIRS_MOORE,
    ) -> np.ndarray:
        """
        Focal sum across neighbours for a static (y, x) array.

        Parameters
        ----------
        name : str
        neighborhood : list[tuple[int, int]]

        Returns
        -------
        np.ndarray
        """
        arr    = self.arrays[name]
        if arr.ndim != 2:
            raise ValueError(
                f"focal_sum requires a static 2D array. "
                f"'{name}' has shape {arr.shape}. "
                f"Use get('{name}', time=t) to select a slice first."
            )
        result = np.zeros_like(arr, dtype=float)
        for dr, dc in neighborhood:
            result += self.shift2d(arr, dr, dc)
        return result

    def focal_sum_mask(
        self,
        mask: np.ndarray,
        neighborhood: list[tuple[int, int]] = DIRS_MOORE,
    ) -> np.ndarray:
        """
        Count neighbours where ``mask`` is ``True``.

        Parameters
        ----------
        mask : np.ndarray  shape (y, x)
        neighborhood : list[tuple[int, int]]

        Returns
        -------
        np.ndarray  int
        """
        result: np.ndarray = np.zeros(self.shape, dtype=int)
        m: np.ndarray = mask.astype(np.int8)
        for dr, dc in neighborhood:
            result += self.shift2d(m, dr, dc)
        return result

    # ── utilities ─────────────────────────────────────────────────────────────

    def __repr__(self) -> str:
        static   = [f"{k}:{v.dtype}{list(v.shape)}"
                    for k, v in self.arrays.items() if k not in self.time_coords]
        temporal = [f"{k}:{v.dtype}{list(v.shape)}@{list(self.time_coords[k])}"
                    for k, v in self.arrays.items() if k in self.time_coords]
        parts = static + temporal
        return f"RasterBackend(shape={self.shape}, arrays=[{', '.join(parts)}])"

nodata_mask property

Boolean mask: True = valid cell, False = outside extent / nodata.

Derived in priority order: 1. arrays["mask"] — explicit mask band. 2. nodata_value — applied over the first available static array. 3. None — no information.

band_names()

Return the names of all arrays currently stored in the backend.

Source code in dissmodel/geo/raster/backend.py

def band_names(self) -> list[str]:
    """Return the names of all arrays currently stored in the backend."""
    return list(self.arrays.keys())

focal_sum(name, neighborhood=DIRS_MOORE)

Focal sum across neighbours for a static (y, x) array.

Parameters:

Name	Type	Description	Default
name	str		required
neighborhood	list[tuple[int, int]]		DIRS_MOORE

Returns:

Type	Description
ndarray

Source code in dissmodel/geo/raster/backend.py

def focal_sum(
    self,
    name: str,
    neighborhood: list[tuple[int, int]] = DIRS_MOORE,
) -> np.ndarray:
    """
    Focal sum across neighbours for a static (y, x) array.

    Parameters
    ----------
    name : str
    neighborhood : list[tuple[int, int]]

    Returns
    -------
    np.ndarray
    """
    arr    = self.arrays[name]
    if arr.ndim != 2:
        raise ValueError(
            f"focal_sum requires a static 2D array. "
            f"'{name}' has shape {arr.shape}. "
            f"Use get('{name}', time=t) to select a slice first."
        )
    result = np.zeros_like(arr, dtype=float)
    for dr, dc in neighborhood:
        result += self.shift2d(arr, dr, dc)
    return result

focal_sum_mask(mask, neighborhood=DIRS_MOORE)

Count neighbours where mask is True.

Parameters:

Name	Type	Description	Default
mask	np.ndarray shape (y, x)		required
neighborhood	list[tuple[int, int]]		DIRS_MOORE

Returns:

Type	Description
np.ndarray int

Source code in dissmodel/geo/raster/backend.py

def focal_sum_mask(
    self,
    mask: np.ndarray,
    neighborhood: list[tuple[int, int]] = DIRS_MOORE,
) -> np.ndarray:
    """
    Count neighbours where ``mask`` is ``True``.

    Parameters
    ----------
    mask : np.ndarray  shape (y, x)
    neighborhood : list[tuple[int, int]]

    Returns
    -------
    np.ndarray  int
    """
    result: np.ndarray = np.zeros(self.shape, dtype=int)
    m: np.ndarray = mask.astype(np.int8)
    for dr, dc in neighborhood:
        result += self.shift2d(m, dr, dc)
    return result

from_xarray(ds, nodata_value=None) classmethod

Build a RasterBackend from an xr.Dataset or xr.DataArray.

Variables with dimensions (y, x) are imported as static arrays. Variables with dimensions (time, y, x) are imported as temporal arrays with their time coordinates stored in self.time_coords.

Parameters:

Name	Type	Description	Default
ds	Dataset | DataArray		required
nodata_value	float | int | None		None

Returns:

Type	Description
RasterBackend

Raises:

Type	Description
ValueError	If ds contains no spatial variables.

Source code in dissmodel/geo/raster/backend.py

@classmethod
def from_xarray(cls, ds, nodata_value: float | int | None = None) -> "RasterBackend":
    """
    Build a ``RasterBackend`` from an ``xr.Dataset`` or ``xr.DataArray``.

    Variables with dimensions ``(y, x)`` are imported as static arrays.
    Variables with dimensions ``(time, y, x)`` are imported as temporal
    arrays with their time coordinates stored in ``self.time_coords``.

    Parameters
    ----------
    ds : xr.Dataset | xr.DataArray
    nodata_value : float | int | None

    Returns
    -------
    RasterBackend

    Raises
    ------
    ValueError
        If ``ds`` contains no spatial variables.
    """
    try:
        import xarray as xr
    except ImportError:
        raise ImportError(
            "xarray is required for RasterBackend.from_xarray(). "
            "Install it with: pip install xarray"
        )

    if isinstance(ds, xr.DataArray):
        name = ds.name or "data"
        ds = ds.to_dataset(name=name)

    spatial_vars = {
        name: var
        for name, var in ds.data_vars.items()
        if set(var.dims) >= {"y", "x"} and var.ndim in (2, 3)
    }

    if not spatial_vars:
        raise ValueError(
            "No spatial (y, x) or (time, y, x) variables found in Dataset."
        )

    # infer shape from first 2D variable, or spatial slice of first 3D
    shape_2d = None
    for var in spatial_vars.values():
        if var.ndim == 2:
            yi = var.dims.index("y")
            xi = var.dims.index("x")
            shape_2d = (var.shape[yi], var.shape[xi])
            break
    if shape_2d is None:
        var = next(iter(spatial_vars.values()))
        yi = var.dims.index("y")
        xi = var.dims.index("x")
        shape_2d = (var.shape[yi], var.shape[xi])

    rows, cols = shape_2d

    # recover transform
    transform = None
    try:
        import rasterio.transform
        ys = ds.coords["y"].values
        xs = ds.coords["x"].values
        if len(ys) >= 2 and len(xs) >= 2:
            res_y = float(ys[1] - ys[0])
            res_x = float(xs[1] - xs[0])
            origin_x = float(xs[0]) - res_x / 2
            origin_y = float(ys[0]) - res_y / 2
            transform = rasterio.transform.from_origin(
                origin_x, origin_y - res_y * (rows - 1), res_x, abs(res_y)
            ) if res_y < 0 else rasterio.transform.Affine(
                res_x, 0, origin_x, 0, res_y, origin_y
            )
    except Exception:
        pass

    # recover CRS
    crs = None
    if "spatial_ref" in ds.coords:
        try:
            from pyproj import CRS as ProjCRS
            wkt = ds.coords["spatial_ref"].attrs.get("crs_wkt", "")
            if wkt:
                crs = ProjCRS.from_wkt(wkt)
        except Exception:
            pass

    backend = cls(
        shape=(rows, cols),
        nodata_value=nodata_value,
        transform=transform,
        crs=crs,
    )

    for name, var in spatial_vars.items():
        if var.ndim == 3 and "time" in var.dims:
            # temporal variable
            arr = var.transpose("time", "y", "x").values
            t   = ds.coords["time"].values
            backend.arrays[name]      = arr.copy()
            backend.time_coords[name] = t
        else:
            # static variable
            arr = var.transpose("y", "x").values
            backend.arrays[name] = arr.copy()

    return backend

get(name, time=None)

Return array for name.

Behaviour

• Static variable (no time axis): always returns (y, x). The time argument is silently ignored, so CA models can call get(name, time=step) uniformly without checking variable type.
• Temporal variable + time=None: returns full (time, y, x) series.
• Temporal variable + time=t: returns the (y, x) slice selected by a ceiling lookup — exact match returns that slice, a t between coordinates returns the next (later) slice, and out-of-range values clamp to the first/last slice without raising.

Parameters:

Name	Type	Description	Default
name	str		required
time	int | str | None	Time value to select. Uses np.searchsorted with a clamped index (ceiling rule).	None

Raises:

Type	Description
KeyError	If name is not in self.arrays.

Source code in dissmodel/geo/raster/backend.py

def get(
    self,
    name: str,
    time: int | str | None = None,
) -> np.ndarray:
    """
    Return array for ``name``.

    Behaviour
    ---------
    - Static variable (no time axis): always returns ``(y, x)``.
      The ``time`` argument is silently ignored, so CA models can call
      ``get(name, time=step)`` uniformly without checking variable type.
    - Temporal variable + ``time=None``: returns full ``(time, y, x)`` series.
    - Temporal variable + ``time=t``: returns the ``(y, x)`` slice selected
      by a ceiling lookup — exact match returns that slice, a ``t`` between
      coordinates returns the next (later) slice, and out-of-range values
      clamp to the first/last slice without raising.

    Parameters
    ----------
    name : str
    time : int | str | None
        Time value to select. Uses ``np.searchsorted`` with a clamped
        index (ceiling rule).

    Raises
    ------
    KeyError
        If ``name`` is not in ``self.arrays``.
    """
    arr = self.arrays[name]

    if time is None or name not in self.time_coords:
        return arr

    idx = int(np.searchsorted(self.time_coords[name], time))
    # clamp to valid range
    idx = max(0, min(idx, arr.shape[0] - 1))
    return arr[idx]

is_temporal(name)

Return True if name has an associated time axis.

Source code in dissmodel/geo/raster/backend.py

def is_temporal(self, name: str) -> bool:
    """Return ``True`` if ``name`` has an associated time axis."""
    return name in self.time_coords

neighbor_contact(condition, neighborhood=None) staticmethod

Return a boolean mask where each cell has at least one neighbour satisfying condition.

Parameters:

Name	Type	Description	Default
condition	np.ndarray shape (y, x)		required
neighborhood	list[tuple[int, int]] | None		None

Returns:

Type	Description
np.ndarray bool

Source code in dissmodel/geo/raster/backend.py

@staticmethod
def neighbor_contact(
    condition: np.ndarray,
    neighborhood: list[tuple[int, int]] | None = None,
) -> np.ndarray:
    """
    Return a boolean mask where each cell has at least one neighbour
    satisfying ``condition``.

    Parameters
    ----------
    condition : np.ndarray  shape (y, x)
    neighborhood : list[tuple[int, int]] | None

    Returns
    -------
    np.ndarray  bool
    """
    if neighborhood is None:
        return binary_dilation(condition.astype(bool), structure=np.ones((3, 3)))
    result = np.zeros_like(condition, dtype=bool)
    for dr, dc in neighborhood:
        result |= RasterBackend.shift2d(condition.astype(np.int8), dr, dc) > 0
    return result

rename_band(old, new)

Rename an array in-place. No-op if old does not exist. Time coordinates are renamed alongside the array.

Source code in dissmodel/geo/raster/backend.py

def rename_band(self, old: str, new: str) -> None:
    """
    Rename an array in-place. No-op if ``old`` does not exist.
    Time coordinates are renamed alongside the array.
    """
    if old in self.arrays:
        self.arrays[new] = self.arrays.pop(old)
        if old in self.time_coords:
            self.time_coords[new] = self.time_coords.pop(old)

set(name, array, time=None)

Store array under name.

Parameters:

Name	Type	Description	Default
name	str	Variable name.	required
array	ndarray	Shape (y, x) for static variables, (time, y, x) for temporal.	required
time	array - like | None	1D sequence of time coordinate values (int or str) matching the first dimension of array. If provided, the variable is marked as temporal and get(name, time=t) will return a 2D slice. Must be None for static (2D) arrays.	None

Raises:

Type	Description
ValueError	If time length does not match the first dimension of array.

Source code in dissmodel/geo/raster/backend.py

def set(
    self,
    name: str,
    array: np.ndarray,
    time: np.ndarray | list | None = None,
) -> None:
    """
    Store ``array`` under ``name``.

    Parameters
    ----------
    name : str
        Variable name.
    array : np.ndarray
        Shape ``(y, x)`` for static variables, ``(time, y, x)`` for temporal.
    time : array-like | None
        1D sequence of time coordinate values (int or str) matching the
        first dimension of ``array``. If provided, the variable is marked
        as temporal and ``get(name, time=t)`` will return a 2D slice.
        Must be ``None`` for static (2D) arrays.

    Raises
    ------
    ValueError
        If ``time`` length does not match the first dimension of ``array``.
    """
    arr = np.asarray(array).copy()

    if time is not None:
        t = np.asarray(time)
        if arr.ndim != 3:
            raise ValueError(
                f"Expected 3D array (time, y, x) when time is given, "
                f"got shape {arr.shape}"
            )
        if len(t) != arr.shape[0]:
            raise ValueError(
                f"time length ({len(t)}) must match array first dim ({arr.shape[0]})"
            )
        self.time_coords[name] = t
    else:
        # remove any stale time axis when overwriting with static array
        self.time_coords.pop(name, None)

    self.arrays[name] = arr

shift2d(arr, dr, dc) staticmethod

Shift arr by (dr, dc) rows/columns without wrap-around. Edges are filled with zero.

Parameters:

Name	Type	Description	Default
arr	np.ndarray shape (y, x)		required
dr	int		required
dc	int		required

Returns:

Type	Description
ndarray

Source code in dissmodel/geo/raster/backend.py

@staticmethod
def shift2d(arr: np.ndarray, dr: int, dc: int) -> np.ndarray:
    """
    Shift ``arr`` by ``(dr, dc)`` rows/columns without wrap-around.
    Edges are filled with zero.

    Parameters
    ----------
    arr : np.ndarray  shape (y, x)
    dr : int
    dc : int

    Returns
    -------
    np.ndarray
    """
    rows, cols = arr.shape
    out = np.zeros_like(arr)
    rs  = slice(max(0, -dr), min(rows, rows - dr))
    rd  = slice(max(0,  dr), min(rows, rows + dr))
    cs_ = slice(max(0, -dc), min(cols, cols - dc))
    cd  = slice(max(0,  dc), min(cols, cols + dc))
    out[rd, cd] = arr[rs, cs_]
    return out

snapshot()

Return a deep copy of all arrays — equivalent to TerraME's .past mechanism.

For temporal variables the full (time, y, x) array is copied. Use get(name, time=t) on the backend directly to snapshot a single time slice.

Returns:

Type	Description
dict[str, ndarray]

Source code in dissmodel/geo/raster/backend.py

def snapshot(self) -> dict[str, np.ndarray]:
    """
    Return a deep copy of all arrays — equivalent to TerraME's ``.past`` mechanism.

    For temporal variables the full ``(time, y, x)`` array is copied.
    Use ``get(name, time=t)`` on the backend directly to snapshot a single
    time slice.

    Returns
    -------
    dict[str, np.ndarray]
    """
    return {k: v.copy() for k, v in self.arrays.items()}

static_band_names()

Return the names of static (y, x) arrays.

Source code in dissmodel/geo/raster/backend.py

def static_band_names(self) -> list[str]:
    """Return the names of static (y, x) arrays."""
    return [k for k in self.arrays if k not in self.time_coords]

temporal_band_names()

Return the names of temporal (time, y, x) arrays.

Source code in dissmodel/geo/raster/backend.py

def temporal_band_names(self) -> list[str]:
    """Return the names of temporal (time, y, x) arrays."""
    return list(self.time_coords.keys())

time_axis(name)

Return the time coordinate array for name, or None.

Source code in dissmodel/geo/raster/backend.py

def time_axis(self, name: str) -> np.ndarray | None:
    """Return the time coordinate array for ``name``, or ``None``."""
    return self.time_coords.get(name)

to_xarray(time=None)

Convert the backend to an xr.Dataset.

Static variables become DataArray(y, x). Temporal variables become DataArray(time, y, x) with explicit time coordinates from self.time_coords.

If time is given (simulation step), a scalar time coordinate is added to static variables — useful when assembling multi-step outputs.

Parameters:

Name	Type	Description	Default
time	int | None	Optional simulation step to attach as a scalar coordinate (applies to static variables only).	None

Returns:

Type	Description
Dataset

Source code in dissmodel/geo/raster/backend.py

def to_xarray(self, time: int | None = None):
    """
    Convert the backend to an ``xr.Dataset``.

    Static variables become ``DataArray(y, x)``.
    Temporal variables become ``DataArray(time, y, x)`` with explicit
    time coordinates from ``self.time_coords``.

    If ``time`` is given (simulation step), a scalar ``time`` coordinate
    is added to static variables — useful when assembling multi-step outputs.

    Parameters
    ----------
    time : int | None
        Optional simulation step to attach as a scalar coordinate
        (applies to static variables only).

    Returns
    -------
    xr.Dataset
    """
    try:
        import xarray as xr
    except ImportError:
        raise ImportError(
            "xarray is required for RasterBackend.to_xarray(). "
            "Install it with: pip install xarray"
        )

    rows, cols = self.shape

    if self.transform is not None:
        try:
            xs = np.array([self.transform.c + (c + 0.5) * self.transform.a
                            for c in range(cols)])
            ys = np.array([self.transform.f + (r + 0.5) * self.transform.e
                            for r in range(rows)])
        except AttributeError:
            xs = np.arange(cols, dtype=float)
            ys = np.arange(rows, dtype=float)
    else:
        xs = np.arange(cols, dtype=float)
        ys = np.arange(rows, dtype=float)

    base_coords: dict = {"y": ys, "x": xs}

    # CRS as spatial_ref coordinate (CF / rioxarray convention)
    if self.crs is not None:
        try:
            from pyproj import CRS as ProjCRS
            import xarray as xr
            crs_obj = ProjCRS.from_user_input(self.crs)
            base_coords["spatial_ref"] = xr.DataArray(
                0,
                attrs={
                    "crs_wkt":      crs_obj.to_wkt(),
                    "grid_mapping": "spatial_ref",
                },
            )
        except Exception:
            pass

    data_vars = {}
    for name, arr in self.arrays.items():
        attrs: dict = {}
        if self.nodata_value is not None:
            attrs["_FillValue"] = self.nodata_value
        if self.crs is not None and "spatial_ref" in base_coords:
            attrs["grid_mapping"] = "spatial_ref"

        if name in self.time_coords:
            # temporal variable — emit (time, y, x)
            data_vars[name] = xr.DataArray(
                arr.copy(),
                dims=["time", "y", "x"],
                coords={
                    "time": self.time_coords[name],
                    "y": ys,
                    "x": xs,
                },
                attrs=attrs,
            )
        else:
            # static variable — emit (y, x)
            # holds both coordinate arrays and the optional scalar time
            coords: dict[str, Any] = {"y": ys, "x": xs}
            if time is not None:
                coords["time"] = time
            data_vars[name] = xr.DataArray(
                arr.copy(),
                dims=["y", "x"],
                coords=coords,
                attrs=attrs,
            )

    ds = xr.Dataset(data_vars, coords=base_coords)
    ds.attrs["Conventions"] = "CF-1.8"
    return ds

dissmodel.geo.raster.raster_model

dissmodel/geo/raster/model.py

Base class for models backed by RasterBackend (NumPy 2D arrays).

Analogous to SpatialModel for the raster substrate — provides infrastructure without imposing a transition rule contract.

Class hierarchy

Model  (dissmodel.core)
  ├── SpatialModel     GeoDataFrame + Queen/Rook neighbourhood  (vector)
  └── RasterModel      RasterBackend + shift2d                  (raster)  ← this file
        ├── FloodRasterModel
        └── MangroveRasterModel

Usage

class MyRasterModel(RasterModel):
    def setup(self, backend, my_param=1.0):
        super().setup(backend)
        self.my_param = my_param

    def execute(self):
        uso = self.backend.get("uso").copy()
        ...
        self.backend.arrays["uso"] = new_uso

RasterModel

Bases: Model

Model backed by a RasterBackend.

Subclass of Model that adds raster infrastructure without imposing a transition rule contract. Can be subclassed directly by any model that operates on NumPy 2D arrays.

Parameters (setup)

backend : RasterBackend Backend shared across all models in the same Environment.

Attributes available in subclasses

backend : RasterBackend The shared array store. shape : tuple[int, int] Grid shape (rows, cols) — shortcut for self.backend.shape. shift : callable Shortcut for RasterBackend.shift2d (static method). dirs : list[tuple[int, int]] DIRS_MOORE — the 8 directions of the Moore neighbourhood.

Examples:

>>> class HeatDiffusion(RasterModel):
...     def execute(self):
...         temp = self.backend.get("temp").copy()
...         for dr, dc in self.dirs:
...             temp += 0.1 * self.shift(temp, dr, dc)
...         self.backend.arrays["temp"] = temp

Source code in dissmodel/geo/raster/raster_model.py

class RasterModel(Model):
    """
    Model backed by a RasterBackend.

    Subclass of ``Model`` that adds raster infrastructure without imposing
    a transition rule contract. Can be subclassed directly by any model
    that operates on NumPy 2D arrays.

    Parameters (setup)
    ------------------
    backend : RasterBackend
        Backend shared across all models in the same ``Environment``.

    Attributes available in subclasses
    ------------------------------------
    backend : RasterBackend
        The shared array store.
    shape : tuple[int, int]
        Grid shape ``(rows, cols)`` — shortcut for ``self.backend.shape``.
    shift : callable
        Shortcut for ``RasterBackend.shift2d`` (static method).
    dirs : list[tuple[int, int]]
        ``DIRS_MOORE`` — the 8 directions of the Moore neighbourhood.

    Examples
    --------
    >>> class HeatDiffusion(RasterModel):
    ...     def execute(self):
    ...         temp = self.backend.get("temp").copy()
    ...         for dr, dc in self.dirs:
    ...             temp += 0.1 * self.shift(temp, dr, dc)
    ...         self.backend.arrays["temp"] = temp
    """

    # Narrowing the base Model.setup(**kwargs) contract is intentional:
    # Model.__init__ forwards extra kwargs to setup, and each component
    # declares the keywords it accepts.
    def setup(self, backend: RasterBackend) -> None:  # type: ignore[override]
        self.backend = backend
        self.shape   = backend.shape
        self.shift   = RasterBackend.shift2d
        self.dirs    = DIRS_MOORE

dissmodel.geo.raster.sync_model

dissmodel/geo/raster/sync_model.py

SyncRasterModel — RasterModel with automatic _past snapshot semantics.

This module provides the snapshot mechanism equivalent to TerraME's cs:synchronize(), as a reusable base class for any raster model that needs per-step state history (LUCC, fire spread, epidemic models, etc.).

How it works

Before the first execute(), pre_execute() calls synchronize() once to capture the initial state. After each execute(), post_execute() calls synchronize() again to freeze the current state as <name>_past for the next step.

Models can call self.backend.get("f_past") inside execute() to access the state at the beginning of the current step — equivalent to TerraME's cell.past[attr].

Usage

Subclass SyncRasterModel instead of RasterModel and declare self.land_use_types in setup():

class MyRasterModel(SyncRasterModel):
    def setup(self, backend, rate=0.01):
        super().setup(backend)
        self.land_use_types = ["forest", "defor"]
        self.rate = rate

    def execute(self):
        forest_past = self.backend.get("forest_past")
        gain = forest_past * self.rate
        self.backend.arrays["forest"] = forest_past + gain

Relationship to domain libraries

dissluc uses this class as the base for its raster LUCC components, exposing it under the domain-specific alias LUCRasterModel:

# dissluc/raster/core.py
from dissmodel.geo.raster.sync_model import SyncRasterModel as LUCRasterModel

SyncRasterModel

Bases: RasterModel

RasterModel with automatic _past snapshot semantics.

Extends :class:~dissmodel.geo.raster.model.RasterModel with pre_execute() / post_execute() hooks that copy each array listed in self.land_use_types to a <name>_past array in the :class:~dissmodel.geo.raster.backend.RasterBackend before and after every simulation step.

This is the raster analogue of :class:~dissmodel.geo.vector.sync_model.SyncSpatialModel and the Python equivalent of TerraME's cs:synchronize().

Subclass contract

Declare self.land_use_types (list of array names) in setup(). SyncRasterModel will manage all <name>_past arrays automatically. Subclasses must not create or update _past arrays manually.

The backend argument is passed through to :class:~dissmodel.geo.raster.model.RasterModel; any additional keyword arguments are forwarded to the parent class.

Examples:

>>> class ForestRaster(SyncRasterModel):
...     def setup(self, backend, rate=0.01):
...         super().setup(backend)
...         self.land_use_types = ["forest", "defor"]
...         self.rate = rate
...
...     def execute(self):
...         forest_past = self.backend.get("forest_past")
...         gain = forest_past * self.rate
...         self.backend.arrays["forest"] = forest_past + gain

Source code in dissmodel/geo/raster/sync_model.py

class SyncRasterModel(RasterModel):
    """
    ``RasterModel`` with automatic ``_past`` snapshot semantics.

    Extends :class:`~dissmodel.geo.raster.model.RasterModel` with
    ``pre_execute()`` / ``post_execute()`` hooks that copy each array
    listed in ``self.land_use_types`` to a ``<name>_past`` array in the
    :class:`~dissmodel.geo.raster.backend.RasterBackend` before and after
    every simulation step.

    This is the raster analogue of
    :class:`~dissmodel.geo.vector.sync_model.SyncSpatialModel` and the
    Python equivalent of TerraME's ``cs:synchronize()``.

    Subclass contract
    -----------------
    Declare ``self.land_use_types`` (list of array names) in ``setup()``.
    ``SyncRasterModel`` will manage all ``<name>_past`` arrays automatically.
    Subclasses must **not** create or update ``_past`` arrays manually.

    The ``backend`` argument is passed through to
    :class:`~dissmodel.geo.raster.model.RasterModel`; any additional keyword
    arguments are forwarded to the parent class.

    Examples
    --------
    >>> class ForestRaster(SyncRasterModel):
    ...     def setup(self, backend, rate=0.01):
    ...         super().setup(backend)
    ...         self.land_use_types = ["forest", "defor"]
    ...         self.rate = rate
    ...
    ...     def execute(self):
    ...         forest_past = self.backend.get("forest_past")
    ...         gain = forest_past * self.rate
    ...         self.backend.arrays["forest"] = forest_past + gain
    """

    def pre_execute(self) -> None:
        """
        Snapshot arrays before the first step.

        On the first call, freezes the initial state into ``<name>_past``
        arrays so that ``execute()`` can read them. Subsequent calls are
        no-ops — post_execute() handles ongoing snapshots.
        """
        if not getattr(self, "_first_sync_done", False):
            self.synchronize()
            self._first_sync_done = True

    def post_execute(self) -> None:
        """
        Snapshot arrays after each step.

        Freezes the current state into ``<name>_past`` arrays so that
        the next ``execute()`` call reads the state at step start —
        equivalent to TerraME's ``cs:synchronize()``.
        """
        self.synchronize()

    def synchronize(self) -> None:
        """
        Copy each array in ``land_use_types`` to ``<name>_past`` in the backend.

        Equivalent to ``cs:synchronize()`` in TerraME. Called automatically
        via ``pre_execute()`` before the first step and via ``post_execute()``
        after each ``execute()``. Can also be called manually when an explicit
        mid-step snapshot is needed.

        Does nothing if ``land_use_types`` has not been set yet (safe to
        call before ``setup()`` completes).

        If a state variable was loaded from a temporal catalog entry (shape
        ``(time, y, x)``), the first slice is used as the initial state.
        After the first step the model always writes back 2D arrays, so
        subsequent snapshots are 2D unconditionally.
        """
        if not hasattr(self, "land_use_types"):
            return
        for name in self.land_use_types:
            arr = self.backend.get(name)
            if arr.ndim > 2:
                arr = arr[0]
            self.backend.set(name + "_past", arr.copy())

post_execute()

Snapshot arrays after each step.

Freezes the current state into <name>_past arrays so that the next execute() call reads the state at step start — equivalent to TerraME's cs:synchronize().

Source code in dissmodel/geo/raster/sync_model.py

def post_execute(self) -> None:
    """
    Snapshot arrays after each step.

    Freezes the current state into ``<name>_past`` arrays so that
    the next ``execute()`` call reads the state at step start —
    equivalent to TerraME's ``cs:synchronize()``.
    """
    self.synchronize()

pre_execute()

Snapshot arrays before the first step.

On the first call, freezes the initial state into <name>_past arrays so that execute() can read them. Subsequent calls are no-ops — post_execute() handles ongoing snapshots.

Source code in dissmodel/geo/raster/sync_model.py

def pre_execute(self) -> None:
    """
    Snapshot arrays before the first step.

    On the first call, freezes the initial state into ``<name>_past``
    arrays so that ``execute()`` can read them. Subsequent calls are
    no-ops — post_execute() handles ongoing snapshots.
    """
    if not getattr(self, "_first_sync_done", False):
        self.synchronize()
        self._first_sync_done = True

synchronize()

Copy each array in land_use_types to <name>_past in the backend.

Equivalent to cs:synchronize() in TerraME. Called automatically via pre_execute() before the first step and via post_execute() after each execute(). Can also be called manually when an explicit mid-step snapshot is needed.

Does nothing if land_use_types has not been set yet (safe to call before setup() completes).

If a state variable was loaded from a temporal catalog entry (shape (time, y, x)), the first slice is used as the initial state. After the first step the model always writes back 2D arrays, so subsequent snapshots are 2D unconditionally.

Source code in dissmodel/geo/raster/sync_model.py

def synchronize(self) -> None:
    """
    Copy each array in ``land_use_types`` to ``<name>_past`` in the backend.

    Equivalent to ``cs:synchronize()`` in TerraME. Called automatically
    via ``pre_execute()`` before the first step and via ``post_execute()``
    after each ``execute()``. Can also be called manually when an explicit
    mid-step snapshot is needed.

    Does nothing if ``land_use_types`` has not been set yet (safe to
    call before ``setup()`` completes).

    If a state variable was loaded from a temporal catalog entry (shape
    ``(time, y, x)``), the first slice is used as the initial state.
    After the first step the model always writes back 2D arrays, so
    subsequent snapshots are 2D unconditionally.
    """
    if not hasattr(self, "land_use_types"):
        return
    for name in self.land_use_types:
        arr = self.backend.get(name)
        if arr.ndim > 2:
            arr = arr[0]
        self.backend.set(name + "_past", arr.copy())

dissmodel.geo.raster.cellular_automaton

dissmodel/geo/raster_cellular_automaton.py

Base class for cellular automata backed by RasterBackend (NumPy 2D arrays).

Analogous to CellularAutomaton (GeoDataFrame), but for the raster substrate.

Hierarchy

Model
  ├── SpatialModel
  │     └── CellularAutomaton       rule(idx) → value       (vector, pull)
  └── RasterModel
        └── RasterCellularAutomaton rule(arrays) → arrays   (raster, vectorized)

Why a different rule() contract

CellularAutomaton.rule(idx) returns a single value for one cell — it is called once per cell per step (O(n) Python calls). This is correct for the vector substrate where neighborhood lookup is the bottleneck.

For the raster substrate, the bottleneck is the Python loop itself. RasterCellularAutomaton.rule() receives the full snapshot of all arrays and returns a dict of updated arrays — one NumPy call covers the entire grid. This is the natural pattern for NumPy-based CA.

Comparison

# vector CA — rule called n times per step
class GameOfLife(CellularAutomaton):
    def rule(self, idx):
        alive = self.neighbor_values(idx, "state").sum()
        ...
        return new_state

# raster CA — rule called once per step
class GameOfLife(RasterCellularAutomaton):
    def rule(self, arrays):
        state = arrays["state"]
        alive = backend.focal_sum_mask(state == 1)
        ...
        return {"state": new_state}

Usage

from dissmodel.geo.raster_cellular_automaton import RasterCellularAutomaton
from dissmodel.geo.raster.backend import RasterBackend
from dissmodel.core import Environment
import numpy as np

class GameOfLife(RasterCellularAutomaton):
    def rule(self, arrays):
        state     = arrays["state"]
        neighbors = self.backend.focal_sum_mask(state == 1)
        born      = (state == 0) & (neighbors == 3)
        survive   = (state == 1) & np.isin(neighbors, [2, 3])
        return {"state": np.where(born | survive, 1, 0)}

b = RasterBackend(shape=(50, 50))
b.set("state", np.random.randint(0, 2, (50, 50)))

env = Environment(start_time=1, end_time=100)
GameOfLife(backend=b)
env.run()

RasterCellularAutomaton

Bases: RasterModel, ABC

Base class for NumPy-based cellular automata.

Extends :class:~dissmodel.geo.raster.model.RasterModel with a vectorized transition rule — rule() receives all arrays as a snapshot and returns a dict of updated arrays.

See :meth:setup for the keyword arguments accepted by this class.

Examples:

>>> class MyCA(RasterCellularAutomaton):
...     def rule(self, arrays):
...         state = arrays["state"]
...         # ... NumPy operations over full grid ...
...         return {"state": new_state}

Source code in dissmodel/geo/raster/cellular_automaton.py

class RasterCellularAutomaton(RasterModel, ABC):
    """
    Base class for NumPy-based cellular automata.

    Extends :class:`~dissmodel.geo.raster.model.RasterModel` with a
    vectorized transition rule — ``rule()`` receives all arrays as a
    snapshot and returns a dict of updated arrays.

    See :meth:`setup` for the keyword arguments accepted by this class.

    Examples
    --------
    >>> class MyCA(RasterCellularAutomaton):
    ...     def rule(self, arrays):
    ...         state = arrays["state"]
    ...         # ... NumPy operations over full grid ...
    ...         return {"state": new_state}
    """

    # Narrowing the base Model.setup(**kwargs) contract is intentional:
    # Model.__init__ forwards extra kwargs to setup, and each component
    # declares the keywords it accepts.
    def setup(  # type: ignore[override]
        self,
        backend:    RasterBackend,
        state_attr: str = "state"
    ) -> None:
        """
        Configure the cellular automaton.

        Called automatically by ``Model.__init__`` with any keyword
        arguments passed to the constructor.

        Parameters
        ----------
        backend : RasterBackend
            Shared backend with the simulation arrays.
        state_attr : str, optional
            Primary state array name, by default ``"state"``.
            Used only for introspection/logging — rule() can update any array.
        """
        super().setup(backend)
        self.state_attr = state_attr

    @abstractmethod
    def rule(self, arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]:
        """
        Vectorized transition rule applied to the full grid.

        Receives a snapshot of all arrays (equivalent to celula.past[] in
        TerraME) and returns a dict with the arrays to update.

        Only the arrays present in the returned dict are written back —
        arrays not returned are left unchanged.

        Parameters
        ----------
        arrays : dict[str, np.ndarray]
            Snapshot of backend arrays at the start of the step.
            Modifying these arrays does NOT affect the backend — they are
            copies (equivalent to .past semantics).

        Returns
        -------
        dict[str, np.ndarray]
            Dict mapping array name → new array. Partial updates allowed.

        Examples
        --------
        >>> def rule(self, arrays):
        ...     state     = arrays["state"]           # read from snapshot
        ...     neighbors = self.backend.focal_sum_mask(state == 1)
        ...     new_state = np.where(neighbors > 3, 0, state)
        ...     return {"state": new_state}           # write back
        """
        raise NotImplementedError("Subclasses must implement rule().")

    def execute(self) -> None:
        """
        Execute one simulation step by calling rule() once over the full grid.

        Takes a snapshot of all arrays (past state), passes it to rule(),
        and writes the returned arrays back to the backend.
        """
        past    = self.backend.snapshot()   # equivale a celula.past[]
        updates = self.rule(past)
        for name, arr in updates.items():
            self.backend.arrays[name] = arr

execute()

Execute one simulation step by calling rule() once over the full grid.

Takes a snapshot of all arrays (past state), passes it to rule(), and writes the returned arrays back to the backend.

Source code in dissmodel/geo/raster/cellular_automaton.py

def execute(self) -> None:
    """
    Execute one simulation step by calling rule() once over the full grid.

    Takes a snapshot of all arrays (past state), passes it to rule(),
    and writes the returned arrays back to the backend.
    """
    past    = self.backend.snapshot()   # equivale a celula.past[]
    updates = self.rule(past)
    for name, arr in updates.items():
        self.backend.arrays[name] = arr

rule(arrays) abstractmethod

Vectorized transition rule applied to the full grid.

Receives a snapshot of all arrays (equivalent to celula.past[] in TerraME) and returns a dict with the arrays to update.

Only the arrays present in the returned dict are written back — arrays not returned are left unchanged.

Parameters:

Name	Type	Description	Default
arrays	dict[str, ndarray]	Snapshot of backend arrays at the start of the step. Modifying these arrays does NOT affect the backend — they are copies (equivalent to .past semantics).	required

Returns:

Type	Description
dict[str, ndarray]	Dict mapping array name → new array. Partial updates allowed.

Examples:

>>> def rule(self, arrays):
...     state     = arrays["state"]           # read from snapshot
...     neighbors = self.backend.focal_sum_mask(state == 1)
...     new_state = np.where(neighbors > 3, 0, state)
...     return {"state": new_state}           # write back

Source code in dissmodel/geo/raster/cellular_automaton.py

@abstractmethod
def rule(self, arrays: dict[str, np.ndarray]) -> dict[str, np.ndarray]:
    """
    Vectorized transition rule applied to the full grid.

    Receives a snapshot of all arrays (equivalent to celula.past[] in
    TerraME) and returns a dict with the arrays to update.

    Only the arrays present in the returned dict are written back —
    arrays not returned are left unchanged.

    Parameters
    ----------
    arrays : dict[str, np.ndarray]
        Snapshot of backend arrays at the start of the step.
        Modifying these arrays does NOT affect the backend — they are
        copies (equivalent to .past semantics).

    Returns
    -------
    dict[str, np.ndarray]
        Dict mapping array name → new array. Partial updates allowed.

    Examples
    --------
    >>> def rule(self, arrays):
    ...     state     = arrays["state"]           # read from snapshot
    ...     neighbors = self.backend.focal_sum_mask(state == 1)
    ...     new_state = np.where(neighbors > 3, 0, state)
    ...     return {"state": new_state}           # write back
    """
    raise NotImplementedError("Subclasses must implement rule().")

setup(backend, state_attr='state')

Configure the cellular automaton.

Called automatically by Model.__init__ with any keyword arguments passed to the constructor.

Parameters:

Name	Type	Description	Default
backend	RasterBackend	Shared backend with the simulation arrays.	required
state_attr	str	Primary state array name, by default "state". Used only for introspection/logging — rule() can update any array.	'state'

Source code in dissmodel/geo/raster/cellular_automaton.py

def setup(  # type: ignore[override]
    self,
    backend:    RasterBackend,
    state_attr: str = "state"
) -> None:
    """
    Configure the cellular automaton.

    Called automatically by ``Model.__init__`` with any keyword
    arguments passed to the constructor.

    Parameters
    ----------
    backend : RasterBackend
        Shared backend with the simulation arrays.
    state_attr : str, optional
        Primary state array name, by default ``"state"``.
        Used only for introspection/logging — rule() can update any array.
    """
    super().setup(backend)
    self.state_attr = state_attr

dissmodel.geo.raster.raster_grid

dissmodel/geo/raster_grid.py

Utilitário para criar RasterBackend sintético.

Análogo a vector_grid() (GeoDataFrame), mas para o substrato NumPy.

Uso

from dissmodel.geo.raster_grid import raster_grid
import numpy as np

# grade vazia com arrays zerados
b = raster_grid(rows=50, cols=50, attrs={"state": 0})

# grade com array inicial customizado
b = raster_grid(
    rows=50, cols=50,
    attrs={"state": np.random.randint(0, 2, (50, 50))}
)

raster_grid(rows, cols, attrs=None, dtype=None)

Create a RasterBackend with optional pre-filled arrays.

Analogous to :func:~dissmodel.geo.vector_grid for the raster substrate. Useful for tests, examples, and synthetic benchmarks.

Parameters:

Name	Type	Description	Default
rows	int	Number of rows in the grid.	required
cols	int	Number of columns in the grid.	required
attrs	dict	Mapping of array name → initial value. - scalar (int or float): fills the entire grid with that value. - np.ndarray of shape (rows, cols): used directly (a copy is stored). If not provided, an empty backend is returned.	None
dtype	numpy dtype	Default dtype for scalar-initialized arrays. If None, inferred from the scalar type (int → np.int32, float → np.float64).	None

Returns:

Type	Description
RasterBackend	Backend with shape (rows, cols) and the requested arrays.

Examples:

>>> b = raster_grid(10, 10, attrs={"state": 0})
>>> b.shape
(10, 10)
>>> b.get("state").shape
(10, 10)

>>> import numpy as np
>>> state = np.random.randint(0, 2, (10, 10))
>>> b = raster_grid(10, 10, attrs={"state": state})

Source code in dissmodel/geo/raster/raster_grid.py

def raster_grid(
    rows:  int,
    cols:  int,
    attrs: dict[str, AttrValue] | None = None,
    dtype: np.dtype | None             = None,
) -> RasterBackend:
    """
    Create a RasterBackend with optional pre-filled arrays.

    Analogous to :func:`~dissmodel.geo.vector_grid` for the raster
    substrate. Useful for tests, examples, and synthetic benchmarks.

    Parameters
    ----------
    rows : int
        Number of rows in the grid.
    cols : int
        Number of columns in the grid.
    attrs : dict, optional
        Mapping of array name → initial value.
        - scalar (int or float): fills the entire grid with that value.
        - np.ndarray of shape (rows, cols): used directly (a copy is stored).
        If not provided, an empty backend is returned.
    dtype : numpy dtype, optional
        Default dtype for scalar-initialized arrays. If None, inferred
        from the scalar type (int → np.int32, float → np.float64).

    Returns
    -------
    RasterBackend
        Backend with shape (rows, cols) and the requested arrays.

    Examples
    --------
    >>> b = raster_grid(10, 10, attrs={"state": 0})
    >>> b.shape
    (10, 10)
    >>> b.get("state").shape
    (10, 10)

    >>> import numpy as np
    >>> state = np.random.randint(0, 2, (10, 10))
    >>> b = raster_grid(10, 10, attrs={"state": state})
    """
    b = RasterBackend(shape=(rows, cols))

    for name, value in (attrs or {}).items():
        if isinstance(value, np.ndarray):
            if value.shape != (rows, cols):
                raise ValueError(
                    f"Array '{name}' has shape {value.shape}, "
                    f"expected ({rows}, {cols})."
                )
            b.set(name, value)
        else:
            # scalar → infer dtype
            if dtype is not None:
                arr_dtype = dtype
            elif isinstance(value, float):
                arr_dtype = np.dtype(np.float64)
            else:
                arr_dtype = np.dtype(np.int32)
            b.set(name, np.full((rows, cols), value, dtype=arr_dtype))

    return b

dissmodel.geo.raster.band_spec

BandSpec dataclass

Specification of a raster band in a GeoTIFF.

Attributes:

Name	Type	Description
name	str	Name used inside RasterBackend (e.g. 'uso', 'alt', 'soil').
dtype	str	NumPy dtype used to store the band.
nodata	float | int	Value representing missing data.

Source code in dissmodel/geo/raster/band_spec.py

@dataclass
class BandSpec:
    """
    Specification of a raster band in a GeoTIFF.

    Attributes
    ----------
    name : str
        Name used inside RasterBackend (e.g. 'uso', 'alt', 'soil').
    dtype : str
        NumPy dtype used to store the band.
    nodata : float | int
        Value representing missing data.
    """

    name: str
    dtype: str
    nodata: float | int