Grid operations

The grid operations module provides functions for processing and transforming gridded raster data, including resolution decimation and flow direction conversions.

Overview

Grid operations are essential for:

Multi-resolution modeling: Coarsen high-resolution DEMs, flow direction, and accumulation grids
Flow direction processing: Handle different flow direction notation systems (Grass r.watershed vs ArcGIS)
Preprocessing: Prepare gridded data for hydrological modeling

All functions handle NaN values.

Flow direction notations

MOBIDICpy supports two flow direction notation systems. All diagrams below are shown in geographic orientation (north up).

Grass notation (1-8)

Standard GRASS r.watershed convention: codes 1-8 placed counter-clockwise starting from NE. To use this notation, the option raster_settings.flow_dir_type in the YAML configuration file must be set to “Grass”.

Grass flow direction notation

Arc notation (powers of 2)

Standard ESRI ArcGIS convention: powers of 2 placed clockwise starting from E. To use this notation, the option raster_settings.flow_dir_type in the YAML configuration file must be set to “Arc”.

Arc flow direction notation

MOBIDIC notation (1-8)

The flow directions, either in Grass or Arc notation, are then internally converted to MOBIDIC’s D8 encoding as follows:

MOBIDIC flow direction notation

Technical details

Degradation algorithm

Divides the input grid into blocks of size factor × factor
For regular rasters: computes mean of valid cells in each block
For flow direction: finds the cell with maximum flow accumulation in each block, determines the dominant flow direction
Applies min_valid_fraction threshold to avoid blocks with too few valid cells

Flow direction degradation

The algorithm preserves drainage patterns by:

Finding the fine cell with maximum flow accumulation in each coarse block
Determining which coarse neighbor the drainage flows to
Assigning the appropriate flow direction code
Normalizing flow accumulation by factor × factor to maintain consistent scaling

Notes

All functions return new arrays and transforms without modifying inputs
NaN values are properly propagated and excluded from calculations
Flow direction values must be in the valid range for the specified notation
Invalid flow direction values (e.g., not in Grass 1-8 or Arc powers-of-2) are converted to NaN

Functions

Resolution decimation

`mobidic.preprocessing.grid_operations.decimate_raster(data, factor, min_valid_fraction=0.125)`

Coarsen raster resolution by aggregating cells.

Aggregates a fine-resolution grid into a coarser grid by averaging values within each block of size factor x factor. Blocks with too few valid cells (less than min_valid_fraction) are marked as NaN.

Parameters:

Name	Type	Description	Default
`data`	`ndarray`	2D numpy array with raster values (NaN for nodata).	required
`factor`	`int`	Decimation factor (e.g., 2 means 2x2 blocks -> 1 cell).	required
`min_valid_fraction`	`float`	Minimum fraction of valid cells required in each block to compute the mean. Default is 0.125 (1/8 of cells).	`0.125`

Returns:

Type	Description
`ndarray`	Decimated 2D numpy array with shape (floor(nr/factor), floor(nc/factor)).

Examples:

>>> import numpy as np
>>> data = np.random.rand(100, 100)
>>> decimated = decimate_raster(data, factor=2)
>>> decimated.shape
(50, 50)

Source code in mobidic/preprocessing/grid_operations.py

def decimate_raster(
    data: np.ndarray,
    factor: int,
    min_valid_fraction: float = 0.125,
) -> np.ndarray:
    """Coarsen raster resolution by aggregating cells.

    Aggregates a fine-resolution grid into a coarser grid by averaging values
    within each block of size factor x factor. Blocks with too few valid cells
    (less than min_valid_fraction) are marked as NaN.

    Args:
        data: 2D numpy array with raster values (NaN for nodata).
        factor: Decimation factor (e.g., 2 means 2x2 blocks -> 1 cell).
        min_valid_fraction: Minimum fraction of valid cells required in each block
            to compute the mean. Default is 0.125 (1/8 of cells).

    Returns:
        Decimated 2D numpy array with shape (floor(nr/factor), floor(nc/factor)).

    Examples:
        >>> import numpy as np
        >>> data = np.random.rand(100, 100)
        >>> decimated = decimate_raster(data, factor=2)
        >>> decimated.shape
        (50, 50)
    """
    if factor < 1:
        logger.error(f"Decimation factor must be >= 1, got {factor}")
        raise ValueError(f"Decimation factor must be >= 1, got {factor}")

    if factor == 1:
        logger.debug("Decimation factor is 1, returning original data")
        return data.copy()

    nr, nc = data.shape
    nrv = nr // factor
    ncv = nc // factor

    logger.info(f"Decimating raster from {data.shape} to ({nrv}, {ncv}) with factor {factor}")

    # Initialize output array
    decimated = np.full((nrv, ncv), np.nan, dtype=float)

    # Minimum number of valid cells required
    min_valid_cells = int(factor * factor * min_valid_fraction)

    # Process each block
    for i in range(nrv):
        iv_start = i * factor
        iv_end = (i + 1) * factor

        for j in range(ncv):
            jv_start = j * factor
            jv_end = (j + 1) * factor

            # Extract block
            block = data[iv_start:iv_end, jv_start:jv_end]

            # Count valid (non-NaN) cells
            n_valid = np.sum(np.isfinite(block))

            # Compute mean if enough valid cells
            if n_valid >= min_valid_cells:
                decimated[i, j] = np.nanmean(block)

    logger.success(f"Raster decimation complete: {np.sum(np.isfinite(decimated))} valid cells")

    return decimated

`mobidic.preprocessing.grid_operations.decimate_flow_direction(flow_dir, flow_acc, factor, min_valid_fraction=0.125)`

Coarsen flow direction and flow accumulation grids.

Aggregates fine-resolution flow direction and flow accumulation grids into coarser grids. For each coarse cell, identifies the fine cell with maximum flow accumulation and determines the coarse flow direction based on where that cell drains to.

Parameters:

Name	Type	Description	Default
`flow_dir`	`ndarray`	2D numpy array with flow directions (1-8 notation, NaN for nodata).	required
`flow_acc`	`ndarray`	2D numpy array with flow accumulation values (NaN for nodata).	required
`factor`	`int`	Decimation factor (e.g., 2 means 2x2 blocks -> 1 cell).	required
`min_valid_fraction`	`float`	Minimum fraction of valid cells required in each block. Default is 0.125 (1/8 of cells).	`0.125`

Returns:

Type	Description
`ndarray`	Tuple of (decimated_flow_dir, decimated_flow_acc) as 2D numpy arrays.
`ndarray`	Flow accumulation is normalized by factor **2 to account for cell size change.

Notes

This function operates on flow_dir BEFORE conversion to MOBIDIC notation, so input is expected in standard GRASS r.watershed 1-8 notation: 1=NE, 2=N, 3=NW, 4=W, 5=SW, 6=S, 7=SE, 8=E

Examples:

>>> flow_dir = np.array([[2, 2], [2, 2]])  # All cells flow north (Grass code 2)
>>> flow_acc = np.array([[1, 1], [3, 4]])  # Different accumulation
>>> dec_dir, dec_acc = decimate_flow_direction(flow_dir, flow_acc, factor=2)

Source code in mobidic/preprocessing/grid_operations.py

def decimate_flow_direction(
    flow_dir: np.ndarray,
    flow_acc: np.ndarray,
    factor: int,
    min_valid_fraction: float = 0.125,
) -> tuple[np.ndarray, np.ndarray]:
    """Coarsen flow direction and flow accumulation grids.

    Aggregates fine-resolution flow direction and flow accumulation grids into
    coarser grids. For each coarse cell, identifies the fine cell with maximum
    flow accumulation and determines the coarse flow direction based on where
    that cell drains to.

    Args:
        flow_dir: 2D numpy array with flow directions (1-8 notation, NaN for nodata).
        flow_acc: 2D numpy array with flow accumulation values (NaN for nodata).
        factor: Decimation factor (e.g., 2 means 2x2 blocks -> 1 cell).
        min_valid_fraction: Minimum fraction of valid cells required in each block.
            Default is 0.125 (1/8 of cells).

    Returns:
        Tuple of (decimated_flow_dir, decimated_flow_acc) as 2D numpy arrays.
        Flow accumulation is normalized by factor **2 to account for cell size change.

    Notes:
        This function operates on flow_dir BEFORE conversion to MOBIDIC notation,
        so input is expected in standard GRASS r.watershed 1-8 notation:
            1=NE, 2=N, 3=NW, 4=W, 5=SW, 6=S, 7=SE, 8=E

    Examples:
        >>> flow_dir = np.array([[2, 2], [2, 2]])  # All cells flow north (Grass code 2)
        >>> flow_acc = np.array([[1, 1], [3, 4]])  # Different accumulation
        >>> dec_dir, dec_acc = decimate_flow_direction(flow_dir, flow_acc, factor=2)
    """
    if factor < 1:
        logger.error(f"Decimation factor must be >= 1, got {factor}")
        raise ValueError(f"Decimation factor must be >= 1, got {factor}")

    if factor == 1:
        logger.debug("Decimation factor is 1, returning original data")
        return flow_dir.copy(), flow_acc.copy()

    if flow_dir.shape != flow_acc.shape:
        logger.error("Flow direction and accumulation grids must have the same shape")
        raise ValueError("Flow direction and accumulation grids must have the same shape")

    nr, nc = flow_dir.shape
    nrv = nr // factor
    ncv = nc // factor

    logger.info(f"Decimating flow direction/accumulation from {flow_dir.shape} to ({nrv}, {ncv})")

    # Initialize output arrays
    flow_dir_decimated = np.full((nrv, ncv), np.nan, dtype=float)
    flow_acc_decimated = np.full((nrv, ncv), np.nan, dtype=float)

    # Direction offsets indexed by code-1, in standard GRASS r.watershed notation
    # (1=NE, 2=N, 3=NW, 4=W, 5=SW, 6=S, 7=SE, 8=E). Offsets are matrix-coord
    # (row, col) where row+1 = down (image orientation, before any flipud).
    di = np.array([-1, -1, -1, 0, 1, 1, 1, 0])  # row offset
    dj = np.array([1, 0, -1, -1, -1, 0, 1, 1])  # column offset

    # Minimum number of valid cells required
    min_valid_cells = int(factor * factor * min_valid_fraction)

    # Process each coarse block
    for i in range(nrv):
        iv_start = i * factor
        iv_end = (i + 1) * factor

        for j in range(ncv):
            jv_start = j * factor
            jv_end = (j + 1) * factor

            # Extract blocks
            block_dir = flow_dir[iv_start:iv_end, jv_start:jv_end]
            block_acc = flow_acc[iv_start:iv_end, jv_start:jv_end]

            # Count valid cells
            n_valid = np.sum(np.isfinite(block_dir))

            # Skip if not enough valid cells
            if n_valid < min_valid_cells:
                continue

            # Set NaN in accumulation where direction is NaN
            block_acc_masked = block_acc.copy()
            block_acc_masked[np.isnan(block_dir)] = np.nan

            # Find cell with maximum flow accumulation
            if not np.any(np.isfinite(block_acc_masked)):
                continue

            max_idx = np.nanargmax(block_acc_masked)
            max_i, max_j = np.unravel_index(max_idx, block_acc_masked.shape)

            # Store normalized flow accumulation
            flow_acc_decimated[i, j] = np.nanmax(block_acc_masked) / (factor * factor)

            # Get flow direction of the cell with max accumulation
            direction = int(block_dir[max_i, max_j])

            if not (1 <= direction <= 8):
                # Invalid direction, try to find a valid neighbor
                flow_dir_decimated[i, j] = -999
                continue

            # Calculate which coarse cell this fine cell drains to
            target_i = max_i + di[direction - 1]
            target_j = max_j + dj[direction - 1]

            # Determine coarse cell index offset
            coarse_di = target_i // factor
            coarse_dj = target_j // factor

            # Calculate target coarse cell coordinates
            target_coarse_i = i + coarse_di
            target_coarse_j = j + coarse_dj

            # Check if target is within bounds
            if not (0 <= target_coarse_i < nrv and 0 <= target_coarse_j < ncv):
                flow_dir_decimated[i, j] = -999
                continue

            # Map coarse cell offset to direction (1-8) in standard GRASS r.watershed
            # convention (image orientation, row+1 = down).
            offset_to_dir = {
                (-1, 1): 1,  # NE
                (-1, 0): 2,  # N
                (-1, -1): 3,  # NW
                (0, -1): 4,  # W
                (1, -1): 5,  # SW
                (1, 0): 6,  # S
                (1, 1): 7,  # SE
                (0, 1): 8,  # E
            }

            flow_dir_decimated[i, j] = offset_to_dir.get((coarse_di, coarse_dj), -999)

    # Fix invalid directions by finding valid neighbors
    invalid_mask = flow_dir_decimated == -999
    if np.any(invalid_mask):
        logger.debug(f"Fixing {np.sum(invalid_mask)} invalid flow directions")

        invalid_indices = np.argwhere(invalid_mask)
        for idx in invalid_indices:
            i, j = idx

            # Try each direction to find a valid neighbor
            for d in range(8):
                ii = i + di[d]
                jj = j + dj[d]

                if 0 <= ii < nrv and 0 <= jj < ncv:
                    if flow_dir_decimated[ii, jj] > 0:
                        flow_dir_decimated[i, j] = d + 1
                        break

    logger.success(
        f"Flow direction decimation complete: "
        f"{np.sum(np.isfinite(flow_dir_decimated))} valid cells, "
        f"{np.sum(flow_dir_decimated == -999)} invalid directions"
    )

    return flow_dir_decimated, flow_acc_decimated

Flow direction conversion

`mobidic.preprocessing.grid_operations.convert_to_mobidic_notation(flow_dir, from_notation='Grass')`

Convert flow direction from Grass or Arc notation to MOBIDIC notation.

Translates standard GRASS r.watershed (1-8) or ESRI ArcGIS (powers of 2) flow direction codes to MOBIDIC’s internal D8 encoding. The renumbering is direction-preserving (each input code maps to the MOBIDIC code with the same physical compass direction). The mapping matches MATLAB’s buildgis_mysql_include.m: GRASS: AI=[1 2 3 4 5 6 7 8]; MD=[5 6 7 8 1 2 3 4]; ArcGIS: AI=[1 2 4 8 16 32 64 128]; MD=[4 3 2 1 8 7 6 5];

Parameters:

Name	Type	Description	Default
`flow_dir`	`ndarray`	2D numpy array with flow directions (NaN for nodata).	required
`from_notation`	`Literal['Grass', 'Arc']`	Source notation (‘Grass’ for 1-8 or ‘Arc’ for power-of-2). Default is ‘Grass’.	`'Grass'`

Returns:

Type	Description
`ndarray`	Flow direction array converted to MOBIDIC notation (1-8).

Notes

Standard GRASS r.watershed convention (1-8, CCW from NE; north up):

3 2 1            1 = NE       5 = SW
4 . 8            2 = N        6 = S
5 6 7            3 = NW       7 = SE
                 4 = W        8 = E

Standard ESRI ArcGIS convention (powers of 2, CW from E; north up):

 32  64 128       1   = E       16  = W
 16   .   1       2   = SE      32  = NW
  8   4   2       4   = S       64  = N
                  8   = SW      128 = NE

MOBIDIC convention (1-8, CCW from SW; north up):

7 6 5            1 = SW       5 = NE
8 . 4            2 = S        6 = N
1 2 3            3 = SE       7 = NW
                 4 = E        8 = W

Direction-preserving mapping: GRASS -> MOBIDIC: 1->5, 2->6, 3->7, 4->8, 5->1, 6->2, 7->3, 8->4 ArcGIS -> MOBIDIC: 1->4, 2->3, 4->2, 8->1, 16->8, 32->7, 64->6, 128->5

Array orientation: rasters are loaded by grid_to_matrix() with np.flipud so that increasing row index corresponds to moving NORTH. With that orientation, MOBIDIC code k drains to the neighbor at offset (di, dj) where di = [-1,-1,-1,0,1,1,1,0][k-1] and dj = [-1,0,1,1,1,0,-1,-1][k-1]. For example, code 5 (NE) -> (+1, +1).

Examples:

>>> flow_dir_grass = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 1]])
>>> flow_dir_mobidic = convert_to_mobidic_notation(flow_dir_grass, 'Grass')
>>> flow_dir_mobidic
array([[5, 6, 7],
       [8, 1, 2],
       [3, 4, 5]])

Source code in mobidic/preprocessing/grid_operations.py

def convert_to_mobidic_notation(
    flow_dir: np.ndarray,
    from_notation: Literal["Grass", "Arc"] = "Grass",
) -> np.ndarray:
    """Convert flow direction from Grass or Arc notation to MOBIDIC notation.

    Translates standard GRASS r.watershed (1-8) or ESRI ArcGIS (powers of 2) flow
    direction codes to MOBIDIC's internal D8 encoding. The renumbering is
    direction-preserving (each input code maps to the MOBIDIC code with the same
    physical compass direction). The mapping matches MATLAB's
    buildgis_mysql_include.m:
        GRASS:  AI=[1 2 3 4 5 6 7 8];        MD=[5 6 7 8 1 2 3 4];
        ArcGIS: AI=[1 2 4 8 16 32 64 128];   MD=[4 3 2 1 8 7 6 5];

    Args:
        flow_dir: 2D numpy array with flow directions (NaN for nodata).
        from_notation: Source notation ('Grass' for 1-8 or 'Arc' for power-of-2). Default is 'Grass'.

    Returns:
        Flow direction array converted to MOBIDIC notation (1-8).

    Notes:
        Standard GRASS r.watershed convention (1-8, CCW from NE; north up):

            3 2 1            1 = NE       5 = SW
            4 . 8            2 = N        6 = S
            5 6 7            3 = NW       7 = SE
                             4 = W        8 = E

        Standard ESRI ArcGIS convention (powers of 2, CW from E; north up):

             32  64 128       1   = E       16  = W
             16   .   1       2   = SE      32  = NW
              8   4   2       4   = S       64  = N
                              8   = SW      128 = NE

        MOBIDIC convention (1-8, CCW from SW; north up):

            7 6 5            1 = SW       5 = NE
            8 . 4            2 = S        6 = N
            1 2 3            3 = SE       7 = NW
                             4 = E        8 = W

        Direction-preserving mapping:
            GRASS  -> MOBIDIC:  1->5, 2->6, 3->7, 4->8, 5->1, 6->2, 7->3, 8->4
            ArcGIS -> MOBIDIC:  1->4, 2->3, 4->2, 8->1, 16->8, 32->7, 64->6, 128->5

        Array orientation: rasters are loaded by `grid_to_matrix()` with `np.flipud`
        so that increasing row index corresponds to moving NORTH. With that
        orientation, MOBIDIC code k drains to the neighbor at offset
        (di, dj) where di = `[-1,-1,-1,0,1,1,1,0][k-1]` and
        dj = `[-1,0,1,1,1,0,-1,-1][k-1]`. For example, code 5 (NE) -> (+1, +1).

    Examples:
        >>> flow_dir_grass = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 1]])
        >>> flow_dir_mobidic = convert_to_mobidic_notation(flow_dir_grass, 'Grass')
        >>> flow_dir_mobidic
        array([[5, 6, 7],
               [8, 1, 2],
               [3, 4, 5]])
    """

    # GRASS to MOBIDIC transformation (from buildgis_mysql_include.m)
    grass_to_mobidic = {1: 5, 2: 6, 3: 7, 4: 8, 5: 1, 6: 2, 7: 3, 8: 4}

    # Arc to GRASS mapping
    arc_to_grass = {1: 8, 2: 7, 4: 6, 8: 5, 16: 4, 32: 3, 64: 2, 128: 1}

    # Create output array
    converted = flow_dir.copy()

    if from_notation == "Arc":
        # Arc -> Grass -> MOBIDIC (two-step conversion)
        logger.debug("Converting Arc -> Grass -> MOBIDIC")
        for arc_val, grass_val in arc_to_grass.items():
            mask = flow_dir == arc_val
            converted[mask] = grass_to_mobidic[grass_val]
    elif from_notation == "Grass":
        # Grass -> MOBIDIC (direct conversion)
        logger.debug("Converting Grass -> MOBIDIC")
        for grass_val, mobidic_val in grass_to_mobidic.items():
            mask = flow_dir == grass_val
            converted[mask] = mobidic_val
    else:
        raise ValueError(f"Invalid from_notation: {from_notation}. Must be 'Grass' or 'Arc'")

    valid_cells = np.sum(np.isfinite(converted))
    logger.info(f"Flow direction conversion to MOBIDIC complete: {valid_cells} cells converted")

    return converted

Examples

Coarsening a raster

from mobidic import grid_to_matrix, decimate_raster
import numpy as np

# Read high-resolution DTM (e.g., 10m)
dtm = grid_to_matrix("dtm_10m.tif")

# Degrade to 50m resolution (factor = 5)
dtm_decimated = decimate_raster(
    data=dtm['data'],
    factor=5,
    min_valid_fraction=0.125  # Require at least 1/8 valid cells
)

print(f"Original shape: {dtm['data'].shape}")
print(f"Decimated shape: {dtm_decimated.shape}")
print(f"Original cellsize: {dtm['cellsize']} m")
print(f"New cellsize: {dtm['cellsize'] * 5} m")

Coarsening flow direction

from mobidic import grid_to_matrix, decimate_flow_direction

# Read flow direction and accumulation grids
flow_dir_data = grid_to_matrix("flow_direction.tif")
flow_acc_data = grid_to_matrix("flow_accumulation.tif")

# Degrade both grids together
flow_dir_coarse, flow_acc_coarse = decimate_flow_direction(
    flow_dir=flow_dir_data['data'],
    flow_acc=flow_acc_data['data'],
    factor=5,
    min_valid_fraction=0.5
)

print(f"Original shape: {flow_dir_data['data'].shape}")
print(f"Decimated shape: {flow_dir_coarse.shape}")

Converting flow direction notation

from mobidic import grid_to_matrix, convert_to_mobidic_notation

# Read flow direction in Grass notation (1-8)
flow_dir_grass = grid_to_matrix("flow_dir_grass.tif")

# Convert to MOBIDIC notation (used internally by the model)
flow_dir_mobidic = convert_to_mobidic_notation(
    flow_dir=flow_dir_grass['data'],
    from_notation="Grass"
)

# Or convert from Arc notation directly
flow_dir_arc = grid_to_matrix("flow_dir_arc.tif")
flow_dir_mobidic = convert_to_mobidic_notation(
    flow_dir=flow_dir_arc['data'],
    from_notation="Arc"
)