"""Shape and geometry utilities for working with fiona and shapely objects.
Note this module contains a lot of other simple functions that are commonly
used by other functions, but are not included in documentation because they are
likely not useful to users.
"""
from __future__ import annotations
import typing
if typing.TYPE_CHECKING:
from typing import Iterable, Tuple, List, Any, Optional, Literal
import numpy.typing as npt
import geopandas as gpd
import xarray.core.types
import matplotlib.axes
from watershed_workflow.crs import CRS
import datetime, cftime
import logging
import numpy as np
import math
import scipy.interpolate
import shapely.geometry
from shapely.geometry.base import BaseGeometry
import shapely.ops
import shapely.prepared
import shapely.affinity
import xarray
import watershed_workflow.warp
_tol = 1.e-7
#
# Geometric utilities
#
[docs]
def computeDistance(p1: Tuple[float, float], p2: Tuple[float, float]) -> float:
"""Compute Euclidean distance between two points.
Parameters
----------
p1 : tuple of float
First point as (x, y) coordinates.
p2 : tuple of float
Second point as (x, y) coordinates.
Returns
-------
float
Euclidean distance between the two points.
"""
return math.sqrt((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)
[docs]
def computeTriangleArea(p1: Tuple[float, float], p2: Tuple[float, float],
p3: Tuple[float, float]) -> float:
"""Compute the area of a triangle given three vertices.
Parameters
----------
p1 : tuple of float
First vertex as (x, y) coordinates.
p2 : tuple of float
Second vertex as (x, y) coordinates.
p3 : tuple of float
Third vertex as (x, y) coordinates.
Returns
-------
float
Area of the triangle.
"""
return 0.5 * (p2[0] * p3[1] - p3[0] * p2[1] - p1[0] * p3[1] + p3[0] * p1[1] + p1[0] * p2[1]
- p2[0] * p1[1])
[docs]
def computeTriangleCentroid(p1: Tuple[float, float], p2: Tuple[float, float],
p3: Tuple[float, float]) -> Tuple[float, float]:
"""Compute the centroid of a triangle given three vertices.
Parameters
----------
p1 : tuple of float
First vertex as (x, y) coordinates.
p2 : tuple of float
Second vertex as (x, y) coordinates.
p3 : tuple of float
Third vertex as (x, y) coordinates.
Returns
-------
tuple of float
Centroid of the triangle as (x, y) coordinates.
"""
return np.array([p1, p2, p3]).mean(axis=0)
[docs]
def isCollinear(p1: Tuple[float, float],
p2: Tuple[float, float],
p3: Tuple[float, float],
tol: float = 1e-6) -> bool:
"""Check if three points are collinear within a given tolerance.
Parameters
----------
p1 : tuple of float
First point as (x, y) coordinates.
p2 : tuple of float
Second point as (x, y) coordinates.
p3 : tuple of float
Third point as (x, y) coordinates.
tol : float, optional
Tolerance for collinearity test. Default is 1e-6.
Returns
-------
bool
True if the three points are collinear within tolerance.
"""
x1, y1 = p2[0] - p1[0], p2[1] - p1[1]
x2, y2 = p3[0] - p1[0], p3[1] - p1[1]
return abs(x1*y2 - x2*y1) < tol
[docs]
def computeArea(vertices: Iterable[Tuple[float]]) -> float:
"""Compute the area of a polygon given its vertices.
Parameters
----------
vertices : iterable of tuple
Vertices of the polygon as coordinate tuples.
Returns
-------
float
Area of the polygon.
"""
area = shapely.geometry.Polygon(vertices).area
return area
[docs]
def computeAngle(v1: Tuple[float, float] | shapely.geometry.LineString,
v2: Tuple[float, float] | shapely.geometry.LineString) -> float:
"""Compute the angle (in degrees) of v2 relative to v1 in clockwise direction.
Parameters
----------
v1 : tuple of float or shapely.geometry.LineString
First vector as (x, y) coordinates or downstream-oriented LineString.
v2 : tuple of float or shapely.geometry.LineString
Second vector as (x, y) coordinates or downstream-oriented LineString.
Returns
-------
float
Angle in degrees of v2 relative to v1 in clockwise direction (0-360).
"""
if isinstance(v1, shapely.geometry.LineString):
if isinstance(v2, shapely.geometry.LineString):
assert isClose(v1.coords[0], v2.coords[-1])
c1 = np.array(v1.coords[0:2])[:, 0:2]
return computeAngle(c1[1] - c1[0], v2)
if isinstance(v2, shapely.geometry.LineString):
c2 = np.array(v2.coords[-2:])[:, 0:2]
return computeAngle(v1, c2[0] - c2[1])
x1, y1 = v1
x2, y2 = v2
# Compute the angle of each vector with respect to the positive x-axis
angle1 = math.atan2(y1, x1) # Angle of vec1
angle2 = math.atan2(y2, x2) # Angle of vec2
# Compute the difference in angles, clockwise
delta_angle = angle1 - angle2
# Convert to clockwise angle in degrees
clockwise_angle = math.degrees(delta_angle)
# Normalize the angle to be between 0 and 360 degrees
if clockwise_angle < 0:
clockwise_angle += 360
return clockwise_angle
[docs]
def projectVectorAtAngle(v1: Tuple[float, float] | shapely.geometry.LineString, angle: float,
distance: float) -> Tuple[float, float]:
"""Project a vector from v1 at a specified angle and distance.
Parameters
----------
v1 : tuple of float or shapely.geometry.LineString
Reference vector as (x, y) coordinates or downstream-oriented LineString.
angle : float
Angle in degrees for the projection.
distance : float
Distance (magnitude) of the projected vector.
Returns
-------
tuple of float
Projected vector as (x, y) coordinates with specified angle and distance.
"""
if isinstance(v1, shapely.geometry.LineString):
c1 = np.array(v1.coords[0:2])
return projectVectorAtAngle(c1[1] - c1[0], angle, distance)
x1, y1 = v1
# Compute the angle of vec1 with respect to the positive x-axis
angle1 = math.atan2(y1, x1)
# Convert the clockwise angle to radians
angle_offset = math.radians(angle)
# Compute the resulting angle of vec2
angle2 = angle1 - angle_offset
# Compute the components of vec2 using the distance and angle2
x2 = distance * math.cos(angle2)
y2 = distance * math.sin(angle2)
return (x2, y2)
[docs]
def computeMidpoint(p1: Tuple[float, float] | np.ndarray,
p2: Tuple[float, float] | np.ndarray) -> Tuple[float, float]:
"""Compute the midpoint between two points.
Parameters
----------
p1 : tuple of float
First point as (x, y) coordinates.
p2 : tuple of float
Second point as (x, y) coordinates.
Returns
-------
tuple of float
Midpoint as (x, y) coordinates.
"""
return ((p1[0] + p2[0]) / 2., (p1[1] + p2[1]) / 2.)
[docs]
def findClosestPointIndex(point: Tuple[float, float], points: npt.ArrayLike) -> int:
"""Find the index of the closest point in an array of points.
Parameters
----------
point : tuple of float
Reference point as (x, y) coordinates.
points : array_like
Array of points to search, shape (n, 2).
Returns
-------
int
Index of the closest point in the points array.
"""
points2 = np.asarray(points)
point2 = np.array(point)
dist = np.sum((points2 - point2)**2, axis=1)
return int(np.argmin(dist))
[docs]
def cluster(points: np.ndarray, tol: float) -> Tuple[np.ndarray, np.ndarray]:
"""Cluster points based on distance tolerance.
Parameters
----------
points : np.ndarray
Array of points to cluster, shape (n, 2).
tol : float
Distance tolerance for clustering.
Returns
-------
tuple of np.ndarray
Tuple containing (cluster_indices, cluster_centroids) where cluster_indices
is an array of cluster assignments and cluster_centroids are the centroid
coordinates of each cluster.
"""
import scipy.cluster.hierarchy as hcluster
if type(points) is list:
points = np.array(points)
if len(points) > 1:
indices = hcluster.fclusterdata(points, tol, criterion='distance')
centroids = np.array(
[points[indices == (i + 1)].mean(axis=0) for i in range(indices.max())])
else:
indices = np.array([1, ] * len(points))
centroids = points
return indices - 1, centroids
#
# Shape utilities
#
[docs]
def flatten(list_of_shps: Any) -> List[BaseGeometry]:
"""Flatten a list of shapes by expanding Multi-objects into individual geometries.
Parameters
----------
list_of_shps : iterable
List of shapely geometry objects, may contain Multi-objects.
Returns
-------
list of BaseGeometry
Flattened list containing only single geometry objects.
"""
new_list = []
for shp in list_of_shps:
if isinstance(shp, shapely.geometry.MultiLineString) or \
isinstance(shp, shapely.geometry.MultiPoint) or \
isinstance(shp, shapely.geometry.MultiPolygon):
new_list.extend(list(shp.geoms))
else:
new_list.append(shp)
return new_list
[docs]
def recenter(
objects: Iterable[BaseGeometry],
centering: bool | str | shapely.geometry.Point = True
) -> Tuple[List[BaseGeometry], shapely.geometry.Point]:
"""Center a collection of objects by translating to remove their centroid.
Parameters
----------
objects : iterable of BaseGeometry
Collection of shapely geometry objects to center.
centering : bool, str, or shapely.geometry.Point, optional
Centering method: True or 'geometric' for geometric center,
'mass' for center of mass, or Point for custom center.
Returns
-------
tuple
Tuple containing (centered_objects, centroid_point) where centered_objects
is the list of translated geometries and centroid_point is the original centroid.
Raises
------
ValueError
If centering method is not recognized.
"""
if type(centering) is shapely.geometry.Point:
centroid = centering
elif centering is True or centering == 'geometric':
union = shapely.ops.unary_union(objects)
centroid = shapely.geometry.Point([(union.bounds[0] + union.bounds[2]) / 2.,
(union.bounds[1] + union.bounds[3]) / 2.])
elif centering == 'mass':
union = shapely.ops.unary_union(objects)
centroid = union.centroid
else:
raise ValueError('Centering: option centering = "{}" unknown'.format(centering))
new_objs = [
shapely.affinity.translate(obj, -centroid.coords[0][0], -centroid.coords[0][1])
for obj in objects
]
for new, old in zip(new_objs, objects):
if hasattr(old, 'properties'):
new.properties = old.properties
return new_objs, centroid
[docs]
def intersects(shp1: BaseGeometry, shp2: BaseGeometry) -> bool:
"""Check whether two geometries intersect.
Parameters
----------
shp1 : BaseGeometry
First shapely geometry object.
shp2 : BaseGeometry
Second shapely geometry object.
Returns
-------
bool
True if geometries intersect, False otherwise.
Notes
-----
This function computes the actual intersection rather than using
shapely.intersects() for more reliable results.
"""
inter = shp1.intersection(shp2)
return not isEmpty(inter)
[docs]
def isNonPointIntersection(shp1: BaseGeometry, shp2: BaseGeometry) -> bool:
"""Check whether two geometries intersect with more than just a point.
Parameters
----------
shp1 : BaseGeometry
First shapely geometry object.
shp2 : BaseGeometry
Second shapely geometry object.
Returns
-------
bool
True if geometries intersect with line or area, False if no intersection
or point intersection only.
Notes
-----
This function computes the actual intersection to check its dimensionality.
"""
inter = shp1.intersection(shp2)
return not (isEmpty(inter) or \
isinstance(inter, shapely.geometry.Point))
[docs]
def isVolumetricIntersection(shp1: BaseGeometry, shp2: BaseGeometry) -> bool:
"""Check whether two geometries have a volumetric (area) intersection.
Parameters
----------
shp1 : BaseGeometry
First shapely geometry object.
shp2 : BaseGeometry
Second shapely geometry object.
Returns
-------
bool
True if intersection has positive area, False otherwise.
"""
inter = shp1.intersection(shp2)
return inter.area > 0
[docs]
def filterToShape(df: gpd.GeoDataFrame,
shape: BaseGeometry,
shape_crs: CRS,
method: str = 'contains',
tol: Optional[float] = None) -> gpd.GeoDataFrame:
"""Filter GeoDataFrame features based on spatial relationship with shape.
Parameters
----------
df : gpd.GeoDataFrame
GeoDataFrame to filter.
shape : BaseGeometry
Shapely geometry to use as filter.
shape_crs : CRS
Coordinate reference system of the shape.
method : str, optional
Spatial relationship method: 'contains', 'intersects', or
'non_point_intersection'. Default is 'contains'.
tol : float, optional
Tolerance for spatial operations. Uses default if None.
Returns
-------
gpd.GeoDataFrame
Filtered GeoDataFrame containing features that satisfy the spatial relationship.
Raises
------
ValueError
If method is not one of the supported options.
"""
if shape_crs != df.crs:
shape = watershed_workflow.warp.shply(shape, shape_crs, df.crs)
if tol is None: tol = _tol
if method == 'contains':
shape = shapely.prepared.prep(shape.buffer(2 * tol))
op = shape.contains
elif method == 'intersects':
shape = shapely.prepared.prep(shape.buffer(2 * tol))
op = shape.intersects
elif method == 'non_point_intersection':
op = lambda a: isNonPointIntersection(shape, a)
else:
raise ValueError(
"method must be one of 'intersects', 'contains', or 'non_point_intersection'")
return df[df.geometry.apply(op)]
[docs]
def isEmpty(shply: BaseGeometry | None) -> bool:
"""Check if a shapely geometry is None or empty.
Parameters
----------
shply : BaseGeometry or None
Shapely geometry object to check.
Returns
-------
bool
True if geometry is None or empty.
"""
return shply is None or shply.is_empty
[docs]
def isConvex(points: Iterable[Tuple[float, float]]) -> bool:
"""Check if a set of points forms a convex polygon.
Parameters
----------
points : iterable of tuple
Points as coordinate tuples to check for convexity.
Returns
-------
bool
True if points form a convex polygon.
"""
poly = shapely.geometry.Polygon(points)
return math.isclose(poly.area, poly.convex_hull.area, rel_tol=1e-4)
[docs]
def breakLineStringCollinearity(linestring_coords: np.ndarray, tol: float = 1e-5) -> np.ndarray:
"""Remove collinearity from linestring by adding small orthogonal perturbations.
Parameters
----------
linestring_coords : np.ndarray
Array of linestring coordinates, shape (n, 2).
tol : float, optional
Tolerance for collinearity detection and perturbation size. Default is 1e-5.
Returns
-------
np.ndarray
Modified coordinates with collinearity removed.
"""
# traversing along the linestring, checking 3 consecutive points at a time
for i, (p0, p1,
p2) in enumerate(zip(linestring_coords, linestring_coords[1:], linestring_coords[2:])):
# treating collinearity through a small pertubation
if isCollinear(p0, p1, p2, tol=tol):
dp = p2 - p0
ortho = 10 * tol * np.array([-dp[1], dp[0]]) / np.linalg.norm(dp)
linestring_coords[i + 1] = p1 + ortho
return linestring_coords
[docs]
def reverseLineString(ls: shapely.geometry.LineString) -> shapely.geometry.LineString:
"""Reverse the direction of a LineString.
Parameters
----------
ls : shapely.geometry.LineString
LineString to reverse.
Returns
-------
shapely.geometry.LineString
LineString with reversed coordinate order.
"""
return shapely.geometry.LineString(reversed(ls.coords))
[docs]
def isClose(s1: BaseGeometry, s2: BaseGeometry, tol: float = _tol) -> bool:
"""Check if two shapely geometries are topologically equivalent and geometrically close.
Parameters
----------
s1 : BaseGeometry
First shapely geometry object.
s2 : BaseGeometry
Second shapely geometry object.
tol : float, optional
Tolerance for geometric comparison. Default uses module tolerance.
Returns
-------
bool
True if geometries are topologically equivalent and geometrically close.
Notes
-----
This function handles complex cases like polygon rotations and coordinate
ordering that simple coordinate comparison would miss.
Raises
------
NotImplementedError
If geometry type is not supported.
"""
# deal with Multi* or list objects
def is_multi(thing):
if isinstance(thing, shapely.geometry.MultiPoint):
return True
if isinstance(thing, shapely.geometry.MultiLineString):
return True
if isinstance(thing, shapely.geometry.MultiPolygon):
return True
if isinstance(thing, list):
return True
return False
def local_len(thing):
try:
return len(thing)
except AttributeError:
return len(thing.geoms)
def iter(thing):
assert (is_multi(thing))
if isinstance(thing, list):
for t in thing:
yield t
else:
for t in thing.geoms:
yield t
if is_multi(s1):
if local_len(s1) == 1:
return isClose(next(iter(s1)), s2, tol)
if is_multi(s2):
if local_len(s2) == 1:
return isClose(s1, next(iter(s2)), tol)
if is_multi(s1) and is_multi(s2):
if local_len(s1) != local_len(s2):
return False
return all(isClose(i1, i2, tol) for (i1, i2) in zip(s1, s2))
# points get compared as tuples
if isinstance(s1, shapely.geometry.Point):
return isClose(s1.coords[0], s2, tol)
elif isinstance(s1, np.ndarray) and len(s1.shape) == 1:
return isClose((s1[0], s1[1]), s2)
if isinstance(s2, shapely.geometry.Point):
return isClose(s1, s2.coords[0], tol)
elif isinstance(s2, np.ndarray) and len(s2.shape) == 1:
return isClose(s1, (s2[0], s2[1]))
# types should be the same now
if type(s1) != type(s2):
return False
# compare tuples
if isinstance(s1, tuple):
if len(s1) != len(s2):
return False
return sum((p1 - p2)**2 for p1, p2 in zip(s1, s2)) < tol**2
# compare lines
elif isinstance(s1, shapely.geometry.LineString):
if len(s1.coords) != len(s2.coords):
return False
if np.allclose(np.array(s1.coords), np.array(s2.coords), tol, tol):
return True
if np.allclose(np.array(s1.coords), np.array(list(reversed(s2.coords))), tol, tol):
return True
return False
# compare polygons
elif type(s1) is shapely.geometry.Polygon:
# note, this does not correctly deal with nonequal holes...
if len(s1.boundary.coords) != len(s2.boundary.coords):
return False
ls1 = s1.boundary.coords[:-1]
ls2 = np.array(s2.boundary.coords[:-1])
ls2f = np.flipud(ls2)
return any(np.allclose(ls1, np.roll(ls2, i, 0), tol, tol) for i in range(len(ls2))) or \
any(np.allclose(ls1, np.roll(ls2f, i, 0), tol, tol) for i in range(len(ls2)))
# compare multi-shapes by checking if each one has a match in the
# other and lengths are the same
elif isinstance(s1, (shapely.geometry.MultiPoint, shapely.geometry.MultiLineString,
shapely.geometry.MultiPolygon)):
if len(s1) != len(s2):
return False
good2 = [False, ] * len(s2)
for l1 in s1:
found = False
for i, l2 in enumerate(s2):
if not good2[i]:
if isClose(l1, l2, tol):
good2[i] = True
found = True
break
if not found:
return False
return True
else:
raise NotImplementedError("Not implemented for type '%r'" % type(s1))
[docs]
def contains(s1: BaseGeometry, s2: BaseGeometry, tol: float = _tol) -> bool:
"""Check if one geometry contains another with tolerance for roundoff issues.
Parameters
----------
s1 : BaseGeometry
Container geometry.
s2 : BaseGeometry
Geometry to test for containment.
tol : float, optional
Tolerance buffer for containment test. Default uses module tolerance.
Returns
-------
bool
True if s1 contains s2 within tolerance.
"""
return s1.buffer(tol, 2).contains(s2)
[docs]
class CutError(Exception):
"""Exception raised when cutting geometries fails.
Parameters
----------
message : str
Error message.
line : shapely.geometry.LineString
The line being cut.
seg : shapely.geometry.LineString
The segment causing issues.
cutline : shapely.geometry.LineString
The cutting line.
"""
def __init__(self, message: str, line: shapely.geometry.LineString,
seg: shapely.geometry.LineString, cutline: shapely.geometry.LineString) -> None:
super(Exception, self).__init__(message)
self.line = line
self.seg = seg
self.cutline = cutline
[docs]
def cut(
line1: shapely.geometry.LineString, line2: shapely.geometry.LineString
) -> Tuple[List[shapely.geometry.LineString], List[shapely.geometry.LineString]]:
"""Cut two LineStrings at their intersection point.
Parameters
----------
line1 : shapely.geometry.LineString
First LineString to cut.
line2 : shapely.geometry.LineString
Second LineString to cut.
Returns
-------
tuple of list
Tuple containing (line1_segments, line2_segments) where each is a list
of LineString segments created by the cut operation.
"""
return list(shapely.ops.split(line1, line2).geoms), \
list(shapely.ops.split(line2, line1).geoms)
# def cut(line : shapely.geometry.LineString,
# cutline : shapely.geometry.LineString,
# tol : float = 1.e-5) -> List[shapely.geometry.LineString]:
# def plot():
# from matplotlib import pyplot as plt
# plt.plot(cutline.xy[0], cutline.xy[1], 'k-x', linewidth=3)
# plt.plot(line.xy[0], line.xy[1], 'g-+', linewidth=3)
# assert type(line) is shapely.geometry.LineString
# assert type(cutline) is shapely.geometry.LineString
# assert line.intersects(cutline)
# segs = []
# coords = list(line.coords)
# segcoords = [coords[0], ]
# i = 0
# while i < len(coords) - 1:
# seg = shapely.geometry.LineString(coords[i:i + 2])
# #logging.debug("Intersecting seg %d"%i)
# point = seg.intersection(cutline)
# if type(point) is shapely.geometry.LineString and len(point.coords) == 0:
# #logging.debug("Cut seg no intersection")
# segcoords.append(seg.coords[-1])
# i += 1
# elif type(point) is shapely.geometry.Point:
# #logging.debug("Cut intersected at point")
# #logging.debug(" inter point: %r"%list(point.coords[0]))
# #logging.debug(" seg final point: %r"%list(seg.coords[-1]))
# #logging.debug(" close? = %r"%(isClose(point, seg.coords[-1], tol)))
# if isClose(point, seg.coords[-1], tol):
# # intersects at the far point
# segs.append(shapely.geometry.LineString(segcoords + [point, ]))
# #logging.debug(" appended full linestring: %r"%(list(segs[-1].coords)))
# if (i < len(coords) - 2):
# #logging.debug(" (not the end)")
# segcoords = [point, coords[i + 2]]
# else:
# #logging.debug(" (the end)")
# segcoords = [point, ]
# i += 2 # also skip the next seg, which would also
# # intersect at that seg's start point
# elif isClose(point, seg.coords[0], tol):
# # intersects at the near point
# if i != 0:
# assert (len(segcoords) > 1)
# segs.append(shapely.geometry.LineString(segcoords[:-1] + [point, ]))
# segcoords = [point, ]
# else:
# assert (len(segcoords) == 1)
# segcoords[0] = point
# segcoords.append(seg.coords[-1])
# i += 1
# else:
# # intersects in the middle
# segs.append(shapely.geometry.LineString(segcoords + [point, ]))
# #logging.debug(" appended partial linestring: %r"%(list(segs[-1].coords)))
# segcoords = [point, seg.coords[-1]]
# i += 1
# else:
# print("Dual/multiple section: type = {}".format(type(point)))
# print(" point = {}".format(point))
# raise CutError(
# "Dual/multiple intersection in a single seg... ugh! "
# + "Intersection is of type '{}'".format(type(point)), line, seg, cutline)
# if len(segcoords) > 1:
# segs.append(shapely.geometry.LineString(segcoords))
# return segs
[docs]
def inNeighborhood(shp1: BaseGeometry, shp2: BaseGeometry, tol: float = 0.1) -> bool:
"""Check if two geometries can possibly intersect using bounding box test.
Parameters
----------
shp1 : BaseGeometry
First shapely geometry object.
shp2 : BaseGeometry
Second shapely geometry object.
tol : float, optional
Tolerance for bounding box expansion. Default is 0.1.
Returns
-------
bool
True if bounding boxes indicate possible intersection.
"""
minx1, miny1, maxx1, maxy1 = shp1.bounds
minx2, miny2, maxx2, maxy2 = shp2.bounds
if maxx2 < minx1 - tol or \
maxy2 < miny1 - tol or \
minx2 > maxx1 + tol or \
miny2 > maxy1 + tol:
return False
return True
[docs]
def intersectPointToSegment(point: shapely.geometry.Point, seg_start: shapely.geometry.Point,
seg_end: shapely.geometry.Point) -> shapely.geometry.Point:
"""Find the nearest point on a line segment to a given point.
Parameters
----------
point : shapely.geometry.Point
Point to find nearest location for.
seg_start : shapely.geometry.Point
Start point of the line segment.
seg_end : shapely.geometry.Point
End point of the line segment.
Returns
-------
shapely.geometry.Point
Nearest point on the line segment to the input point.
"""
seg_magnitude = seg_end.distance(seg_start)
assert (seg_magnitude > _tol)
u = ((point.x - seg_start.x) * (seg_end.x - seg_start.x) +
(point.y - seg_start.y) * (seg_end.y - seg_start.y)) \
/ (seg_magnitude ** 2)
# closest point does not fall within the line linestring,
# take the shorter distance to an endpoint
if u < 0.:
return seg_start
elif u > 1.:
return seg_end
else:
ix = seg_start.x + u * (seg_end.x - seg_start.x)
iy = seg_start.y + u * (seg_end.y - seg_start.y)
return shapely.geometry.Point([ix, iy])
[docs]
def findNearestPoint(point: shapely.geometry.Point | Tuple[float, float],
line: shapely.geometry.LineString,
tol: Optional[float] = None) -> Tuple[float, float] | None:
"""Find the nearest point on a LineString to a given point.
Parameters
----------
point : shapely.geometry.Point or tuple of float
Point to find nearest location for.
line : shapely.geometry.LineString
LineString to search along.
tol : float, optional
Distance tolerance. If provided, returns None if nearest point
is farther than tolerance.
Returns
-------
tuple of float or None
Nearest point coordinates as (x, y) tuple, or None if outside tolerance.
"""
if tol is None:
if isinstance(point, tuple):
point = shapely.geometry.Point(point)
return shapely.ops.nearest_points(point, line)[1].coords[0]
else:
if inNeighborhood(shapely.geometry.Point(point), line, tol):
logging.debug(" - in neighborhood")
nearest_p = findNearestPoint(point, line)
assert nearest_p is not None
dist = computeDistance(nearest_p, point)
logging.debug(" - nearest p = {0}, dist = {1}, tol = {2}".format(nearest_p, dist, tol))
if dist < tol:
return nearest_p
return None
[docs]
def removeThirdDimension(
geom: shapely.geometry.base.BaseGeometry) -> shapely.geometry.base.BaseGeometry:
"""Remove the third dimension (Z-coordinate) from a shapely geometry.
Parameters
----------
geom : shapely.geometry.base.BaseGeometry
Input geometry that may have Z-coordinates.
Returns
-------
shapely.geometry.base.BaseGeometry
Geometry with only X and Y coordinates.
"""
def _drop_z(*args):
return tuple(filter(None, [args[0], args[1]]))
return shapely.ops.transform(_drop_z, geom)
[docs]
def computeSegmentLengths(ls: shapely.geometry.LineString) -> np.ndarray:
"""Compute the length of each segment in a LineString.
Parameters
----------
ls : shapely.geometry.LineString
LineString to analyze.
Returns
-------
np.ndarray
Array of segment lengths, shape (n-1,) for n coordinates.
"""
coords = np.array(ls.coords)
return np.linalg.norm((coords[1:] - coords[:-1]), axis=1)
[docs]
def computeArclengths(ls: shapely.geometry.LineString) -> np.ndarray:
"""Compute cumulative arc length at each coordinate in a LineString.
Parameters
----------
ls : shapely.geometry.LineString
LineString to analyze.
Returns
-------
np.ndarray
Array of cumulative arc lengths, shape (n,) for n coordinates.
"""
ds = computeSegmentLengths(ls)
return np.cumsum(np.concatenate([[0.0, ], ds]))
#
# Dataset utilities
#
#
# fiona utilities -- probably need to go away?
#
[docs]
def generateRings(obj: Any) -> Iterable[List[Tuple[float, float]]]:
"""Generate coordinate rings from a fiona shape object.
Parameters
----------
obj : dict
Fiona shape dictionary with 'coordinates' or 'geometry' key.
Yields
------
list of tuple
Each ring as a list of coordinate tuples.
"""
def _generateRings(coords):
for e in coords:
if isinstance(e[0], (float, int)):
yield coords
break
else:
for r in _generateRings(e):
yield r
if 'geometry' in obj:
obj = obj['geometry']
for r in _generateRings(obj['coordinates']):
yield r
[docs]
def generateCoords(obj: Any) -> Iterable[Tuple[float, float]]:
"""Generate all coordinates from a fiona geometry object.
Parameters
----------
obj : dict
Fiona shape dictionary with 'coordinates' or 'geometry' key.
Yields
------
tuple of float
Individual coordinate tuples.
"""
if 'geometry' in obj:
obj = obj['geometry']
if obj['type'] == 'Point':
yield obj['coordinates']
else:
for ring in generateRings(obj):
for c in ring:
yield c