import heapq
import logging
from concurrent.futures import ThreadPoolExecutor
from contextlib import contextmanager
from os import cpu_count
from pathlib import Path
from typing import TYPE_CHECKING, Any, ClassVar
import geopandas as gpd
import numpy as np
import pandas as pd
import shapely
from pyproj import CRS, Geod
from pyproj.aoi import AreaOfInterest
from pyproj.database import CRSInfo, query_utm_crs_info
from shapely.geometry.linestring import LineString
from shapely.geometry.point import Point
from shapely.geometry.polygon import Polygon
from shapelib.core.cloud.cloud_path import CloudPath
from shapelib.core.geometries.operations_mixin import OperationsMixin
from shapelib.core.transformations import datum_shift, transform_geometries
from shapelib.exceptions.geometries_exceptions import InvalidGeometryError
logger = logging.getLogger(__name__)
if TYPE_CHECKING:
from shapelib.core.ports.file_writer_port import FileWriterStrategy
[docs]
class BaseShapes(OperationsMixin, gpd.GeoDataFrame):
"""
A base class that extends GeoDataFrame with additional spatial and CRS-related
functionality.
This class provides utility methods for handling geometries, projections, and CRS
(Coordinate Reference System) transformations. It extends the functionality of
GeoDataFrame by incorporating additional metadata and methods for handling different
types of spatial data.
Parameters
----------
data : GeoDataFrame or DataFrame, optional
The input data for the GeoDataFrame. Can be either a GeoDataFrame or a standard
DataFrame.
geometry : LineString, Point, or Polygon, optional
The geometry type to be applied for the spatial data.
crs : CRS or str, optional
The coordinate reference system for the spatial data. Can be provided as a
`pyproj.CRS`
object or an EPSG code string.
path : str or Path, optional
Path to the file from which the data is loaded, if applicable.
driver : str, optional
The driver to use for reading the file, e.g., 'ESRI Shapefile', 'GeoJSON'.
This is useful when loading spatial data from external sources.
keep_z : bool, optional
Keep the Z coordinate of 3-D geometries. Defaults to False, which flattens
every geometry to 2-D on load. Turn it on when Z is data the caller needs,
such as the elevation carried by every vertex of a contour line.
*args : tuple
Additional positional arguments passed to `geopandas.GeoDataFrame`.
**kwargs : dict
Additional keyword arguments passed to `geopandas.GeoDataFrame`.
Attributes
----------
geometry_column_name : str, optional
The name of the active geometry column in the GeoDataFrame.
original_crs : CRS, optional
The original CRS of the data, if provided. Useful for preserving the original
projection of the dataset.
utm_crs : CRS, optional
The UTM CRS (Universal Transverse Mercator) for the dataset, if calculated.
wgs84_crs : CRS
The WGS84 CRS (EPSG:4326), used as a default for geographic data.
geod : Geod
A `pyproj.Geod` object initialized with the WGS84 ellipsoid, used for performing
geodetic calculations.
path : Path, optional
The file path from which the data was loaded or where it should be saved, if
applicable.
driver : str, optional
The driver used for reading the file, if provided.
legacy_datum : bool
True when the data is in a legacy datum (e.g. SAD69, Córrego Alegre) that
needs a datum shift to reach WGS 84. Its positional accuracy is limited by
that shift: see `datum_shift_accuracy`.
datum_shift_accuracy : float, optional
For data in a legacy datum, the published accuracy in metres of the shift
to WGS 84 (e.g. 2.0 for SAD69 with the IBGE grid); None otherwise, or when
no accuracy is published.
keep_z : bool
Whether 3-D geometries keep their Z coordinate. Carried over to the objects
derived from this one (copies, filters, reprojections).
Notes
-----
This class is designed for spatial data operations and extends GeoDataFrame with
additional capabilities such as CRS handling and custom metadata. The metadata fields,
such as `geometry_column_name` and `original_crs`, are populated during initialization
to track important spatial information throughout operations.
"""
_metadata = [
"_geometry_column_name",
"geometry_column_name",
"original_crs",
"utm_crs",
"wgs84_crs",
"geod",
"path",
"driver",
"legacy_datum",
"datum_shift_accuracy",
"keep_z",
]
# Registry of file writers (Strategy pattern)
# Writers are registered at import time, checked in order
_writers: ClassVar[list["FileWriterStrategy"]] = []
[docs]
@classmethod
def register_writer(cls, writer: "FileWriterStrategy") -> None:
"""
Register a file writer strategy.
Writers are checked in registration order. The first writer
that supports the destination path will be used.
Parameters
----------
writer : FileWriterStrategy
The writer strategy to register.
Examples
--------
>>> from shapelib.core.cloud.writers import CloudFileWriter
>>> BaseShapes.register_writer(CloudFileWriter())
"""
cls._writers.insert(0, writer)
def __init__(
self,
data: gpd.GeoDataFrame | pd.DataFrame | None = None,
*args: Any,
geometry: LineString | Point | Polygon | None = None,
crs: CRS | str | None = None,
path: str | Path | None = None,
driver: str | None = None,
keep_z: bool = False,
**kwargs: Any,
) -> None:
"""
Initialize the BaseShapes class.
This constructor initializes the base class with a set of additional parameters
specific to spatial data handling. It sets up the CRS, geometry, and file handling
attributes needed for spatial data operations.
Parameters
----------
data : GeoDataFrame or DataFrame, optional
The input data for the GeoDataFrame. Can be either a GeoDataFrame or a
standard DataFrame.
geometry : LineString, Point, or Polygon, optional
The geometry type to be applied for the spatial data.
crs : CRS or str, optional
The coordinate reference system for the spatial data. Can be provided as a
`pyproj.CRS` object or an EPSG code string.
path : str or Path, optional
Path to the file from which the data is loaded, if applicable.
driver : str, optional
The driver to use for reading the file, e.g., 'ESRI Shapefile', 'GeoJSON'.
This is useful when loading spatial data from external sources.
keep_z : bool, optional
Keep the Z coordinate of 3-D geometries instead of flattening them to
2-D (default False).
*args : tuple
Additional positional arguments passed to `geopandas.GeoDataFrame`.
**kwargs : dict
Additional keyword arguments passed to `geopandas.GeoDataFrame`.
"""
super().__init__(data, *args, geometry=geometry, crs=crs, **kwargs)
self.geometry_column_name: str | None = None
self.original_crs: CRS | None = None
self.utm_crs: CRS | None = None
self.wgs84_crs: CRS = CRS.from_epsg(4326)
self.geod: Geod = Geod(ellps="WGS84")
self.path: Path | None = Path(path) if path else None
self.driver: str | None = driver
self.legacy_datum: bool = False
self.datum_shift_accuracy: float | None = None
# Set before the initial state: that step is where Z would be dropped.
self.keep_z: bool = keep_z
self._prepare_initial_state()
# TODO: Lazy initialization, to improve performance after calling special methods
def _prepare_initial_state(self) -> None:
"""
Prepare the initial state of the BaseShapes object by performing various checks
and transformations.
This method validates the file path, checks the CRS (Coordinate Reference System),
removes invalid geometries, converts Z-coordinates (3D) to 2D if present (unless
`keep_z` is set), and estimates the UTM CRS for the dataset.
Steps performed:
- Validate the file path (if provided).
- Ensure the CRS is correctly set.
- Remove rows with invalid geometries.
- Convert geometries with Z-coordinates to 2D, unless `keep_z` is set.
- Validate the encoding of the columns.
- Check the datum: refuse data that cannot be transformed precisely, flag
data in a legacy datum.
- Estimate the UTM CRS for the data.
Returns
-------
None
This method prepares the object in-place without returning anything.
"""
self._check_path() # Ensure the path is valid.
self._check_original_crs() # Ensure the CRS is correctly set.
self.remove_missing_geometries() # Remove rows with null/empty geometries.
self._remove_z_coordinates() # Convert Z-coordinates to 2D if present.
self.check_columns_encoding() # Validate and process column encodings.
self._check_datum() # Refuse imprecise datum shifts, flag legacy datums.
self.estimate_utm_crs() # Estimate UTM CRS for the data.
if hasattr(self, "geometry"):
self.geometry_column_name: str = self.geometry.name
def _check_path(self) -> None:
"""
Validate and convert the file path to a Path object if it exists.
If a path is provided, this method ensures it is converted to a `Path` object.
This conversion simplifies file operations in subsequent methods.
Returns
-------
None
This method modifies the `self.path` attribute in place.
"""
if self.path is not None:
self.path = Path(self.path)
def _remove_z_coordinates(self) -> None:
"""
Convert geometries with Z-coordinates (3D) to 2D geometries.
This method checks if the geometries have Z-coordinates (3D). If Z-coordinates are
found, the method converts all geometries to 2D, ensuring that no Z-coordinate
information is preserved. Nothing is changed when `keep_z` is set.
Returns
-------
None
This method modifies the geometries in place by removing Z-coordinates.
"""
if getattr(self, "keep_z", False) or not hasattr(self, "geometry"):
return
if self.has_z.any():
self.geometry = self.force_2d()
def _check_original_crs(self) -> None:
"""
Validate and set the original CRS for the GeoDataFrame.
This method ensures that the CRS (Coordinate Reference System) is correctly set
for the GeoDataFrame. If the CRS is not set, it raises a ValueError. If no CRS
attribute exists (likely an empty GeoDataFrame), it sets `original_crs` to None.
Returns
-------
None
Modifies the `original_crs` attribute in place or raises a ValueError.
Raises
------
ValueError
If the GeoDataFrame has no CRS set.
"""
if not hasattr(self, "crs"): # Probably an empty geodataframe
self.original_crs = None
return
if self.crs is None:
raise ValueError("The provided data has no CRS information.")
self.original_crs = self.crs
def _check_datum(self) -> None:
"""
Refuse data that cannot be transformed precisely and flag data in a legacy datum.
Every shapelib measurement and reprojection goes through WGS 84, so the
shift from the data's datum to WGS 84 is checked where the data lies. If the
best transformation known cannot be used — it needs a grid that is not
installed, or the datum is unknown — a `DatumTransformationError` is raised
now, instead of the data being moved by metres later. If the datum needs a
shift (a legacy datum such as SAD69 or Córrego Alegre), `legacy_datum` and
`datum_shift_accuracy` are set and a warning is logged: no transformation
makes such data more accurate than the shift itself.
Raises
------
DatumTransformationError
If the data cannot be transformed to WGS 84 with the best accuracy known.
"""
if not hasattr(self, "crs") or self.crs is None or not len(self):
return
shift = datum_shift(self.crs, self.geometry.to_numpy())
if shift is None:
return
self.legacy_datum = True
self.datum_shift_accuracy = shift.accuracy
accuracy = f"{shift.accuracy:g} m" if shift.accuracy is not None else "an unpublished accuracy"
logger.warning(
"%s is in a legacy datum (%s): reaching WGS 84 needs a datum shift (%s) accurate to %s, so the "
"positional accuracy of this data cannot be better than that.",
self.path.name if self.path else "The data",
self.crs.name,
shift.operation,
accuracy,
)
[docs]
def remove_missing_geometries(self) -> None:
"""
Remove rows with missing (null or empty) geometries from the GeoDataFrame.
This method drops rows that contain null or empty geometries, and then resets
the index of the GeoDataFrame to ensure continuity of row indices.
Returns
-------
None
The method modifies the GeoDataFrame in place by removing invalid geometries.
"""
if not hasattr(self, "geometry"):
return
null_mask = self.geometry.isna()
self.drop(index=self[null_mask].index, inplace=True)
empty_mask = self.geometry.is_empty
self.drop(index=self[empty_mask].index, inplace=True)
self.reset_index(drop=True, inplace=True)
[docs]
def remove_column(self, column_name: str) -> None:
"""
Remove a column from the GeoDataFrame if it exists.
Parameters
----------
column_name : str
The name of the column to remove.
Returns
-------
None
The method modifies the GeoDataFrame in place, removing the specified column.
"""
if column_name in self.columns:
self.drop(columns=[column_name], inplace=True)
[docs]
def remove_row(self, index: int) -> None:
"""
Remove a row from the GeoDataFrame based on its index.
Parameters
----------
index : int
The index of the row to remove.
Returns
-------
None
The method modifies the GeoDataFrame in place, removing the specified row.
"""
if index in self.index:
self.drop(index=index, inplace=True)
[docs]
def remove_rows_by_column_values(self, criteria: dict[str, Any]) -> None:
"""
Remove rows where any column matches its corresponding value (OR logic).
This method filters the GeoDataFrame by removing all rows where at least one
column-value pair from the criteria dictionary matches. Multiple criteria are
combined with OR logic: a row is removed if it matches ANY of the conditions.
Both column names and string values are matched case-insensitively.
Parameters
----------
criteria : dict[str, Any]
Dictionary mapping column names to values. Rows matching any key-value
pair will be removed. Column names and string values are case-insensitive.
Columns that do not exist in the GeoDataFrame are silently ignored.
Returns
-------
None
The method modifies the GeoDataFrame in place.
Examples
--------
Remove rows where status is "inactive" OR type is "A" (case-insensitive):
>>> shapes.remove_rows_by_column_values({"STATUS": "Inactive", "Type": "a"})
Remove rows with a specific numeric value:
>>> shapes.remove_rows_by_column_values({"year": 2020})
Ignore non-existent columns (useful for flexible filtering):
>>> # 'chavesig' may not exist in all GeoDataFrames - silently skipped
>>> shapes.remove_rows_by_column_values({"chavesig": "ABC", "status": "inactive"})
Big-O
-----
O(n * k), where n is the number of rows and k is the number of criteria.
Notes
-----
The method uses OR logic: rows matching ANY criterion are removed.
For AND logic (remove only rows matching ALL criteria), chain multiple
calls or use pandas boolean indexing directly.
Case-insensitive matching applies to both column names and string values.
Numeric values are compared with exact equality.
Columns specified in criteria that do not exist in the GeoDataFrame are
silently ignored, making this method flexible for use with varying schemas.
"""
if not criteria:
return
normalized_criteria = self._normalize_criteria(criteria)
mask = self._build_removal_mask(normalized_criteria)
self.drop(index=self[mask].index, inplace=True)
self.reset_index(drop=True, inplace=True)
def _normalize_criteria(self, criteria: dict[str, Any]) -> dict[str, Any]:
"""
Normalize criteria by resolving case-insensitive column names.
Parameters
----------
criteria : dict[str, Any]
Dictionary mapping column names to values.
Returns
-------
dict[str, Any]
Dictionary with actual column names from the DataFrame.
Only includes columns that exist in the GeoDataFrame.
Notes
-----
If a column specified in criteria does not exist, it is silently ignored.
This allows for flexible filtering without requiring exact column matches.
"""
normalized = {}
for column_name, value in criteria.items():
actual_column_name = self._find_column_case_insensitive(column_name)
if actual_column_name is not None:
normalized[actual_column_name] = value
return normalized
def _find_column_case_insensitive(self, column_name: str) -> str | None:
"""
Find the actual column name in the DataFrame using case-insensitive matching.
Parameters
----------
column_name : str
The column name to search for (case-insensitive).
Returns
-------
str | None
The actual column name as it appears in the DataFrame,
or None if the column does not exist.
Notes
-----
Returns None instead of raising an error to allow flexible filtering
where not all columns need to exist in every GeoDataFrame.
"""
column_lower = column_name.lower()
for actual_column in self.columns:
if actual_column.lower() == column_lower:
return actual_column
return None
def _build_removal_mask(self, criteria: dict[str, Any]) -> pd.Series:
"""
Build a boolean mask for rows matching any criterion (OR logic).
String values are compared case-insensitively, numeric values use exact equality.
Parameters
----------
criteria : dict[str, Any]
Dictionary mapping column names to values.
Returns
-------
pd.Series
Boolean Series where True indicates rows to be removed.
"""
mask = pd.Series([False] * len(self), index=self.index)
for column_name, value in criteria.items():
column_mask = self._create_column_mask(column_name, value)
mask |= column_mask
return mask
def _create_column_mask(self, column_name: str, value: Any) -> pd.Series:
"""
Create a boolean mask for a single column-value comparison.
For string values, comparison is case-insensitive.
For other types, exact equality is used.
Parameters
----------
column_name : str
The column name to compare.
value : Any
The value to match against.
Returns
-------
pd.Series
Boolean Series where True indicates matching rows.
"""
column_data = self[column_name]
if isinstance(value, str):
# Case-insensitive comparison for strings
value_lower = value.lower()
return column_data.astype(str).str.lower() == value_lower
else:
# Exact comparison for non-string types
return column_data == value
[docs]
def add_index(self, column_name: str = "index") -> None:
"""
Add an index column to the GeoDataFrame.
This method creates a new column in the GeoDataFrame that stores the row index
values as integers, effectively adding a new index column.
Parameters
----------
column_name : str, optional
The name of the index column to add. Defaults to 'index'.
Returns
-------
None
The method modifies the GeoDataFrame in place by adding a new column with
index values.
Examples
--------
>>> gdf.add_index()
>>> gdf.head()
index geometry
0 0 POINT (1 1)
1 1 POINT (2 2)
"""
self[column_name] = range(0, len(self))
[docs]
def has_column(self, column_name: str, case_sensitive: bool = False) -> bool:
"""
Check if a column exists in the GeoDataFrame.
This method checks whether the GeoDataFrame contains a column with the specified
name, either case-sensitive or case-insensitive, depending on the parameter.
Parameters
----------
column_name : str
The name of the column to check.
case_sensitive : bool, optional
Whether the column name check should be case-sensitive. Defaults to False.
Returns
-------
bool
True if the column exists in the GeoDataFrame, otherwise False.
Examples
--------
>>> gdf.has_column("geometry")
True
>>> gdf.has_column("Geometry", case_sensitive=True)
False
"""
if not case_sensitive:
return column_name.lower() in map(str.lower, self.columns)
return column_name in self.columns
[docs]
def get_element(self, index: int) -> gpd.GeoDataFrame:
"""
Get a specific element (row) from the GeoDataFrame by index.
This method retrieves a specific row from the GeoDataFrame based on its index
and returns it as a new GeoDataFrame.
Parameters
----------
index : int
The index of the row to retrieve.
Returns
-------
gpd.GeoDataFrame
A GeoDataFrame containing the selected row.
Examples
--------
>>> gdf.get_element(1)
geometry
1 POINT (2 2)
Notes
-----
The `index` parameter corresponds to the row index of the GeoDataFrame.
"""
return self.iloc[index]
[docs]
def check_columns_encoding(self) -> None:
"""
Check and fix the encoding of column names and cell values in the GeoDataFrame.
This method identifies and fixes encoding issues with column names and individual
cell values in the GeoDataFrame. If column names or values are in `bytes` format,
the method decodes them into `str` format by first decoding from 'ascii' to hex
and then decoding the resulting string using the 'iso8859_15' encoding.
If the GeoDataFrame contains invalid indices (e.g., NaN values in the index), the
method returns without performing any further actions.
Raises
------
ValueError
If a column name or cell value cannot be decoded properly.
Examples
--------
Fix encoding issues in a GeoDataFrame:
>>> gdf.check_columns_encoding()
This method will automatically detect and correct any columns or cells in the
GeoDataFrame that are incorrectly encoded as `bytes`.
Notes
-----
This method is particularly useful for fixing encoding issues in data imported
from files that may have encoding errors, such as CSV files with mis-encoded
text fields.
"""
# Early return if index is invalid
if self.index.max() is np.nan:
return
self._fix_column_names_encoding()
self._fix_cell_values_encoding()
def _fix_bytes(self, value: bytes) -> str:
"""
Decode a byte string into a string using 'iso8859_15'.
This function first decodes the byte string from 'ascii' to hex, then decodes
the resulting string using the 'iso8859_15' encoding. If decoding fails, it
raises a `ValueError` providing the value that caused the error.
Parameters
----------
value : bytes
The byte string to decode.
Returns
-------
str
The decoded string.
Raises
------
ValueError
If decoding fails due to an invalid byte string or encoding issue.
Examples
--------
Convert bytes to string:
>>> fix_bytes(b'616263')
'abc'
This function is used within the `check_columns_encoding` method to ensure
that byte-encoded column names and values are properly decoded.
"""
try:
return bytes.fromhex(value.decode("ascii")).decode("iso8859_15")
except (ValueError, UnicodeDecodeError) as e:
raise ValueError(f"Failed to fix encoding for value {value}. Error: {e}") from e
def _fix_column_names_encoding(self) -> None:
"""
Fix encoding issues in column names of the GeoDataFrame.
This method identifies column names that are in `bytes` format and converts
them to properly encoded `str` format using the `_fix_bytes` method.
"""
# Fix column names if they are bytes
for column_name in self.columns:
if isinstance(column_name, bytes):
new_column_name = self._fix_bytes(column_name)
self.rename(columns={column_name: new_column_name}, inplace=True)
def _fix_cell_values_encoding(self) -> None:
"""
Fix encoding issues in cell values of the GeoDataFrame.
This method identifies cell values that are in `bytes` format and converts
them to properly encoded `str` format using the `_fix_bytes` method. It
processes each column and each row in the GeoDataFrame.
"""
# Skip processing if index is invalid
if self.index.max() is np.nan:
return
# Fix each cell's encoding if it contains bytes
for column_name in self.columns:
for index in self.index:
value = self.loc[index, column_name]
if isinstance(value, bytes):
self.loc[index, column_name] = self._fix_bytes(value)
[docs]
def multi_to_single_geometries(self) -> "BaseShapes":
"""
Convert multi-part geometries (e.g., MultiPolygon, MultiLineString) into single
geometries.
This method uses the 'explode' function to split multi-part geometries into
individual geometries and returns a new BaseShapes with each part as a separate
geometry.
Returns
-------
BaseShapes
A new BaseShapes with multi-part geometries split into single geometries.
Examples
--------
>>> gdf.multi_to_single_geometries()
The method will convert geometries like MultiPolygons into separate Polygon
objects.
Notes
-----
The `explode` method used in this function will create a new row for each
individual geometry, allowing for easier manipulation of each part of multi-part
geometries.
"""
return self.explode(
column=self.geometry_column_name,
ignore_index=True,
index_parts=False,
)
[docs]
def intersect_with_geom(
self,
geom: Point | Polygon | LineString,
predicate: str = "intersects",
sort: bool = False,
distance: float | None = None,
) -> gpd.GeoDataFrame:
"""
Query the spatial index to find geometries that satisfy the given predicate with
respect to a provided geometry.
This method uses the spatial index (sindex) of the GeoDataFrame to query for
geometries that satisfy the specified spatial predicate (e.g., 'intersects',
'within', 'contains'). It returns a subset of the GeoDataFrame containing only
those geometries that meet the criteria when compared to the input geometry.
Parameters
----------
geom : Point, Polygon, or LineString
The geometry to query against.
predicate : str, optional
The spatial predicate to use for the query. Defaults to 'intersects'.
sort : bool, optional
Whether to sort the results by distance from the input geometry.
Defaults to False.
distance : float, optional
A distance parameter used for certain predicates (e.g., 'dwithin').
Defaults to None.
Returns
-------
gpd.GeoDataFrame
A subset of the GeoDataFrame containing geometries that satisfy the predicate.
Raises
------
ValueError
If the spatial index (sindex) is not available or the query fails.
Examples
--------
>>> gdf.intersect_with_geom(Point(1, 1), predicate="contains")
This will return all geometries in the GeoDataFrame that contain the Point(1, 1).
Notes
-----
This method relies on the spatial index (sindex) being available. If the sindex
is missing, a `ValueError` is raised. This method is useful for efficient spatial
queries.
"""
try:
return self.iloc[self.sindex.query(geometry=geom, predicate=predicate, sort=sort, distance=distance)]
except AttributeError as e:
raise ValueError("Spatial index (sindex) not available for this GeoDataFrame.") from e
[docs]
def estimate_utm_crs(self) -> None:
"""
Estimate the UTM CRS for the geometries in the GeoDataFrame.
This method estimates the appropriate UTM (Universal Transverse Mercator) CRS for
the geometries in the GeoDataFrame by querying possible UTM CRS candidates based
on the total bounds (bounding box) of the geometries. The method:
- Converts the geometries to WGS84 for accurate CRS estimation.
- Determines the appropriate datum based on the geometry's location.
- Queries the UTM CRS candidates for the bounding box using the selected datum.
- Automatically sets the UTM CRS if exactly one candidate is found.
- If 2 to 4 candidates are found (e.g., data on a zone border or across the
equator), lists them and picks the zone holding most of the data (see
`_dominant_utm_candidate`).
- Prints an error if no UTM CRS can be determined or if the geometries cover too
many zones.
After estimation, the method reverts the GeoDataFrame back to its original CRS.
Returns
-------
None
The method sets the `utm_crs` attribute or prints messages if estimation
fails.
Raises
------
ValueError
If no appropriate datum can be determined for the geometry's location or
if no UTM CRS can be determined for the area.
Notes
-----
- This method depends on `pyproj.query_utm_crs_info` for UTM CRS estimation.
- The method automatically selects the most appropriate datum based on the
geometry's location.
- Supported datums include: SIRGAS 2000 (South America), NAD83(2011) (USA),
NAD83(PA11) (Hawaii), NAD83(CSRS) (Canada), ITRF2008 (Mexico).
"""
if pd.isna(self.index.max()) or not hasattr(self, "geometry") or not len(self.index):
return
# Convert to WGS84 for accurate CRS estimation
self.to_wgs84()
# Get bounds for datum selection
min_x, min_y, max_x, max_y = self.total_bounds
centroid_lon = (min_x + max_x) / 2
centroid_lat = (min_y + max_y) / 2
try:
# Select appropriate datum based on geometry location
datum_name = self._select_appropriate_datum(centroid_lon, centroid_lat)
# Query possible UTM CRS candidates based on the total bounds (bounding box)
utm_candidates = query_utm_crs_info(
datum_name=datum_name,
area_of_interest=AreaOfInterest(*self.total_bounds),
)
except ValueError as e:
# Revert to the original CRS before raising the error
self.to_original_crs()
raise ValueError(f"Failed to determine appropriate datum: {e}") from e
# Check the number of UTM CRS candidates
n_results = len(utm_candidates)
if n_results == 0:
self.to_original_crs()
raise ValueError(
f"No UTM CRS found for the provided geometries using datum "
f"'{datum_name}'."
"The geometries may be outside the coverage area of this datum."
)
elif n_results == 1:
# Set the UTM CRS if exactly one candidate is found
self.utm_crs = CRS.from_epsg(utm_candidates[0].code)
elif 1 < n_results < 5:
print(f"The provided geometries are covering {n_results} UTM zones:")
for index, crs in enumerate(utm_candidates, start=1):
print(f"Option {index}: {crs.name}, EPSG:{crs.code}")
dominant = self._dominant_utm_candidate(utm_candidates)
self.utm_crs = CRS.from_epsg(dominant.code)
print(f"Picking {dominant.name}, EPSG:{dominant.code}, the zone holding most of the data.")
else:
print(f"The provided geometries are too large, covering more than {n_results} zones.")
print("Could not automatically determine an UTM CRS. Please provide it manually.")
self.utm_crs = None
# Revert to the original CRS after estimation
self.to_original_crs()
def _dominant_utm_candidate(self, candidates: list[CRSInfo]) -> CRSInfo:
"""
Pick the UTM candidate whose zone holds most of the data.
Each geometry is assigned to the UTM zone of its representative point (the
zone number from the longitude, the hemisphere from the sign of the latitude)
and weighted by its size: area for polygons, length for lines, one for points.
The candidate whose zone has the largest total weight wins. Ties, and zones
holding no data, keep the candidates' order, so the first one is still picked
when the data gives no answer.
Must be called while the data is in WGS84. Areas and lengths are taken in
degrees: they only rank zones that lie side by side, so the planar
approximation does not change the outcome.
Parameters
----------
candidates : list[CRSInfo]
The UTM CRSs returned by `pyproj.database.query_utm_crs_info`.
Returns
-------
CRSInfo
The candidate covering most of the data.
Examples
--------
A field mostly in zone 23S, crossing the 48°W border into zone 22S:
>>> polygons = Polygons(gpd.GeoDataFrame(geometry=[box(-48.2, -21.3, -46.5, -21.0)], crs=4326))
>>> polygons.utm_crs.to_epsg()
31983
Big-O
-----
O(n) in the number of geometries, plus O(c) CRS lookups for the c candidates.
"""
geometries = self.geometry.to_numpy()
points = shapely.point_on_surface(geometries)
lon, lat = shapely.get_x(points), shapely.get_y(points)
located = ~(np.isnan(lon) | np.isnan(lat))
areas, lengths = shapely.area(geometries), shapely.length(geometries)
if np.nansum(areas) > 0:
weights = areas
elif np.nansum(lengths) > 0:
weights = lengths
else:
weights = np.ones(len(geometries))
zone_numbers = np.clip(np.floor((lon[located] + 180) / 6).astype(int) + 1, 1, 60)
hemispheres = np.where(lat[located] >= 0, "N", "S")
zones = np.char.add(zone_numbers.astype(str), hemispheres)
weight_per_zone = pd.Series(weights[located]).groupby(zones).sum()
# max() returns the first of equal weights, so ties keep the candidates' order
return max(
candidates,
key=lambda candidate: weight_per_zone.get(CRS.from_epsg(candidate.code).utm_zone, 0.0),
)
def _select_appropriate_datum(self, lon: float, lat: float) -> str:
"""
Select the most appropriate datum based on the geometry's location.
This method determines which geographic datum is most appropriate for the
given coordinates by checking which region they fall within. The evaluation
follows a priority order to resolve overlapping regions.
Parameters
----------
lon : float
Longitude of the geometry's centroid.
lat : float
Latitude of the geometry's centroid.
Returns
-------
str
Name of the most appropriate datum for the location.
Raises
------
ValueError
If the location doesn't fall within any of the defined regions with
specialized datums.
Notes
-----
The method selects from the following datums based on location priority:
1. Mexico ITRF2008: Mexico territories
2. NAD83(PA11): Hawaii, American Samoa and Pacific islands
3. NAD83(2011): Continental USA, Alaska, Puerto Rico, Virgin Islands
4. NAD83(CSRS) v8: Canada
5. SIRGAS 2000: South America and Central America
"""
# Define datum regions with official EPSG bounds and test functions
# The order here defines priority (earlier entries take precedence)
datum_regions = [
# Mexico ITRF2008 (EPSG:6365)
{
"name": "Mexico ITRF2008",
"epsg": 6365,
"datum": "ITRF2008",
"bounds": [-122.19, 12.10, -87.5, 32.72],
"test": lambda lon, lat: -122.19 <= lon <= -87.5 and 12.10 <= lat <= 32.72,
},
# NAD83(PA11) (EPSG:6320) - Pacific Plate (Hawaii, American Samoa)
# This region crosses the International Date Line
{
"name": "NAD83(PA11)",
"epsg": 6320,
"datum": "NAD83(PA11)",
"bounds": [157.47, -17.56, -151.27, 31.80],
"test": lambda lon, lat: ((lon >= 157.47) or (lon <= -151.27)) and -17.56 <= lat <= 31.80,
},
# NAD83(2011) (EPSG:6318) - USA, Alaska, Puerto Rico, Virgin Islands
# This region crosses the International Date Line
{
"name": "NAD83(2011)",
"epsg": 6318,
"datum": "NAD83(2011)",
"bounds": [167.65, 14.92, -63.88, 74.71],
"test": lambda lon, lat: ((lon >= 167.65) or (lon <= -63.88)) and 14.92 <= lat <= 74.71,
},
# NAD83(CSRS) v8 (EPSG:10413) - Canada
{
"name": "NAD83(CSRS) v8",
"epsg": 10413,
"datum": "NAD83(CSRS) v8",
"bounds": [-141.01, 38.21, -40.73, 86.46],
"test": lambda lon, lat: -141.01 <= lon <= -40.73 and 38.21 <= lat <= 86.46,
},
# SIRGAS 2000 (EPSG:4674) - South America and Central America
{
"name": "SIRGAS 2000",
"epsg": 4674,
"datum": "SIRGAS 2000",
"bounds": [-122.19, -59.87, -25.28, 32.72],
"test": lambda lon, lat: -122.19 <= lon <= -25.28 and -59.87 <= lat <= 32.72,
},
# ETRS89 (EPSG:4258) – Europa continental e off-shore
{
"name": "ETRS89",
"epsg": 4258,
"datum": "ETRS89",
"bounds": [-16.10, 33.26, 38.01, 84.73],
"test": lambda lon, lat: -16.10 <= lon <= 38.01 and 33.26 <= lat <= 84.73,
},
]
# Check each region in priority order
for region in datum_regions:
if region["test"](lon, lat):
return region["datum"]
# If no matching region is found, raise an error with details
datum_list = ", ".join([f"{r['name']} (EPSG:{r['epsg']})" for r in datum_regions])
raise ValueError(
f"No appropriate datum found for location (lon={lon}, lat={lat}). "
f"The geometry is outside all defined regions with specialized datums: "
f"{datum_list}."
)
[docs]
def save(
self,
path: str | Path | None = None,
driver: str | None = None,
schema: Any | None = None,
index: Any | None = None,
**kwargs: Any,
) -> None:
"""
Save the GeoDataFrame to local file or remote cloud storage.
This method automatically detects whether the path is local or remote
and handles the save operation accordingly. For remote paths, the file
is first saved locally then uploaded to the cloud destination.
Parameters
----------
path : str | Path | None
The destination path. Can be:
- A local file path (string or Path object)
- An HTTPS URL (presigned PUT URL for S3)
- An S3 URL (s3://bucket/key)
If not provided, uses the instance's stored path.
driver : str | None
The file driver (e.g., 'ESRI Shapefile', 'FlatGeobuf').
schema : Any | None
The schema definition for the output file.
index : Any | None
Whether to write the DataFrame index to the file.
**kwargs : Any
Additional arguments passed to `to_file`.
Raises
------
FileNotFoundError
If no path is provided and no instance path exists.
ConnectionError
If upload to remote storage fails.
PermissionError
If access to remote destination is denied.
Returns
-------
None
Saves the GeoDataFrame and returns nothing.
Examples
--------
Save to a local file:
>>> gdf.save(path="output.shp", driver="ESRI Shapefile")
Save to S3 using s3:// protocol (requires AWS credentials):
>>> gdf.save(path="s3://my-bucket/output/data.fgb")
Save using a presigned S3 PUT URL:
>>> gdf.save(path="https://bucket.s3.amazonaws.com/file.fgb?...")
Big-O
-----
O(n) where n is data size. Remote saves add network latency.
"""
if not path and not self.path:
raise FileNotFoundError("You need to specify a path to save the file.")
target_path = str(path) if path else str(self.path)
cloud_path = CloudPath.from_string(target_path)
# Try registered writers first (Strategy pattern)
for writer in self._writers:
if writer.supports(cloud_path):
writer.write(self, cloud_path, driver, schema, index, **kwargs)
return
# Default: save to local filesystem
self._save_to_local(Path(target_path), driver, schema, index, **kwargs)
def _save_to_local(
self,
path: Path,
driver: str | None,
schema: Any | None,
index: Any | None,
**kwargs: Any,
) -> None:
"""Save to a local file path (default fallback)."""
self.path = path
path.parent.mkdir(parents=True, exist_ok=True)
self.to_file(filename=path, driver=driver, schema=schema, index=index, **kwargs)
[docs]
def save_empty_file(self) -> None:
"""
Placeholder method to save an empty file (e.g., a blank shapefile).
This method is a placeholder and is not yet implemented. When called, it raises a
`NotImplementedError` to indicate that saving an empty file is not currently
supported.
Raises
------
NotImplementedError
Always, as saving an empty file is not yet implemented.
Notes
-----
The implementation could potentially save a blank shapefile or a similar empty
file format in the future.
"""
# TODO: Implement a way to save an empty file, possibly a blank shapefile.
raise NotImplementedError("Saving an empty file is not yet implemented.")
[docs]
def to_crs(self, crs: Any = None, epsg: int | None = None, inplace: bool = False) -> "BaseShapes | None":
"""
Transform the geometries to a new CRS with the best transformation known, or fail.
Same contract as `geopandas.GeoDataFrame.to_crs`, with one difference that
matters for precision. geopandas silently falls back to a less accurate
transformation when the best one needs a grid that is not installed, or to a
"ballpark" one that ignores the datum shift when the datum is unknown — moving
the data by metres, or tens of metres — and turns coordinates that fail to
transform into ``inf``. This method only uses the best transformation and
checks every coordinate: if that cannot be done, it raises and leaves the
data unchanged. Brazil's legacy datums (SAD69, SAD69(96), Córrego Alegre 1961
and 1970-72) use the IBGE grids bundled with shapelib.
Every shapelib method that reprojects (`to_utm`, `to_wgs84`, `temp_crs`, the
area and length calculations...) goes through here.
Parameters
----------
crs : pyproj.CRS or anything `CRS.from_user_input` accepts, optional
The target CRS. Either ``crs`` or ``epsg`` must be given.
epsg : int, optional
The target CRS as an EPSG code.
inplace : bool, default False
Whether to transform this object instead of returning a transformed copy.
Returns
-------
BaseShapes or None
The transformed copy, or None when ``inplace`` is True.
Raises
------
ValueError
If the data has no CRS, or neither ``crs`` nor ``epsg`` is given.
DatumTransformationError
If the best transformation cannot be used: it needs a grid that is not
installed, or there is no known transformation between the datums.
Examples
--------
>>> lines = read_file("rows_sad69.shp") # SAD69 / UTM zone 23S
>>> lines_sirgas = lines.to_crs("EPSG:31983") # IBGE grid, accurate to 1 m
Big-O
-----
O(V) in the total number of vertices.
"""
if self.crs is None:
raise ValueError("Cannot transform naive geometries. Please set a CRS on the object first.")
if crs is not None:
target = CRS.from_user_input(crs)
elif epsg is not None:
target = CRS.from_epsg(epsg)
else:
raise ValueError("Must pass either crs or epsg.")
geometry = self.geometry
if not self.crs.is_exact_same(target):
transformed = transform_geometries(geometry.to_numpy(), self.crs, target)
geometry = gpd.GeoSeries(transformed, index=geometry.index, crs=target, name=geometry.name)
frame = self if inplace else self.copy()
frame.geometry = geometry
return None if inplace else frame
[docs]
def to_utm(self) -> None:
"""
Transform the GeoDataFrame to its estimated UTM (Universal Transverse Mercator)
CRS.
This method transforms the GeoDataFrame to the UTM CRS, which is estimated based
on the geometries in the dataset. If the UTM CRS is not already set, the method
will estimate it. If the GeoDataFrame is empty, this is a no-op with a warning.
Raises
------
ValueError
If the UTM CRS cannot be determined or estimated on a non-empty GeoDataFrame.
Examples
--------
>>> gdf.to_utm()
The GeoDataFrame is transformed in place to its appropriate UTM CRS.
"""
if len(self) == 0:
print("Warning: Cannot convert empty GeoDataFrame to UTM. Skipping.")
return
if self.utm_crs is None:
self.estimate_utm_crs()
if self.utm_crs is None:
raise ValueError("UTM CRS could not be determined.")
self.to_crs(crs=self.utm_crs, inplace=True)
[docs]
def to_wgs84(self) -> None:
"""
Transform the GeoDataFrame to the WGS84 CRS (EPSG:4326).
This method transforms the GeoDataFrame to the WGS84 CRS (World Geodetic System
1984), which is commonly used for global geographic coordinates (latitude and
longitude). The transformation is performed in place, modifying the
GeoDataFrame's CRS.
Returns
-------
None
The method modifies the GeoDataFrame in place.
Examples
--------
>>> gdf.to_wgs84()
The GeoDataFrame is transformed to WGS84 (EPSG:4326).
"""
self.to_crs(crs=self.wgs84_crs, inplace=True)
[docs]
def to_original_crs(self) -> None:
"""
Revert the GeoDataFrame to its original CRS.
This method transforms the GeoDataFrame back to its original CRS, if it was
defined earlier.
If the original CRS is not defined, a `ValueError` is raised.
Raises
------
ValueError
If the original CRS is not defined.
Returns
-------
None
The method modifies the GeoDataFrame in place.
Examples
--------
>>> gdf.to_original_crs()
The GeoDataFrame is transformed back to its original CRS, if available.
"""
if self.original_crs is not None:
self.to_crs(crs=self.original_crs, inplace=True)
else:
raise ValueError("Original CRS is not defined.")
[docs]
@contextmanager
def temp_crs(self, crs: str | int | CRS) -> "BaseShapes":
"""
Temporarily change the CRS (Coordinate Reference System) of the GeoDataFrame.
This method allows you to temporarily set a new CRS for the GeoDataFrame and
automatically revert back to the original CRS after the context block is exited.
This is useful for performing operations in a specific CRS without permanently
altering the GeoDataFrame's CRS.
Parameters
----------
crs : str, int, or pyproj.CRS
The new CRS to set, which can be specified as a string (e.g., "EPSG:4326"),
an integer representing an EPSG code, or a pyproj.CRS object.
Yields
------
BaseShapes
The GeoDataFrame with the temporarily set CRS.
Raises
------
pyproj.exceptions.CRSError
If the specified CRS is invalid or cannot be applied.
Examples
--------
Temporarily switch the CRS to WGS84:
>>> with geodataframe.temp_crs("EPSG:4326"):
... # Perform operations with the new CRS
... # CRS will automatically revert back to the original after this block
"""
original_crs = self.crs # Store the original CRS
try:
self.to_crs(crs, inplace=True) # Change to the new CRS
yield self # Yield the object with the new CRS
finally:
self.to_crs(original_crs, inplace=True) # Revert to the original CRS
def _validate_geoms(self) -> None:
"""
Validate and fix invalid geometries in the GeoDataFrame.
This method checks whether all geometries are valid. If any geometries
are found to be invalid (e.g., self-intersections, invalid ring order),
it attempts to fix them using Shapely's `make_valid` with the 'structure'
method. If invalid geometries still exist after this attempt, an
`InvalidGeometryError` is raised.
Raises
------
InvalidGeometryError
If there are geometries that remain invalid after attempting to fix them.
Examples
--------
Fix invalid geometries before spatial operations:
>>> polygons = read_file("data.shp")
>>> polygons._validate_geoms()
>>> result = gpd.overlay(polygons, mask, how="intersection")
Notes
-----
Call this method before performing spatial operations like overlay, clip,
or intersection if you suspect your data may contain invalid geometries.
Invalid geometries commonly occur in data from external sources (QGIS,
ArcGIS, GPS devices) or after CRS transformations due to floating-point
precision issues.
The 'structure' method reasons from input geometry structure, treating
exterior rings as area boundaries and interior holes as exclusions.
Big-O
-----
O(n * v) where n is number of geometries and v is average vertices per geometry.
"""
if not hasattr(self, "geometry") or len(self) == 0:
return
if self.is_valid.all():
return
# First attempt: make_valid with structure method
self.geometry = self.make_valid(method="structure", keep_collapsed=False)
if self.is_valid.all():
return
# Fallback: buffer(0) for geometries that make_valid couldn't fix
# This is a classic technique for fixing self-intersections
self.geometry = self.geometry.buffer(0)
if not self.is_valid.all():
raise InvalidGeometryError("There are some invalid geometries, and we could not validate them.")
[docs]
def filter_by_geom_type(self, geom_types: str | list[str]) -> "BaseShapes":
"""
Filter the GeoDataFrame by geometry type (e.g., 'Point', 'Polygon').
Parameters
----------
geom_types : list[str]
The geometry type to filter by (e.g., 'Point', 'Polygon').
Returns
-------
BaseShapes
A new GeoDataFrame with geometries filtered by the specified type.
"""
if isinstance(geom_types, str):
geom_types = [geom_types]
return self[self.geometry.geom_type.isin(geom_types)]
_PARALLEL_SIMPLIFY_THRESHOLD: ClassVar[int] = 500
[docs]
def simplify(self, tolerance: float = 0.001, **kwargs: Any) -> "BaseShapes":
"""
Simplify geometries by reducing vertex count.
Geometries are temporarily projected to UTM for metric accuracy.
When the geometry count exceeds ``_PARALLEL_SIMPLIFY_THRESHOLD``,
simplification runs in parallel threads (GEOS releases the GIL).
Parameters
----------
tolerance : float, optional
Maximum distance between original and simplified geometry
(in UTM metres). Defaults to 0.001.
**kwargs : Any
Forwarded to ``shapely.simplify`` (e.g. ``preserve_topology``).
Returns
-------
BaseShapes
The same object with simplified geometries.
Examples
--------
>>> shapes = BaseShapes(gdf)
>>> shapes.simplify(tolerance=0.01)
Big-O
-----
O(n * v) where *n* = geometry count, *v* = avg vertices per geometry.
"""
if self.utm_crs is None:
self.estimate_utm_crs()
if not len(self.geometry):
return self
with self.temp_crs(self.utm_crs):
if len(self.geometry) < self._PARALLEL_SIMPLIFY_THRESHOLD:
self.geometry = super().simplify(tolerance, **kwargs)
else:
self.geometry = self._parallel_simplify(tolerance, **kwargs)
return self
def _parallel_simplify(self, tolerance: float, **kwargs: Any) -> gpd.GeoSeries:
"""
Simplify geometries across multiple threads with vertex-weight balancing.
Chunks are built using LPT (Longest Processing Time) scheduling so
that heavy polygons land on separate workers instead of clustering
together. Each thread gets its own GEOS context (thread-local in
Shapely) and the GIL is released during the C-level call.
Parameters
----------
tolerance : float
Simplification tolerance in current CRS units.
**kwargs : Any
Forwarded to ``shapely.simplify``.
Returns
-------
gpd.GeoSeries
Simplified geometries with the same index and CRS.
Big-O
-----
O(n log w) for scheduling (n geometries, w workers) + O(V) for
simplification where V = total vertices.
"""
preserve_topology: bool = kwargs.get("preserve_topology", True)
geom_array: np.ndarray = self.geometry.to_numpy()
num_workers: int = min(cpu_count() or 1, len(geom_array))
buckets: list[list[int]] = self._build_balanced_buckets(geom_array, num_workers)
logger.info(
"Parallel simplify | geometries=%d workers=%d tolerance=%.4f preserve_topology=%s",
len(geom_array),
num_workers,
tolerance,
preserve_topology,
)
def _simplify_bucket(indices: list[int]) -> tuple[list[int], np.ndarray]:
chunk = geom_array[indices]
simplified = shapely.simplify(chunk, tolerance, preserve_topology=preserve_topology)
return indices, simplified
result_array: np.ndarray = np.empty(len(geom_array), dtype=object)
with ThreadPoolExecutor(max_workers=num_workers) as pool:
for indices, simplified in pool.map(_simplify_bucket, buckets):
for position, idx in enumerate(indices):
result_array[idx] = simplified[position]
return gpd.GeoSeries(
result_array,
index=self.geometry.index,
crs=self.crs,
)
@staticmethod
def _build_balanced_buckets(
geom_array: np.ndarray,
num_workers: int,
) -> list[list[int]]:
"""
Assign geometry indices to worker buckets balanced by vertex weight.
Uses a greedy LPT (Longest Processing Time first) algorithm: sort
geometries by coordinate count descending, then always assign the
next-heaviest geometry to the least-loaded worker.
Parameters
----------
geom_array : np.ndarray
Array of Shapely geometry objects.
num_workers : int
Number of worker threads.
Returns
-------
list[list[int]]
One list of geometry indices per worker.
Big-O
-----
O(n log w) where n = geometry count, w = num_workers.
"""
weights: np.ndarray = shapely.get_num_coordinates(geom_array)
order: np.ndarray = np.argsort(weights)[::-1]
heap: list[tuple[int, int]] = [(0, i) for i in range(num_workers)]
heapq.heapify(heap)
buckets: list[list[int]] = [[] for _ in range(num_workers)]
for idx in order:
load, worker_id = heapq.heappop(heap)
buckets[worker_id].append(int(idx))
heapq.heappush(heap, (load + int(weights[idx]), worker_id))
non_empty = [b for b in buckets if b]
total_weights = [sum(int(weights[i]) for i in b) for b in non_empty]
if total_weights:
logger.info(
"Bucket balance | max_weight=%d min_weight=%d ratio=%.2f",
max(total_weights),
min(total_weights),
max(total_weights) / max(min(total_weights), 1),
)
return non_empty