GIS data I/O

The GIS I/O module provides functions for reading geospatial data in various formats, including raster (GeoTIFF) and vector (Shapefile) files.

Overview

MOBIDICpy uses industry-standard libraries for geospatial data handling:

• Rasterio for raster data (GeoTIFF)
• GeoPandas for vector data (Shapefiles, GeoJSON, etc.)

All functions provide comprehensive logging and error handling, converting nodata values to NaN for consistent numerical processing.

Functions

Raster I/O

mobidic.preprocessing.gis_reader.grid_to_matrix(gridname)

Read raster grid file and return data array with corner coordinates. This function replicates the behavior of MATLAB’s function for GeoTIFF files.

Reads raster files in GeoTIFF (.tif, .tiff) format. Returns the data array and corner coordinates adjusted to cell centers.

Parameters:

Name	Type	Description	Default
gridname	str | Path	Path to the grid file (.tif, .tiff).	required

Returns:

Type	Description
dict[str, Any]	Dictionary containing: - ‘data’: 2D numpy array with raster values (NaN for nodata) - ‘xllcorner’: X coordinate of lower-left corner (cell center) - ‘yllcorner’: Y coordinate of lower-left corner (cell center) - ‘cellsize’: Cell size in map units

Raises:

Type	Description
FileNotFoundError	If the grid file does not exist.
ValueError	If the file format is not supported.
RuntimeError	If there are errors reading the file.

Notes

• Corner coordinates are adjusted to cell centers (0.5 * cellsize offset)
• Data is flipped vertically to match MATLAB convention
• Very small values (< -1e32) are converted to NaN

Examples:

>>> result = grid_to_matrix('elevation.tif')
>>> data = result['data']
>>> print(f"Shape: {data.shape}, Resolution: {result['cellsize']}m")

Source code in mobidic/preprocessing/gis_reader.py

def grid_to_matrix(gridname: str | Path) -> dict[str, Any]:
    """Read raster grid file and return data array with corner coordinates. This function
    replicates the behavior of MATLAB's function for GeoTIFF files.

    Reads raster files in GeoTIFF (.tif, .tiff)
    format. Returns the data array and corner coordinates adjusted to cell centers.

    Args:
        gridname: Path to the grid file (.tif, .tiff).

    Returns:
        Dictionary containing:
            - 'data': 2D numpy array with raster values (NaN for nodata)
            - 'xllcorner': X coordinate of lower-left corner (cell center)
            - 'yllcorner': Y coordinate of lower-left corner (cell center)
            - 'cellsize': Cell size in map units

    Raises:
        FileNotFoundError: If the grid file does not exist.
        ValueError: If the file format is not supported.
        RuntimeError: If there are errors reading the file.

    Notes:
        - Corner coordinates are adjusted to cell centers (0.5 * cellsize offset)
        - Data is flipped vertically to match MATLAB convention
        - Very small values (< -1e32) are converted to NaN

    Examples:
        >>> result = grid_to_matrix('elevation.tif')
        >>> data = result['data']
        >>> print(f"Shape: {data.shape}, Resolution: {result['cellsize']}m")
    """
    gridname = Path(gridname)
    suffix_lower = gridname.suffix.lower()

    logger.info(f"Reading grid file: {gridname}")

    # GeoTIFF format
    if suffix_lower in [".tif", ".tiff"]:
        if not gridname.exists():
            logger.error(f"File not found: {gridname}")
            raise FileNotFoundError(f"File not found: {gridname}")

        try:
            with rasterio.open(gridname) as src:
                # Read data and convert to float
                matgr = src.read(1).astype(float)

                # Get nodata value and convert to NaN
                nodata = src.nodata
                if nodata is not None:
                    matgr[matgr == nodata] = np.nan
                else:
                    # Without a declared nodata value, common sentinels are read as valid data
                    for sentinel in _NODATA_TYPICAL:
                        n_sentinel = int(np.sum(matgr == sentinel))
                        if n_sentinel > 0:
                            logger.warning(
                                f"{gridname} declares no nodata value but has {n_sentinel} cells equal to "
                                f"{sentinel:g}: they are read as valid data. If they are nodata, set the "
                                "nodata value of the file."
                            )

                # Get spatial information
                transform = src.transform
                cellsize = transform[0]  # Assuming square pixels
                xllcorner = src.bounds.left
                yllcorner = src.bounds.bottom

                # Get CRS
                crs = src.crs

        except Exception as e:
            logger.error(f"Error reading GeoTIFF file {gridname}: {e}")
            raise RuntimeError(f"Error reading GeoTIFF file: {e}") from e

        # Flip vertically (to match MATLAB convention)
        matgr = np.flipud(matgr)

        # Filter nodata values (e.g., -3.4e38)
        matgr[matgr < -1e32] = np.nan

        # Adjust to cell center
        xllcorner += 0.5 * cellsize
        yllcorner += 0.5 * cellsize

        logger.info(f"GeoTIFF read: shape={matgr.shape}, cellsize={cellsize}")

        return {
            "data": matgr,
            "xllcorner": xllcorner,
            "yllcorner": yllcorner,
            "cellsize": cellsize,
            "crs": crs,
        }

    else:
        logger.error(f"Unsupported file format: {suffix_lower}")
        raise ValueError(f"Unsupported file format: {suffix_lower}. Must be .tif or .tiff")

Vector I/O

mobidic.preprocessing.gis_reader.read_shapefile(filepath, crs=None)

Read a shapefile and return a GeoDataFrame.

Parameters:

Name	Type	Description	Default
filepath	str | Path	Path to the shapefile (.shp).	required
crs	str | None	Optional CRS to reproject the data. If None, keeps original CRS.	None

Returns:

Type	Description
GeoDataFrame	GeoDataFrame containing the shapefile data.

Raises:

Type	Description
FileNotFoundError	If the shapefile does not exist.
Exception	If the file cannot be read by geopandas.

Examples:

>>> river_network = read_shapefile("example/shp/Arno_river_network.shp")
>>> reach_ids = river_network['REACH_ID']

Source code in mobidic/preprocessing/gis_reader.py

def read_shapefile(
    filepath: str | Path,
    crs: str | None = None,
) -> gpd.GeoDataFrame:
    """Read a shapefile and return a GeoDataFrame.

    Args:
        filepath: Path to the shapefile (.shp).
        crs: Optional CRS to reproject the data. If None, keeps original CRS.

    Returns:
        GeoDataFrame containing the shapefile data.

    Raises:
        FileNotFoundError: If the shapefile does not exist.
        Exception: If the file cannot be read by geopandas.

    Examples:
        >>> river_network = read_shapefile("example/shp/Arno_river_network.shp")
        >>> reach_ids = river_network['REACH_ID']
    """
    filepath = Path(filepath)

    # Check if file exists
    if not filepath.exists():
        logger.error(f"Shapefile not found: {filepath}")
        raise FileNotFoundError(f"Shapefile not found: {filepath}")

    logger.info(f"Reading shapefile: {filepath}")

    try:
        # Read shapefile
        gdf = gpd.read_file(filepath)

        logger.debug(f"Shapefile loaded: {len(gdf)} features, CRS={gdf.crs}")

        # Reproject if requested
        if crs is not None:
            logger.debug(f"Reprojecting from {gdf.crs} to {crs}")
            gdf = gdf.to_crs(crs)

        logger.success(f"Successfully read shapefile: {len(gdf)} features, CRS={gdf.crs}")

        return gdf

    except Exception as e:
        logger.error(f"Failed to read shapefile {filepath}: {e}")
        raise RuntimeError(f"Failed to read shapefile {filepath}: {e}") from e

Input checks

mobidic.preprocessing.gis_reader.check_crs(crs, source)

Check the CRS of an input file.

Warns when the file has no CRS and rejects geographic CRSs, since MOBIDIC assumes coordinates, lengths and cell sizes in metres.

Parameters:

Name	Type	Description	Default
crs	Any	CRS of the file (rasterio/pyproj CRS, WKT, EPSG code) or None.	required
source	str | Path	Path of the file, used in messages.	required

Returns:

Type	Description
bool	True if the file has a CRS, False otherwise.

Raises:

Type	Description
GridAlignmentError	If the CRS is geographic.

Source code in mobidic/preprocessing/gis_reader.py

def check_crs(crs: Any, source: str | Path) -> bool:
    """Check the CRS of an input file.

    Warns when the file has no CRS and rejects geographic CRSs, since MOBIDIC
    assumes coordinates, lengths and cell sizes in metres.

    Args:
        crs: CRS of the file (rasterio/pyproj CRS, WKT, EPSG code) or None.
        source: Path of the file, used in messages.

    Returns:
        True if the file has a CRS, False otherwise.

    Raises:
        GridAlignmentError: If the CRS is geographic.
    """
    parsed = parse_crs(crs)
    if parsed is None:
        logger.warning(f"{source} has no CRS: assuming it is in the same CRS as the DTM")
        return False

    if parsed.is_geographic:
        msg = (
            f"{source} has a geographic CRS ({_describe_crs(crs)}). "
            "MOBIDIC requires a projected CRS in metres: reproject the file."
        )
        logger.error(msg)
        raise GridAlignmentError(msg)

    return True

mobidic.preprocessing.gis_reader.check_raster_alignment(raster, reference, raster_path, reference_path)

Check that a raster has the same shape, origin, cell size and CRS as the reference (DTM).

No resampling or reprojection is performed: mismatching rasters must be fixed by the user.

Parameters:

Name	Type	Description	Default
raster	dict[str, Any]	Result of grid_to_matrix() for the raster to check.	required
reference	dict[str, Any]	Result of grid_to_matrix() for the reference raster (DTM).	required
raster_path	str | Path	Path of the raster to check, used in messages.	required
reference_path	str | Path	Path of the reference raster, used in messages.	required

Raises:

Type	Description
GridAlignmentError	If the raster has a geographic CRS or does not match the reference grid.

Source code in mobidic/preprocessing/gis_reader.py

def check_raster_alignment(
    raster: dict[str, Any],
    reference: dict[str, Any],
    raster_path: str | Path,
    reference_path: str | Path,
) -> None:
    """Check that a raster has the same shape, origin, cell size and CRS as the reference (DTM).

    No resampling or reprojection is performed: mismatching rasters must be fixed by the user.

    Args:
        raster: Result of grid_to_matrix() for the raster to check.
        reference: Result of grid_to_matrix() for the reference raster (DTM).
        raster_path: Path of the raster to check, used in messages.
        reference_path: Path of the reference raster, used in messages.

    Raises:
        GridAlignmentError: If the raster has a geographic CRS or does not match the reference grid.
    """
    has_crs = check_crs(raster["crs"], raster_path)

    errors = []
    shape = raster["data"].shape
    ref_shape = reference["data"].shape
    if shape != ref_shape:
        errors.append(f"shape {shape} != {ref_shape} (rows, cols)")

    cellsize = reference["cellsize"]
    if not np.isclose(raster["cellsize"], cellsize, rtol=1e-6):
        errors.append(f"cell size {raster['cellsize']} != {cellsize}")

    # Origins are compared at cell centers, with a tolerance for floating point noise
    tol = 1e-3 * cellsize
    if not (
        np.isclose(raster["xllcorner"], reference["xllcorner"], rtol=0, atol=tol)
        and np.isclose(raster["yllcorner"], reference["yllcorner"], rtol=0, atol=tol)
    ):
        errors.append(
            f"lower-left cell center ({raster['xllcorner']}, {raster['yllcorner']}) "
            f"!= ({reference['xllcorner']}, {reference['yllcorner']})"
        )

    ref_crs = reference["crs"]
    if has_crs and parse_crs(ref_crs) is not None and not crs_equals(raster["crs"], ref_crs):
        errors.append(f"CRS {_describe_crs(raster['crs'])} != {_describe_crs(ref_crs)}")

    if errors:
        msg = (
            f"Raster {raster_path} does not match the DTM grid ({reference_path}):\n  - "
            + "\n  - ".join(errors)
            + "\nAll rasters must share the DTM's extent, resolution and CRS."
        )
        logger.error(msg)
        raise GridAlignmentError(msg)

mobidic.preprocessing.gis_reader.check_vector_crs(gdf, source, reference_crs=None)

Check the CRS of a vector dataset against the reference (DTM) CRS.

Warns when the dataset has no CRS. No reprojection is performed.

Parameters:

Name	Type	Description	Default
gdf	GeoDataFrame	Vector dataset to check.	required
source	str | Path	Path of the file, used in messages.	required
reference_crs	Any	CRS of the DTM. If None, only the geographic CRS check is performed.	None

Raises:

Type	Description
GridAlignmentError	If the CRS is geographic or differs from the reference CRS.

Source code in mobidic/preprocessing/gis_reader.py

def check_vector_crs(gdf: gpd.GeoDataFrame, source: str | Path, reference_crs: Any = None) -> None:
    """Check the CRS of a vector dataset against the reference (DTM) CRS.

    Warns when the dataset has no CRS. No reprojection is performed.

    Args:
        gdf: Vector dataset to check.
        source: Path of the file, used in messages.
        reference_crs: CRS of the DTM. If None, only the geographic CRS check is performed.

    Raises:
        GridAlignmentError: If the CRS is geographic or differs from the reference CRS.
    """
    if not check_crs(gdf.crs, source):
        return

    if parse_crs(reference_crs) is not None and not crs_equals(gdf.crs, reference_crs):
        msg = (
            f"{source} has CRS {_describe_crs(gdf.crs)}, but the DTM has CRS {_describe_crs(reference_crs)}. "
            "Reproject the file to the DTM CRS."
        )
        logger.error(msg)
        raise GridAlignmentError(msg)

Examples

Reading a Raster

from mobidic import grid_to_matrix
import numpy as np

# Read a Digital Terrain Model
dtm = grid_to_matrix("path/to/dtm.tif")

# Access raster data and metadata
print(f"Shape: {dtm['data'].shape}")
print(f"Cell size: {dtm['cellsize']} m")
print(f"Lower-left corner: ({dtm['xllcorner']}, {dtm['yllcorner']})")
print(f"CRS: {dtm['crs']}")

# Access the data array (nodata converted to NaN)
elevation = dtm['data']
mean_elevation = np.nanmean(elevation)
print(f"Mean elevation: {mean_elevation:.2f} m")

Reading a shapefile

from mobidic import read_shapefile

# Read river network shapefile
network = read_shapefile("path/to/river_network.shp")

# Access as GeoDataFrame
print(f"Number of reaches: {len(network)}")
print(f"CRS: {network.crs}")
print(network.head())

# Reproject to a different CRS
network_utm = read_shapefile(
    "path/to/river_network.shp",
    crs="EPSG:32632"  # WGS84 / UTM zone 32N
)

Return values

Raster data

grid_to_matrix() returns a dictionary with:

• data (numpy.ndarray): 2D array of raster values, nodata converted to NaN (flipped vertically to match MATLAB convention)
• xllcorner (float): X coordinate of lower-left corner (cell center)
• yllcorner (float): Y coordinate of lower-left corner (cell center)
• cellsize (float): Cell size in map units
• crs (rasterio.crs.CRS): Coordinate reference system

Vector data

read_shapefile() returns a geopandas.GeoDataFrame with geometry and attribute columns.

Error handling

Both functions provide detailed error messages and logging:

• File not found: FileNotFoundError with clear message
• Invalid format: RasterioIOError or RuntimeError with details
• CRS issues: Warning when a file has no CRS; GridAlignmentError when a CRS is geographic, or when a raster or shapefile does not match the DTM grid/CRS (see the input checks above)
• Logging: All operations logged using loguru (DEBUG, INFO, SUCCESS, ERROR levels)