import logging
from pathlib import Path
import geopandas as gpd
import numpy as np
from shapely.geometry import LineString, MultiLineString, Point
logger = logging.getLogger(__name__)
# EPSG code assumed for files that carry no CRS but whose coordinates are lon/lat.
WGS84_EPSG = 4326
[docs]
def multiline_to_single_line(
geometry: LineString | MultiLineString,
) -> LineString:
"""
Flatten a MultiLineString into a single LineString.
The vertices of every part are concatenated, in part order, into one
continuous ``LineString``. If the geometry is already a ``LineString`` it is
returned unchanged.
Parameters
----------
geometry : shapely.geometry.LineString or shapely.geometry.MultiLineString
The geometry to flatten.
Returns
-------
shapely.geometry.LineString
A single LineString built from all the input vertices.
Examples
--------
>>> from shapely.geometry import MultiLineString
>>> mls = MultiLineString([[(0, 0), (1, 0)], [(1, 0), (2, 0)]])
>>> multiline_to_single_line(mls).wkt
'LINESTRING (0 0, 1 0, 1 0, 2 0)'
"""
if isinstance(geometry, LineString):
return geometry
flat_coords = [Point(*point) for part in geometry.geoms for point in part.coords]
return LineString(flat_coords)
[docs]
def get_most_common_geometry_type(gdf: str | Path | gpd.GeoDataFrame) -> str | None:
"""
Identify the most common geometry type in the given shapefile or GeoDataFrame.
This function analyzes the geometries in a GeoDataFrame or a shapefile and identifies
the most frequent geometry type (e.g., "Point", "LineString", "Polygon"). It ignores
invalid or missing geometries during the analysis.
Parameters
----------
gdf : Union[str, Path, gpd.GeoDataFrame]
The path to a shapefile (as a string or Path) or an existing GeoDataFrame.
If a file path is provided, it will be loaded as a GeoDataFrame.
Returns
-------
Union[str, None]
The most common geometry type in the input GeoDataFrame or shapefile. If no valid
geometries are found or the input is invalid, it returns `None`.
Examples
--------
>>> gdf = gpd.read_file("/path/to/shapefile.shp")
>>> most_common_geom = get_most_common_geometry_type(gdf)
>>> print(most_common_geom)
'Polygon'
Notes
-----
The function only processes valid geometries and skips invalid or missing entries.
It uses the `geom_type` property of GeoDataFrame geometries to identify the type.
"""
# Filter out invalid geometries
valid_geometries = gdf.geometry.dropna()
# If there are no valid geometries, return None
if valid_geometries.empty:
print("No valid geometries found in the input!")
return None
# Use value_counts to find the most common geometry type
return valid_geometries.geom_type.value_counts().idxmax()
[docs]
def assume_wgs84_if_lonlat(gdf: gpd.GeoDataFrame, source: str) -> gpd.GeoDataFrame:
"""
Assign WGS84 (EPSG:4326) to data read without a CRS when its coordinates are lon/lat.
A file can come without CRS information — a shapefile whose ``.prj`` is
missing or empty, a FlatGeobuf or GeoPackage written without one. When every
coordinate fits the lon/lat range (x in [-180, 180], y in [-90, 90]) the data
is taken as WGS84 and a warning is logged; otherwise it cannot be placed on
the map and a ``ValueError`` is raised. Data with no coordinates at all (no
non-empty geometry) is taken as WGS84 too, since there is nothing to misplace.
This is a reading policy: it belongs to the file handlers, not to
``BaseShapes``, which keeps rejecting in-memory data built without a CRS.
Parameters
----------
gdf : gpd.GeoDataFrame
The data as read from the file.
source : str
The file the data was read from, used in the log and error messages.
Returns
-------
gpd.GeoDataFrame
``gdf`` unchanged if it already has a CRS (or no active geometry column),
otherwise a copy with EPSG:4326 assigned.
Raises
------
ValueError
If the data has no CRS and its coordinates do not fit the lon/lat range.
Examples
--------
>>> gdf = gpd.read_file("prescription.zip") # shapefile without .prj
>>> gdf.crs is None
True
>>> assume_wgs84_if_lonlat(gdf, "prescription.zip").crs.to_epsg()
4326
"""
if not hasattr(gdf, "crs") or gdf.crs is not None:
return gdf
min_x, min_y, max_x, max_y = gdf.total_bounds
if np.isnan([min_x, min_y, max_x, max_y]).all():
logger.warning("%s has no CRS and no coordinates; assuming EPSG:%d.", source, WGS84_EPSG)
return gdf.set_crs(epsg=WGS84_EPSG)
if not (min_x >= -180 and max_x <= 180 and min_y >= -90 and max_y <= 90):
raise ValueError(
f"{source} has no CRS and its coordinates are not lon/lat "
f"(bounds: {min_x:.2f}, {min_y:.2f}, {max_x:.2f}, {max_y:.2f}), so it cannot be "
"placed on the map. Provide its CRS (for a shapefile, the .prj file)."
)
logger.warning(
"%s has no CRS; assuming EPSG:%d because all coordinates are within the lon/lat range "
"(bounds: %.6f, %.6f, %.6f, %.6f).",
source,
WGS84_EPSG,
min_x,
min_y,
max_x,
max_y,
)
return gdf.set_crs(epsg=WGS84_EPSG)