import warnings
from pathlib import Path
from typing import Any
import geopandas as gpd
import numpy as np
from rasterio.io import DatasetReader
from rasterio.windows import Window
from shapely.geometry.multipolygon import MultiPolygon
from shapely.geometry.polygon import Polygon
from shapelib.core.geometries.base import BaseShapes
from shapelib.tools.converters import mask_to_polygons, polygons_to_mask
[docs]
class Polygons(BaseShapes):
"""
A class for managing polygon geometries in a GeoDataFrame.
This class provides specific functionality for polygon geometries, such as calculating
areas and managing geometry-specific metadata.
Parameters
----------
data : gpd.GeoDataFrame or pd.DataFrame, optional
The input data for the GeoDataFrame, can be a GeoDataFrame or a DataFrame.
geometry : shapely.geometry.Polygon, optional
The geometry type to use (Polygon in this case). Defaults to None.
crs : str, int, or pyproj.CRS, optional
Coordinate reference system (CRS) to use for the GeoDataFrame.
path : str or pathlib.Path, optional
The file path for the data (if reading from a file).
driver : str, optional
The file driver to use when reading/writing data.
Attributes
----------
area_column_name : str
The name of the column storing the calculated area values.
geometry_type : type
The expected geometry type for this class (Polygon).
area_precision : int
The precision for rounding area values when calculating the total area.
Examples
--------
Create a Polygons object:
>>> polys = Polygons(data=gdf)
>>> total_area = polys.get_total_area() # Get total area of all polygons
>>> print(f"Total: {total_area} ha")
"""
_metadata = BaseShapes._metadata + [
"area_column_name",
"perimeter_column_name",
"geometry_type",
"area_precision",
"perimeter_precision",
]
def __init__(
self,
data: Any = None,
*args: Any,
geometry: Any = None,
crs: Any = None,
path: str | Path | None = None,
driver: str | None = None,
**kwargs: Any,
) -> None:
super().__init__(data, *args, geometry=geometry, crs=crs, path=path, driver=driver, **kwargs)
self.area_column_name: str = "area (ha)"
self.perimeter_column_name: str = "perimeter (m)"
self.geometry_type: object = Polygon
self.area_precision: int = 6
self.perimeter_precision: int = 6
# TODO: Add cache
@property
def area(self) -> gpd.GeoSeries:
"""
Return the area of each Polygon geometry (GeoPandas native behavior).
.. deprecated:: 0.1.0
Accessing `.area` now returns a GeoSeries with individual polygon areas,
following GeoPandas' native behavior. Use `.get_total_area()` to obtain
the sum of all polygon areas.
This property returns a GeoSeries containing the area of each individual
Polygon in the GeoDataFrame, computed using the current CRS. This is the
standard GeoPandas behavior.
Returns
-------
gpd.GeoSeries
A GeoSeries with the area of each geometry in the current CRS units.
See Also
--------
get_total_area : Get the sum of all polygon areas with geodetic accuracy.
calculate_area : Calculate geodetic area for each polygon and store in a column.
Examples
--------
>>> # Get individual areas (GeoSeries)
>>> areas = polys.area
>>> print(areas)
0 150.5
1 230.8
dtype: float64
>>> # Get total area of all polygons
>>> total = polys.get_total_area()
>>> print(f"Total: {total} ha")
Notes
-----
This property emits a DeprecationWarning to alert users of the behavior change.
In previous versions, `.area` returned the total sum of all polygon areas.
"""
warnings.warn(
"The `.area` property now returns a GeoSeries with individual polygon "
"areas (GeoPandas native behavior). To get the total sum of all polygon "
"areas, use `.get_total_area()` instead. This warning will be removed "
"in a future version.",
DeprecationWarning,
stacklevel=2,
)
return super().area
[docs]
def calculate_area(self, column_name: str = "area (ha)") -> None:
"""
Calculate the area for each polygon in the GeoDataFrame and store the result in a
new column.
This method transforms the GeoDataFrame to WGS84 to ensure accurate area
calculations, validates the geometries, and then calculates the area of each
polygon. The area is stored in a column with the specified name.
Parameters
----------
column_name : str, optional
The name of the column to store the calculated areas. Defaults to "area".
Returns
-------
None
The method modifies the GeoDataFrame in place by adding a new column with the
area values.
Examples
--------
>>> polys.calculate_area()
>>> print(polys["area"]) # See the calculated area for each polygon.
"""
with self.temp_crs(self.wgs84_crs):
self._validate_geoms()
self[column_name] = self.geometry.apply(self._geom_area)
self.area_column_name = column_name
[docs]
def calculate_perimeter(self, column_name: str = "perimeter (m)") -> None:
"""
Calculate the perimeter for each polygon in the GeoDataFrame and store the result
in a new column.
This method converts the GeoDataFrame to WGS84 to ensure accurate perimeter
calculations, validates the geometries, and calculates the perimeter of each
polygon. The perimeter is stored in a column with the specified name.
Parameters
----------
column_name : str, optional
The name of the column to store the calculated perimeters. Defaults to
"perimeter (m)".
Returns
-------
None
This method modifies the GeoDataFrame in place by adding a new column with
the perimeter values.
Examples
--------
>>> polys.calculate_perimeter()
>>> print(polys["perimeter (m)"]) # See the calculated perimeter for each
polygon.
"""
with self.temp_crs(self.wgs84_crs):
self._validate_geoms()
self[column_name] = self.geometry.apply(self._geom_perimeter)
self.perimeter_column_name = column_name
def _geom_area(self, geometry: Polygon | MultiPolygon) -> float:
"""
Calculate the area of a single polygon geometry.
This method uses the `Geod.geometry_area_perimeter` method to compute the area
of a given polygon geometry in square kilometers.
Parameters
----------
geometry : shapely.geometry.Polygon
The polygon geometry to calculate the area for.
Returns
-------
float
The area of the polygon in hectares.
"""
area = round(
abs(self.geod.geometry_area_perimeter(geometry)[0]) / 10_000,
self.area_precision,
)
return area
def _geom_perimeter(self, geometry: Polygon | MultiPolygon) -> float:
"""
Calculate the perimeter of a single polygon geometry.
This method uses the `Geod.geometry_area_perimeter` method to compute the
perimeter of a given polygon geometry.
Parameters
----------
geometry : shapely.geometry.Polygon
The polygon geometry to calculate the perimeter for.
Returns
-------
float
The perimeter of the polygon in meters.
"""
perimeter = round(
abs(self.geod.geometry_area_perimeter(geometry)[1]),
self.perimeter_precision,
)
return perimeter
[docs]
def get_total_perimeter(self) -> float:
"""
Return the total perimeter of all polygons in the GeoDataFrame.
This method recalculates the perimeter for each polygon, sums the values,
and rounds to the precision specified by `perimeter_precision`.
Returns
-------
float
The total perimeter of all polygons, rounded to the specified precision.
Returns 0.0 if the GeoDataFrame is empty.
Examples
--------
>>> total_perimeter = polys.get_total_perimeter()
>>> print(f"Total perimeter: {total_perimeter} meters")
"""
if len(self) == 0:
return 0.0
self.calculate_perimeter()
return round(float(self[self.perimeter_column_name].sum()), self.perimeter_precision)
[docs]
def get_total_area(self) -> float:
"""
Return the total area of all polygons in the GeoDataFrame.
This method recalculates the area for each polygon, sums the values, and rounds
to the precision specified by `area_precision`.
Returns
-------
float
The total area of all polygons, rounded to the specified precision.
Returns 0.0 if the GeoDataFrame is empty.
"""
if len(self) == 0:
return 0.0
self.calculate_area()
return round(float(self[self.area_column_name].sum()), self.area_precision)
[docs]
def force_make_valid(self) -> None:
"""
Force geometries to be valid by exploding, removing invalid/empty geometries,
dissolving, and exploding again.
This method modifies the GeoDataFrame in place.
"""
self.explode(ignore_index=True, inplace=True)
self.drop(
index=self[self.geometry.isna() | self.geometry.is_empty].index,
inplace=True,
)
# dissolve() doesn't support inplace, so we reassign
dissolved = self.dissolve()
self._update_inplace(dissolved)
self.explode(ignore_index=True, inplace=True)
self.reset_index(drop=True, inplace=True)
[docs]
def filter_by_area(self, area: float) -> None:
"""
Filter the GeoDataFrame to retain only Polygon geometries that are larger than
the specified area.
This method removes any geometries whose area is less than the specified value.
Parameters
----------
area : float
The minimum area required for a geometry to be retained. If the
value is less than or equal to zero, no filtering is performed.
Returns
-------
None
The method modifies the GeoDataFrame in place.
Examples
--------
>>> polygons.filter_by_area(500)
>>> print(f"Retained polygons larger than 500.")
"""
# Skip filtering if disabled or GeoDataFrame is empty
if area <= 0 or len(self) == 0:
return
self.calculate_area()
self.drop(self[self[self.area_column_name] < area].index, inplace=True)
[docs]
def simplify(self, tolerance: float = 0.001, **kwargs: Any) -> "Polygons":
"""
Simplify the geometries in the GeoDataFrame by reducing the number of vertices.
This method delegates to BaseShapes.simplify and invalidates area calculations
since the simplified geometries will have different areas.
Parameters
----------
tolerance : float, optional
The tolerance for simplification. Defaults to 0.001.
**kwargs : Any
Additional keyword arguments passed to BaseShapes.simplify.
Returns
-------
Polygons
The Polygons object with simplified geometries.
"""
super().simplify(tolerance, **kwargs)
return self
[docs]
def filter_by_hole_area(self, area: float) -> None:
"""
Remove interior rings (holes) whose area is smaller than the specified value, preserving all polygons.
The area threshold is expressed in hectares, consistent with `calculate_area`. This method
temporarily transforms geometries to WGS84 for accurate geodesic area computation and validates
geometries before processing. Only holes are removed; exterior boundaries are preserved.
Parameters
----------
area : float
The minimum area in hectares required for a hole to be retained. Holes with
area smaller than this value will be removed.
Returns
-------
None
The method modifies the GeoDataFrame in place.
Raises
------
InvalidGeometryError
If geometries cannot be validated after transformation to WGS84.
Examples
--------
Remove small holes (less than 1 hectare) from polygons:
>>> import shapelib
>>> polygons = shapelib.read_file("path/to/polygons.shp")
>>> print(f"Original: {len(polygons.geometry.iloc[0].interiors)} holes")
>>> polygons.filter_by_hole_area(1.0)
>>> print(f"After filtering: {len(polygons.geometry.iloc[0].interiors)} holes")
Big-O
-----
O(n * h), where n is the number of geometries and h is the average number
of holes per geometry.
Notes
-----
The method uses geodesic area calculation on the WGS84 ellipsoid for accuracy.
Area precision is 6 decimal places (0.000001 ha = 0.01 m²). The original CRS
is automatically restored after filtering.
"""
if not hasattr(self, "geometry") or not len(self.geometry):
return
with self.temp_crs(self.wgs84_crs):
self._validate_geoms()
def _filter_polygon_holes(poly: Polygon) -> Polygon:
if not poly.interiors:
return poly
kept_holes = [ring.coords for ring in poly.interiors if self._geom_area(Polygon(ring)) >= area]
return Polygon(poly.exterior.coords, kept_holes)
def _filter_holes(
geom: Polygon | MultiPolygon | None,
) -> Polygon | MultiPolygon | None:
if geom is None or geom.is_empty:
return geom
if geom.geom_type == "Polygon":
return _filter_polygon_holes(geom)
if geom.geom_type == "MultiPolygon":
parts = [
_filter_polygon_holes(part) for part in geom.geoms if part is not None and not part.is_empty
]
return MultiPolygon(parts) if parts else geom
return geom
self.geometry = self.geometry.apply(_filter_holes)
[docs]
@classmethod
def from_mask(
cls,
mask: np.ndarray[Any, Any],
dataset_reader: DatasetReader,
window: Window | None = None,
target_resolution_cm_per_pixel: float | None = None,
) -> "Polygons":
"""
Create a Polygons object from a labeled or binary mask and a rasterio dataset
reader.
This method extracts polygon geometries from a 2D NumPy array (mask), where
nonzero values represent features of interest (e.g., 1 for foreground, 0 for
background, or any labeled regions).
When the mask has been resampled to a different resolution than the original
dataset, use `target_resolution_cm_per_pixel` to ensure correct spatial
positioning of the extracted polygons.
Parameters
----------
mask : np.ndarray
A 2D NumPy array where 0 represents background and any nonzero (or labeled)
value represents features to be converted into polygons.
dataset_reader : rasterio.io.DatasetReader
An 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.
The mask should correspond to the windowed 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.
Examples:
- Original 6 cm/pixel downsampled to 3 cm/pixel: pass 3.0
- Original 3 cm/pixel upsampled to 6 cm/pixel: pass 6.0
- No resampling: omit this parameter or pass None
Returns
-------
Polygons
A Polygons object containing the extracted polygons as georeferenced
geometries.
Example
-------
>>> with rasterio.open('image.tif') as src:
... # Example 1: Without resampling
... mask = model.predict(src.read(1))
... polys = Polygons.from_mask(mask, src)
... print(polys)
...
... # Example 2: With resampling (6cm -> 3cm resolution)
... resampled_mask = model.predict_resampled(src.read(1))
... polys_resampled = Polygons.from_mask(
... resampled_mask,
... src,
... target_resolution_cm_per_pixel=3.0
... )
... print(polys_resampled)
...
... # Example 3: With window
... window = rasterio.windows.Window(0, 0, 100, 100)
... mask_window = mask[:100, :100]
... polys_window = Polygons.from_mask(mask_window, src, window=window)
... print(polys_window)
"""
gdf = mask_to_polygons(
mask,
dataset_reader,
window=window,
target_resolution_cm_per_pixel=target_resolution_cm_per_pixel,
)
return cls(gdf)
[docs]
def to_mask(
self,
dataset_reader: DatasetReader,
fill_value: int = 1,
background_value: int = 0,
all_touched: bool = True,
window: Window = None,
min_size: tuple[int, int] | None = None,
) -> np.ndarray:
"""
Convert the Polygons object to a binary mask using a rasterio dataset reader.
This method rasterizes the polygon geometries in this GeoDataFrame into a binary
mask with the same dimensions and coordinate reference system as the provided
rasterio dataset reader. If a window is provided, only the region inside the
window is rasterized. If min_size is specified and the resulting mask is smaller,
it will be padded with background_value to match the exact minimum size.
Parameters
----------
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 intersect with polygon geometries.
Defaults to 1.
background_value : int, optional
Value to assign to pixels that do not intersect with any polygon.
Defaults to 0.
all_touched : bool, optional
If True, all pixels touched by geometries will be burned in.
If False, only pixels whose center is within the polygon or that
are selected by Bresenham's line algorithm will be burned in.
Defaults to False for more precise boundary adherence.
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. Defaults to None.
Returns
-------
np.ndarray
A 2D numpy array (mask) with the same dimensions as the input raster or
window (or min_size if padding was applied), where pixels intersecting with
polygons have the fill_value and background pixels have the background_value.
Examples
--------
>>> import rasterio
>>> with rasterio.open('image.tif') as src:
... mask = polys.to_mask(src)
... print(mask.shape) # Should match src.shape
... print(mask.dtype) # uint8
...
... window = rasterio.windows.Window(0, 0, 50, 50)
... mask_window = polys.to_mask(src, window=window)
... print(mask_window.shape) # Should be (50, 50)
...
... # With minimum size padding
... mask_padded = polys.to_mask(src, window=window, min_size=(256, 256))
... print(mask_padded.shape) # Should be (256, 256)
"""
return polygons_to_mask(
self,
dataset_reader,
fill_value=fill_value,
background_value=background_value,
all_touched=all_touched,
window=window,
min_size=min_size,
)