terraflow.geo

The geo module owns ROI clipping and CRS alignment. The pipeline invariant is all output coordinates are EPSG:4326 — geo.py reprojects any input that differs.

Overview

• ROI bounds in any CRS are reprojected to the raster CRS (all four corners, then axis-aligned bounding box) before windowing.
• ROI clipping snaps requested bounds to an intersecting pixel window so very small ROIs avoid oversized raster reads.
• Degenerate windows (NaN dimensions after reprojection) raise ValueError with diagnostic context.

API Reference

geo

clip_raster_to_roi(raster, roi, roi_crs='EPSG:4326', band=1)

Clip a raster to the region of interest and return data & transform.

Parameters:

Name	Type	Description	Default
raster	DatasetReader	Open rasterio dataset.	required
roi	Dict[str, Any]	Dictionary with keys: xmin, ymin, xmax, ymax. Coordinates must be in the CRS specified by roi_crs (default EPSG:4326 / WGS 84).	required
roi_crs	str	EPSG code or WKT string for the CRS of the ROI coordinates. Defaults to "EPSG:4326" (latitude/longitude degrees). When the raster is in a different CRS the ROI is automatically reprojected before clipping.	'EPSG:4326'

Returns:

Name	Type	Description
data	MaskedArray	A masked array (band 1) containing the clipped raster values.
transform	Affine	Affine transform corresponding to the clipped window.

Raises:

Type	Description
ValueError:	If raster does not have band 1, ROI bounds are invalid, or the (reprojected) ROI does not intersect the raster extent.

Source code in terraflow/geo.py

def clip_raster_to_roi(
    raster: DatasetReader,
    roi: Dict[str, Any],
    roi_crs: str = "EPSG:4326",
    band: int = 1,
) -> Tuple[np.ma.MaskedArray, rasterio.Affine]:
    """
    Clip a raster to the region of interest and return data & transform.

    Parameters
    ----------
    raster:
        Open rasterio dataset.
    roi:
        Dictionary with keys: xmin, ymin, xmax, ymax.  Coordinates must be
        in the CRS specified by *roi_crs* (default EPSG:4326 / WGS 84).
    roi_crs:
        EPSG code or WKT string for the CRS of the ROI coordinates.
        Defaults to ``"EPSG:4326"`` (latitude/longitude degrees).  When the
        raster is in a different CRS the ROI is automatically reprojected
        before clipping.

    Returns
    -------
    data:
        A masked array (band 1) containing the clipped raster values.
    transform:
        Affine transform corresponding to the clipped window.

    Raises
    ------
    ValueError:
        If raster does not have band 1, ROI bounds are invalid, or the
        (reprojected) ROI does not intersect the raster extent.
    """
    if raster.count < 1:
        raise ValueError("Raster has no bands. Cannot read band 1.")
    if band < 1 or band > raster.count:
        raise ValueError(
            f"raster_band={band} is out of range; raster has {raster.count} band(s). "
            f"Valid values: 1..{raster.count}."
        )

    if roi["xmin"] >= roi["xmax"] or roi["ymin"] >= roi["ymax"]:
        raise ValueError(
            f"Invalid ROI bounds: xmin={roi['xmin']}, xmax={roi['xmax']}, "
            f"ymin={roi['ymin']}, ymax={roi['ymax']}. "
            f"Require xmin < xmax and ymin < ymax."
        )

    xmin, ymin, xmax, ymax = roi["xmin"], roi["ymin"], roi["xmax"], roi["ymax"]

    raster_crs = raster.crs
    src_crs = CRS.from_user_input(roi_crs)
    if src_crs != raster_crs:
        transformer = Transformer.from_crs(src_crs, raster_crs, always_xy=True)
        # Transform all 4 corners so non-linear projections are handled correctly,
        # then take the axis-aligned bounding box of the results.
        corners_x = [xmin, xmin, xmax, xmax]
        corners_y = [ymin, ymax, ymin, ymax]
        xs, ys = transformer.transform(corners_x, corners_y)
        xmin, ymin, xmax, ymax = min(xs), min(ys), max(xs), max(ys)
        logger.info(
            f"Reprojected ROI from {roi_crs} to raster CRS {raster_crs.to_epsg() or 'custom'}: "
            f"xmin={xmin:.2f} ymin={ymin:.2f} xmax={xmax:.2f} ymax={ymax:.2f}"
        )

    raster_bounds = raster.bounds
    clipped_xmin = max(xmin, raster_bounds.left)
    clipped_ymin = max(ymin, raster_bounds.bottom)
    clipped_xmax = min(xmax, raster_bounds.right)
    clipped_ymax = min(ymax, raster_bounds.top)

    if clipped_xmin >= clipped_xmax or clipped_ymin >= clipped_ymax:
        raise ValueError(
            f"ROI does not intersect the raster (raster extent: {raster.bounds}). "
            f"Verify that roi_crs='{roi_crs}' matches the coordinate system of "
            f"your xmin/ymin/xmax/ymax values."
        )

    try:
        window: Window = from_bounds(
            clipped_xmin,
            clipped_ymin,
            clipped_xmax,
            clipped_ymax,
            transform=raster.transform,
        )
    except Exception as e:
        raise ValueError(
            f"Could not compute raster window for ROI bounds "
            f"({xmin:.4f}, {ymin:.4f}, {xmax:.4f}, {ymax:.4f}): {e}. "
            f"Check roi_crs='{roi_crs}' and that the ROI overlaps the raster."
        ) from e

    # Guard against degenerate windows (NaN dimensions) that arise when
    # reprojected bounds are numerically invalid (e.g. coordinates outside
    # the valid domain of the target CRS).
    if math.isnan(window.width) or math.isnan(window.height):
        raise ValueError(
            f"ROI bounds ({xmin}, {ymin}, {xmax}, {ymax}) produced a "
            f"degenerate window after reprojection. "
            f"Check that roi_crs='{roi_crs}' is correct and that the ROI "
            f"coordinates are valid in that CRS."
        )

    # Snap to pixel boundaries and clamp to the raster extent so small ROIs
    # read only the intersecting window instead of an oversized float window.
    full_window = Window(col_off=0, row_off=0, width=raster.width, height=raster.height)
    try:
        window = window.round_offsets().round_lengths().intersection(full_window)
    except WindowError as e:
        raise ValueError(
            f"ROI does not intersect the raster (raster extent: {raster.bounds}). "
            f"Verify that roi_crs='{roi_crs}' matches the coordinate system of "
            f"your xmin/ymin/xmax/ymax values."
        ) from e

    data = raster.read(band, window=window, masked=True)
    transform = raster.window_transform(window)

    if data.size == 0:
        raise ValueError(
            f"ROI does not intersect the raster (raster extent: {raster.bounds}). "
            f"Verify that roi_crs='{roi_crs}' matches the coordinate system of "
            f"your xmin/ymin/xmax/ymax values."
        )

    return data, transform