import geopandas as gpd
import numpy as np
from rasterio.features import rasterize
from rasterio.io import DatasetReader
from shapely.geometry import Point
from skimage.measure import label, regionprops
from .transform_utils import get_transform_with_resampling, pad_mask_to_min_size
[docs]
def mask_to_points(
mask: np.ndarray,
dataset_reader: DatasetReader,
window=None,
target_resolution_cm_per_pixel: float | None = None,
) -> gpd.GeoDataFrame:
"""
Convert a binary or labeled mask to a GeoDataFrame of points.
Each connected component (blob) in the mask is converted to a point at its
centroid using the dataset's affine transform. If the mask was resampled,
use ``target_resolution_cm_per_pixel`` to ensure correct spatial positioning.
Parameters
----------
mask : np.ndarray
A 2D array with 0 for background and any nonzero value for features.
dataset_reader : rasterio.io.DatasetReader
Rasterio dataset that the mask is derived from, used for transform and CRS.
window : rasterio.windows.Window, optional
If provided, the transform will be adjusted for the window region.
target_resolution_cm_per_pixel : float, optional
Target resolution in centimeters per pixel if the mask was resampled.
If provided and different from the dataset's original resolution, the
transform will be automatically adjusted to account for resampling.
This works for both upsampling and downsampling.
Returns
-------
gpd.GeoDataFrame
GeoDataFrame with the extracted points and their associated values (label).
Examples
--------
>>> import rasterio
>>> with rasterio.open('image.tif') as src:
... # Extract points without resampling
... mask = model.predict(src.read(1))
... gdf = mask_to_points(mask, src)
...
... # Extract points with resampling
... resampled_mask = model.predict_resampled(src.read(1))
... gdf_resampled = mask_to_points(
... resampled_mask,
... src,
... target_resolution_cm_per_pixel=3.0
... )
"""
transform = get_transform_with_resampling(
dataset_reader=dataset_reader,
window=window,
target_resolution_cm_per_pixel=target_resolution_cm_per_pixel,
)
# Label connected components in the binary mask (blobs)
binary = mask > 0
labeled = label(binary)
if labeled.max() == 0:
return gpd.GeoDataFrame(geometry=[], crs=dataset_reader.crs)
# Use original mask as intensity image to derive a representative value per blob
regions = regionprops(labeled, intensity_image=mask.astype(float))
points = []
values = []
for region in regions:
# region.centroid is (row, col) in image coordinates
r, c = region.centroid
x, y = transform * (c, r)
points.append(Point(x, y))
# Representative value: rounded mean intensity within the blob
# Works for binary (-> 1) and labeled masks (-> label value if homogeneous)
rep_val = int(round(float(region.mean_intensity))) if not np.isnan(region.mean_intensity) else 1
values.append(rep_val)
return gpd.GeoDataFrame({"value": values, "geometry": points}, crs=dataset_reader.crs)
[docs]
def points_to_mask(
gdf: gpd.GeoDataFrame,
dataset_reader: DatasetReader,
fill_value: int = 1,
background_value: int = 0,
all_touched: bool = True,
window=None,
min_size: tuple[int, int] | None = None,
) -> np.ndarray:
"""
Convert a GeoDataFrame of points to a binary mask.
Parameters
----------
gdf : gpd.GeoDataFrame
GeoDataFrame containing point geometries to be rasterized.
dataset_reader : rasterio.io.DatasetReader
Rasterio dataset reader that defines the target raster properties
(transform, CRS, height, width) for the output mask.
fill_value : int, optional
Value to assign to pixels that coincide with point geometries. Defaults to 1.
background_value : int, optional
Value to assign to pixels that do not coincide with any point. Defaults to 0.
all_touched : bool, optional
If True, all pixels touched by geometries will be burned in. Defaults to True.
window : rasterio.windows.Window, optional
If provided, only rasterize the region inside the window. The output mask will
have shape (window.height, window.width).
min_size : tuple[int, int], optional
Minimum size (height, width) for the output mask. If the generated mask is
smaller than this size, it will be padded with background_value to match
the exact minimum size. If None or if the mask is already equal to or larger
than min_size, no padding is applied. Defaults to None.
Returns
-------
np.ndarray
A 2D numpy array (mask) where pixels corresponding to input points have the
fill_value and background pixels have the background_value.
"""
if gdf.crs != dataset_reader.crs:
gdf = gdf.to_crs(dataset_reader.crs)
# Keep only point-like geometries
gdf = gdf[gdf.geometry.notnull() & gdf.geometry.type.isin(["Point", "MultiPoint"])]
if window is not None:
raster_shape = (window.height, window.width)
raster_transform = dataset_reader.window_transform(window)
else:
raster_shape = (dataset_reader.height, dataset_reader.width)
raster_transform = dataset_reader.transform
if len(gdf) == 0:
mask = np.full(raster_shape, background_value, dtype="uint8")
else:
mask = rasterize(
[(geom, fill_value) for geom in gdf.geometry],
out_shape=raster_shape,
transform=raster_transform,
fill=background_value,
dtype="uint8",
all_touched=all_touched,
)
# Apply padding if min_size is specified and mask is smaller
if min_size is not None:
mask = pad_mask_to_min_size(mask, min_size, background_value)
return mask