"""
Helper functions for computational geometry in PySAL.
"""
__author__ = "Sergio J. Rey, Xinyue Ye, Charles Schmidt, Andrew Winslow"
__credits__ = "Copyright (c) 2005-2009 Sergio J. Rey"
# ruff: noqa: F403, F405, N803
import copy
import math
import random
from itertools import islice
import numpy as np
import scipy.spatial
from .shapes import *
EPSILON_SCALER = 3
__all__ = [
"bbcommon",
"get_bounding_box",
"get_angle_between",
"is_collinear",
"get_segments_intersect",
"get_segment_point_intersect",
"get_polygon_point_intersect",
"get_rectangle_point_intersect",
"get_ray_segment_intersect",
"get_rectangle_rectangle_intersection",
"get_polygon_point_dist",
"get_points_dist",
"get_segment_point_dist",
"get_point_at_angle_and_dist",
"convex_hull",
"is_clockwise",
"point_touches_rectangle",
"get_shared_segments",
"distance_matrix",
]
[docs]
def bbcommon(bb, bbother):
"""Old Stars method for bounding box overlap testing.
Also defined in ``pysal.weights._cont_binning``.
Parameters
----------
bb : list
A bounding box.
bbother : list
The bounding box to test against.
Returns
-------
chflag : int
``1`` if ``bb`` overlaps ``bbother``, otherwise ``0``.
Examples
--------
>>> b0 = [0, 0, 10, 10]
>>> b1 = [10, 0, 20, 10]
>>> bbcommon(b0, b1)
1
"""
chflag = 0
if (not ((bbother[2] < bb[0]) or (bbother[0] > bb[2]))) and (
not ((bbother[3] < bb[1]) or (bbother[1] > bb[3]))
):
chflag = 1
return chflag
[docs]
def get_bounding_box(items):
"""Find bounding box for a list of geometries.
Parameters
----------
items : list
PySAL shapes.
Returns
-------
rect = libpysal.cg.Rectangle
The bounding box for a list of geometries.
Examples
--------
>>> bb = get_bounding_box([Point((-1, 5)), Rectangle(0, 6, 11, 12)])
>>> bb.left
-1.0
>>> bb.lower
5.0
>>> bb.right
11.0
>>> bb.upper
12.0
"""
def left(o):
# Polygon, Ellipse
if hasattr(o, "bounding_box"):
return o.bounding_box.left
# Rectangle
elif hasattr(o, "left"):
return o.left
# Point
else:
return o[0]
def right(o):
# Polygon, Ellipse
if hasattr(o, "bounding_box"):
return o.bounding_box.right
# Rectangle
elif hasattr(o, "right"):
return o.right
# Point
else:
return o[0]
def lower(o):
# Polygon, Ellipse
if hasattr(o, "bounding_box"):
return o.bounding_box.lower
# Rectangle
elif hasattr(o, "lower"):
return o.lower
# Point
else:
return o[1]
def upper(o):
# Polygon, Ellipse
if hasattr(o, "bounding_box"):
return o.bounding_box.upper
# Rectangle
elif hasattr(o, "upper"):
return o.upper
# Point
else:
return o[1]
rect = Rectangle(
min(list(map(left, items))),
min(list(map(lower, items))),
max(list(map(right, items))),
max(list(map(upper, items))),
)
return rect
[docs]
def get_angle_between(ray1, ray2):
"""Returns the angle formed between a pair of rays which share an origin.
Parameters
----------
ray1 : libpysal.cg.Ray
A ray forming the beginning of the angle measured.
ray2 : libpysal.cg.Ray
A ray forming the end of the angle measured.
Returns
-------
angle : float
The angle between ``ray1`` and ``ray2``.
Raises
------
ValueError
Raised when rays do not have the same origin.
Examples
--------
>>> get_angle_between(
... Ray(Point((0, 0)), Point((1, 0))),
... Ray(Point((0, 0)), Point((1, 0)))
... )
0.0
"""
if ray1.o != ray2.o:
raise ValueError("Rays must have the same origin.")
vec1 = (ray1.p[0] - ray1.o[0], ray1.p[1] - ray1.o[1])
vec2 = (ray2.p[0] - ray2.o[0], ray2.p[1] - ray2.o[1])
rot_theta = -math.atan2(vec1[1], vec1[0])
rot_matrix = [
[math.cos(rot_theta), -math.sin(rot_theta)],
[math.sin(rot_theta), math.cos(rot_theta)],
]
rot_vec2 = (
rot_matrix[0][0] * vec2[0] + rot_matrix[0][1] * vec2[1],
rot_matrix[1][0] * vec2[0] + rot_matrix[1][1] * vec2[1],
)
angle = math.atan2(rot_vec2[1], rot_vec2[0])
return angle
[docs]
def is_collinear(p1, p2, p3):
"""Returns whether a triplet of points is collinear.
Parameters
----------
p1 : libpysal.cg.Point
A point.
p2 : libpysal.cg.Point
A point.
p3 : libpysal.cg.Point
A point.
Returns
-------
collinear : bool
``True`` if ``{p1, p2, p3}`` are collinear, otherwise ``False``.
Examples
--------
>>> is_collinear(Point((0, 0)), Point((1, 1)), Point((5, 5)))
True
>>> is_collinear(Point((0, 0)), Point((1, 1)), Point((5, 0)))
False
"""
eps = np.finfo(type(p1[0])).eps
slope_diff = abs(
(p2[0] - p1[0]) * (p3[1] - p1[1]) - (p2[1] - p1[1]) * (p3[0] - p1[0])
)
very_small_dist = EPSILON_SCALER * eps
collinear = slope_diff < very_small_dist
return collinear
[docs]
def get_segments_intersect(seg1, seg2):
"""Returns the intersection of two segments if one exists.
Parameters
----------
seg1 : libpysal.cg.LineSegment
A segment to check for an intersection.
seg2 : libpysal.cg.LineSegment
The segment to check against ``seg1`` for an intersection.
Returns
-------
intersection : {libpysal.cg.Point, libpysal.cg.LineSegment, None}
The intersecting point or line between ``seg1`` and
``seg2`` if an intersection exists or ``None`` if
``seg1`` and ``seg2`` do not intersect.
Examples
--------
>>> seg1 = LineSegment(Point((0, 0)), Point((0, 10)))
>>> seg2 = LineSegment(Point((-5, 5)), Point((5, 5)))
>>> i = get_segments_intersect(seg1, seg2)
>>> isinstance(i, Point)
True
>>> str(i)
'(0.0, 5.0)'
>>> seg3 = LineSegment(Point((100, 100)), Point((100, 101)))
>>> i = get_segments_intersect(seg2, seg3)
"""
p1 = seg1.p1
p2 = seg1.p2
p3 = seg2.p1
p4 = seg2.p2
a = p2[0] - p1[0]
b = p3[0] - p4[0]
c = p2[1] - p1[1]
d = p3[1] - p4[1]
det = float(a * d - b * c)
intersection = None
if det == 0:
if seg1 == seg2:
intersection = LineSegment(seg1.p1, seg1.p2)
else:
a = get_segment_point_intersect(seg2, seg1.p1)
b = get_segment_point_intersect(seg2, seg1.p2)
c = get_segment_point_intersect(seg1, seg2.p1)
d = get_segment_point_intersect(seg1, seg2.p2)
if a and b: # seg1 in seg2
intersection = LineSegment(seg1.p1, seg1.p2)
if c and d: # seg2 in seg1
intersection = LineSegment(seg2.p1, seg2.p2)
if (a or b) and (c or d):
p1 = a if a else b
p2 = c if c else d
intersection = LineSegment(p1, p2)
else:
a_inv = d / det
b_inv = -b / det
c_inv = -c / det
d_inv = a / det
m = p3[0] - p1[0]
n = p3[1] - p1[1]
x = a_inv * m + b_inv * n
y = c_inv * m + d_inv * n
intersect_exists = 0 <= x <= 1 and 0 <= y <= 1
if intersect_exists:
intersection = Point(
(p1[0] + x * (p2[0] - p1[0]), p1[1] + x * (p2[1] - p1[1]))
)
return intersection
[docs]
def get_segment_point_intersect(seg, pt):
"""Returns the intersection of a segment and point.
Parameters
----------
seg : libpysal.cg.LineSegment
A segment to check for an intersection.
pt : libpysal.cg.Point
A point to check ``seg`` for an intersection.
Returns
-------
pt : {libpysal.cg.Point, None}
The intersection of a ``seg`` and ``pt`` if one exists, otherwise ``None``.
Examples
--------
>>> seg = LineSegment(Point((0, 0)), Point((0, 10)))
>>> pt = Point((0, 5))
>>> i = get_segment_point_intersect(seg, pt)
>>> str(i)
'(0.0, 5.0)'
>>> pt2 = Point((5, 5))
>>> get_segment_point_intersect(seg, pt2)
"""
eps = np.finfo(type(pt[0])).eps
if is_collinear(pt, seg.p1, seg.p2):
if get_segment_point_dist(seg, pt)[0] < EPSILON_SCALER * eps:
pass
else:
pt = None
else:
vec1 = (pt[0] - seg.p1[0], pt[1] - seg.p1[1])
vec2 = (seg.p2[0] - seg.p1[0], seg.p2[1] - seg.p1[1])
if abs(vec1[0] * vec2[1] - vec1[1] * vec2[0]) < eps:
pass
else:
pt = None
return pt
[docs]
def get_polygon_point_intersect(poly, pt):
"""Returns the intersection of a polygon and point.
Parameters
----------
poly : libpysal.cg.Polygon
A polygon to check for an intersection.
pt : libpysal.cg.Point
A point to check ``poly`` for an intersection.
Returns
-------
ret : {libpysal.cg.Point, None}
The intersection of a ``poly`` and ``pt`` if one exists, otherwise ``None``.
Examples
--------
>>> poly = Polygon([Point((0, 0)), Point((1, 0)), Point((1, 1)), Point((0, 1))])
>>> pt = Point((0.5, 0.5))
>>> i = get_polygon_point_intersect(poly, pt)
>>> str(i)
'(0.5, 0.5)'
>>> pt2 = Point((2, 2))
>>> get_polygon_point_intersect(poly, pt2)
"""
def pt_lies_on_part_boundary(p, vx):
vx_range = range(-1, len(vx) - 1)
seg = lambda i: LineSegment(vx[i], vx[i + 1]) # noqa: E731
return [i for i in vx_range if get_segment_point_dist(seg(i), p)[0] == 0] != []
ret = None
# Weed out points that aren't even close
if get_rectangle_point_intersect(poly.bounding_box, pt) is None:
pass
else:
if [
vxs for vxs in poly._vertices if pt_lies_on_part_boundary(pt, vxs)
] != [] or [vxs for vxs in poly._vertices if _point_in_vertices(pt, vxs)] != []:
ret = pt
if poly._holes != [[]]:
if [vxs for vxs in poly.holes if pt_lies_on_part_boundary(pt, vxs)] != []:
# pt lies on boundary of hole.
pass
if [vxs for vxs in poly.holes if _point_in_vertices(pt, vxs)] != []:
# pt lines inside a hole.
ret = None
# raise NotImplementedError,
# 'Cannot compute containment for polygon with holes'
return ret
[docs]
def get_rectangle_point_intersect(rect, pt):
"""Returns the intersection of a rectangle and point.
Parameters
----------
rect : libpysal.cg.Rectangle
A rectangle to check for an intersection.
pt : libpysal.cg.Point
A point to check ``rect`` for an intersection.
Returns
-------
pt : {libpysal.cg.Point, None}
The intersection of a ``rect`` and ``pt`` if one exists, otherwise ``None``.
Examples
--------
>>> rect = Rectangle(0, 0, 5, 5)
>>> pt = Point((1, 1))
>>> i = get_rectangle_point_intersect(rect, pt)
>>> str(i)
'(1.0, 1.0)'
>>> pt2 = Point((10, 10))
>>> get_rectangle_point_intersect(rect, pt2)
"""
if rect.left <= pt[0] <= rect.right and rect.lower <= pt[1] <= rect.upper:
pass
else:
pt = None
return pt
[docs]
def get_ray_segment_intersect(ray, seg):
"""Returns the intersection of a ray and line segment.
Parameters
----------
ray : libpysal.cg.Ray
A ray to check for an intersection.
seg : libpysal.cg.LineSegment
A segment to check for an intersection against ``ray``.
Returns
-------
intersection : {libpysal.cg.Point, libpysal.cg.LineSegment, None}
The intersecting point or line between ``ray`` and
``seg`` if an intersection exists or ``None`` if
``ray`` and ``seg`` do not intersect.
See Also
--------
libpysal.cg.get_segments_intersect
Examples
--------
>>> ray = Ray(Point((0, 0)), Point((0, 1)))
>>> seg = LineSegment(Point((-1, 10)), Point((1, 10)))
>>> i = get_ray_segment_intersect(ray, seg)
>>> isinstance(i, Point)
True
>>> str(i)
'(0.0, 10.0)'
>>> seg2 = LineSegment(Point((10, 10)), Point((10, 11)))
>>> get_ray_segment_intersect(ray, seg2)
"""
# Upper bound on origin to segment dist (+1)
d = (
max(
math.hypot(seg.p1[0] - ray.o[0], seg.p1[1] - ray.o[1]),
math.hypot(seg.p2[0] - ray.o[0], seg.p2[1] - ray.o[1]),
)
+ 1
)
ratio = d / math.hypot(ray.o[0] - ray.p[0], ray.o[1] - ray.p[1])
ray_seg = LineSegment(
ray.o,
Point(
(
ray.o[0] + ratio * (ray.p[0] - ray.o[0]),
ray.o[1] + ratio * (ray.p[1] - ray.o[1]),
)
),
)
intersection = get_segments_intersect(seg, ray_seg)
return intersection
[docs]
def get_rectangle_rectangle_intersection(r0, r1, checkOverlap=True):
"""Returns the intersection between two rectangles.
Parameters
----------
r0 : libpysal.cg.Rectangle
A rectangle to check for an intersection.
r1 : libpysal.cg.Rectangle
A rectangle to check for an intersection against ``r0``.
checkOverlap : bool
Call ``bbcommon(r0, r1)`` prior to complex geometry
checking. Default is ``True``. Prior to setting as
``False`` see the Notes section.
Returns
-------
intersection : {libpysal.cg.Point, libpysal.cg.LineSegment, libpysal.cg.Rectangle, None}
The intersecting point, line, or rectangle between
``r0`` and ``r1`` if an intersection exists or ``None``
if ``r0`` and ``r1`` do not intersect.
Notes
-----
The algorithm assumes the rectangles overlap. The keyword
``checkOverlap=False`` should be used with extreme caution.
Examples
--------
>>> r0 = Rectangle(0,4,6,9)
>>> r1 = Rectangle(4,0,9,7)
>>> ri = get_rectangle_rectangle_intersection(r0,r1)
>>> ri[:]
[4.0, 4.0, 6.0, 7.0]
>>> r0 = Rectangle(0,0,4,4)
>>> r1 = Rectangle(2,1,6,3)
>>> ri = get_rectangle_rectangle_intersection(r0,r1)
>>> ri[:]
[2.0, 1.0, 4.0, 3.0]
>>> r0 = Rectangle(0,0,4,4)
>>> r1 = Rectangle(2,1,3,2)
>>> ri = get_rectangle_rectangle_intersection(r0,r1)
>>> ri[:] == r1[:]
True
""" # noqa: E501
intersection = None
common_bb = True
if checkOverlap and not bbcommon(r0, r1):
# raise ValueError, "Rectangles do not intersect"
common_bb = False
if common_bb:
left = max(r0.left, r1.left)
lower = max(r0.lower, r1.lower)
right = min(r0.right, r1.right)
upper = min(r0.upper, r1.upper)
if upper == lower and left == right:
intersection = Point((left, lower))
elif upper == lower:
intersection = LineSegment(Point((left, lower)), Point((right, lower)))
elif left == right:
intersection = LineSegment(Point((left, lower)), Point((left, upper)))
else:
intersection = Rectangle(left, lower, right, upper)
return intersection
[docs]
def get_polygon_point_dist(poly, pt):
"""Returns the distance between a polygon and point.
Parameters
----------
poly : libpysal.cg.Polygon
A polygon to compute distance from.
pt : libpysal.cg.Point
a point to compute distance from
Returns
-------
dist : float
The distance between ``poly`` and ``point``.
Examples
--------
>>> poly = Polygon([Point((0, 0)), Point((1, 0)), Point((1, 1)), Point((0, 1))])
>>> pt = Point((2, 0.5))
>>> get_polygon_point_dist(poly, pt)
1.0
>>> pt2 = Point((0.5, 0.5))
>>> get_polygon_point_dist(poly, pt2)
0.0
"""
if get_polygon_point_intersect(poly, pt) is not None:
dist = 0.0
else:
part_prox = []
for vertices in poly._vertices:
vx_range = range(-1, len(vertices) - 1)
seg = lambda i: LineSegment( # noqa: E731
vertices[i], # noqa: B023
vertices[i + 1], # noqa: B023
)
_min_dist = min([get_segment_point_dist(seg(i), pt)[0] for i in vx_range])
part_prox.append(_min_dist)
dist = min(part_prox)
return dist
[docs]
def get_points_dist(pt1, pt2):
"""Returns the distance between a pair of points.
Parameters
----------
pt1 : libpysal.cg.Point
A point.
pt2 : libpysal.cg.Point
The other point.
Returns
-------
dist : float
The distance between ``pt1`` and ``pt2``.
Examples
--------
>>> get_points_dist(Point((4, 4)), Point((4, 8)))
4.0
>>> get_points_dist(Point((0, 0)), Point((0, 0)))
0.0
"""
dist = math.hypot(pt1[0] - pt2[0], pt1[1] - pt2[1])
return dist
[docs]
def get_segment_point_dist(seg, pt):
"""Returns (1) the distance between a line segment and point
and (2) the distance along the segment to the closest location on the
segment from the point as a ratio of the length of the segment.
Parameters
----------
seg : libpysal.cg.LineSegment
A line segment to compute distance from.
pt : libpysal.cg.Point
A point to compute distance from.
Returns
-------
dist : float
The distance between ``seg`` and ``pt``.
ratio : float
The distance along ``seg`` to the closest location on
``seg`` from ``pt`` as a ratio of the length of ``seg``.
Examples
--------
>>> seg = LineSegment(Point((0, 0)), Point((10, 0)))
>>> pt = Point((5, 5))
>>> get_segment_point_dist(seg, pt)
(5.0, 0.5)
>>> pt2 = Point((0, 0))
>>> get_segment_point_dist(seg, pt2)
(0.0, 0.0)
"""
src_p = seg.p1
dest_p = seg.p2
# Shift line to go through origin
points_0 = pt[0] - src_p[0]
points_1 = pt[1] - src_p[1]
points_2 = 0
points_3 = 0
points_4 = dest_p[0] - src_p[0]
points_5 = dest_p[1] - src_p[1]
segment_length = get_points_dist(src_p, dest_p)
# Meh, robustness...
# maybe should incorporate this into a more general approach later
if segment_length == 0:
dist, ratio = get_points_dist(pt, src_p), 0
else:
u_x = points_4 / segment_length
u_y = points_5 / segment_length
inter_x = u_x * u_x * points_0 + u_x * u_y * points_1
inter_y = u_x * u_y * points_0 + u_y * u_y * points_1
src_proj_dist = get_points_dist((0, 0), (inter_x, inter_y))
dest_proj_dist = get_points_dist((inter_x, inter_y), (points_4, points_5))
if src_proj_dist > segment_length or dest_proj_dist > segment_length:
src_pt_dist = get_points_dist((points_2, points_3), (points_0, points_1))
dest_pt_dist = get_points_dist((points_4, points_5), (points_0, points_1))
if src_pt_dist < dest_pt_dist:
dist, ratio = src_pt_dist, 0
else:
dist, ratio = dest_pt_dist, 1
else:
dist = get_points_dist((inter_x, inter_y), (points_0, points_1))
ratio = src_proj_dist / segment_length
return dist, ratio
[docs]
def get_point_at_angle_and_dist(ray, angle, dist):
"""Returns the point at a distance and angle relative to the origin of a ray.
Parameters
----------
ray : libpysal.cg.Ray
The ray to which ``angle`` and ``dist`` are relative.
angle : float
The angle relative to ``ray`` at which ``point`` is located.
dist : float
The distance from the origin of ``ray`` at which ``point`` is located.
Returns
-------
point : libpysal.cg.Point
The point at ``dist`` and ``angle`` relative to the origin of ``ray``.
Examples
--------
>>> ray = Ray(Point((0, 0)), Point((1, 0)))
>>> pt = get_point_at_angle_and_dist(ray, math.pi, 1.0)
>>> isinstance(pt, Point)
True
>>> round(pt[0], 8)
-1.0
>>> round(pt[1], 8)
0.0
"""
v = (ray.p[0] - ray.o[0], ray.p[1] - ray.o[1])
cur_angle = math.atan2(v[1], v[0])
dest_angle = cur_angle + angle
point = Point(
(ray.o[0] + dist * math.cos(dest_angle), ray.o[1] + dist * math.sin(dest_angle))
)
return point
[docs]
def convex_hull(points):
"""Returns the convex hull of a set of points.
Parameters
----------
points : list
A list of points for computing the convex hull.
Returns
-------
stack : list
A list of points representing the convex hull.
Examples
--------
>>> points = [Point((0, 0)), Point((4, 4)), Point((4, 0)), Point((3, 1))]
>>> convex_hull(points)
[(0.0, 0.0), (4.0, 0.0), (4.0, 4.0)]
"""
def right_turn(p1, p2, p3) -> bool:
"""Returns if ``p1`` -> ``p2`` -> ``p3`` forms a 'right turn'."""
vec1 = (p2[0] - p1[0], p2[1] - p1[1])
vec2 = (p3[0] - p2[0], p3[1] - p2[1])
_rt = vec2[0] * vec1[1] - vec2[1] * vec1[0] >= 0
return _rt
points = copy.copy(points)
lowest = min(points, key=lambda p: (p[1], p[0]))
points.remove(lowest)
points.sort(key=lambda p: math.atan2(p[1] - lowest[1], p[0] - lowest[0]))
stack = [lowest]
for p in points:
stack.append(p)
while len(stack) > 3 and right_turn(stack[-3], stack[-2], stack[-1]):
stack.pop(-2)
return stack
[docs]
def is_clockwise(vertices):
"""Returns whether a list of points describing
a polygon are clockwise or counterclockwise.
Parameters
----------
vertices : list
A list of points that form a single ring.
Returns
-------
clockwise : bool
``True`` if ``vertices`` are clockwise, otherwise ``False``.
See Also
--------
libpysal.cg.ccw
Examples
--------
>>> is_clockwise([Point((0, 0)), Point((10, 0)), Point((0, 10))])
False
>>> is_clockwise([Point((0, 0)), Point((0, 10)), Point((10, 0))])
True
>>> v = [
... (-106.57798, 35.174143999999998),
... (-106.583412, 35.174141999999996),
... (-106.58417999999999, 35.174143000000001),
... (-106.58377999999999, 35.175542999999998),
... (-106.58287999999999, 35.180543),
... (-106.58263099999999, 35.181455),
... (-106.58257999999999, 35.181643000000001),
... (-106.58198299999999, 35.184615000000001),
... (-106.58148, 35.187242999999995),
... (-106.58127999999999, 35.188243),
... (-106.58138, 35.188243),
... (-106.58108, 35.189442999999997),
... (-106.58104, 35.189644000000001),
... (-106.58028, 35.193442999999995),
... (-106.580029, 35.194541000000001),
... (-106.57974399999999, 35.195785999999998),
... (-106.579475, 35.196961999999999),
... (-106.57922699999999, 35.198042999999998),
... (-106.578397, 35.201665999999996),
... (-106.57827999999999, 35.201642999999997),
... (-106.57737999999999, 35.201642999999997),
... (-106.57697999999999, 35.201543000000001),
... (-106.56436599999999, 35.200311999999997),
... (-106.56058, 35.199942999999998),
... (-106.56048, 35.197342999999996),
... (-106.56048, 35.195842999999996),
... (-106.56048, 35.194342999999996),
... (-106.56048, 35.193142999999999),
... (-106.56048, 35.191873999999999),
... (-106.56048, 35.191742999999995),
... (-106.56048, 35.190242999999995),
... (-106.56037999999999, 35.188642999999999),
... (-106.56037999999999, 35.187242999999995),
... (-106.56037999999999, 35.186842999999996),
... (-106.56037999999999, 35.186552999999996),
... (-106.56037999999999, 35.185842999999998),
... (-106.56037999999999, 35.184443000000002),
... (-106.56037999999999, 35.182943000000002),
... (-106.56037999999999, 35.181342999999998),
... (-106.56037999999999, 35.180433000000001),
... (-106.56037999999999, 35.179943000000002),
... (-106.56037999999999, 35.178542999999998),
... (-106.56037999999999, 35.177790999999999),
... (-106.56037999999999, 35.177143999999998),
... (-106.56037999999999, 35.175643999999998),
... (-106.56037999999999, 35.174444000000001),
... (-106.56037999999999, 35.174043999999995),
... (-106.560526, 35.174043999999995),
... (-106.56478, 35.174043999999995),
... (-106.56627999999999, 35.174143999999998),
... (-106.566541, 35.174144999999996),
... (-106.569023, 35.174157000000001),
... (-106.56917199999999, 35.174157999999998),
... (-106.56938, 35.174143999999998),
... (-106.57061499999999, 35.174143999999998),
... (-106.57097999999999, 35.174143999999998),
... (-106.57679999999999, 35.174143999999998),
... (-106.57798, 35.174143999999998)
... ]
>>> is_clockwise(v)
True
"""
clockwise = True
if not len(vertices) < 3:
area = 0.0
ax, ay = vertices[0]
for bx, by in vertices[1:]:
area += ax * by - ay * bx
ax, ay = bx, by
bx, by = vertices[0]
area += ax * by - ay * bx
clockwise = area < 0.0
return clockwise
def ccw(vertices):
"""Returns whether a list of points is counterclockwise.
Parameters
----------
vertices : list
A list of points that form a single ring.
Returns
-------
counter_clockwise : bool
``True`` if ``vertices`` are counter clockwise, otherwise ``False``.
See Also
--------
libpysal.cg.is_clockwise
Examples
--------
>>> ccw([Point((0, 0)), Point((10, 0)), Point((0, 10))])
True
>>> ccw([Point((0, 0)), Point((0, 10)), Point((10, 0))])
False
"""
counter_clockwise = True
if is_clockwise(vertices):
counter_clockwise = False
return counter_clockwise
def seg_intersect(a, b, c, d):
"""Tests if two segments (a,b) and (c,d) intersect.
Parameters
----------
a : libpysal.cg.Point
The first vertex for the first segment.
b : libpysal.cg.Point
The second vertex for the first segment.
c : libpysal.cg.Point
The first vertex for the second segment.
d : libpysal.cg.Point
The second vertex for the second segment.
Returns
-------
segments_intersect : bool
``True`` if segments ``(a,b)`` and ``(c,d)``, otherwise ``False``.
Examples
--------
>>> a = Point((0,1))
>>> b = Point((0,10))
>>> c = Point((-2,5))
>>> d = Point((2,5))
>>> e = Point((-3,5))
>>> seg_intersect(a, b, c, d)
True
>>> seg_intersect(a, b, c, e)
False
"""
segments_intersect = True
acd_bcd = ccw([a, c, d]) == ccw([b, c, d])
abc_abd = ccw([a, b, c]) == ccw([a, b, d])
if acd_bcd or abc_abd:
segments_intersect = False
return segments_intersect
def _point_in_vertices(pt, vertices):
"""**HELPER METHOD. DO NOT CALL.** Returns whether a point
is contained in a polygon specified by a sequence of vertices.
Parameters
----------
pt : libpysal.cg.Point
A point.
vertices : list
A list of vertices representing as polygon.
Returns
-------
pt_in_poly : bool
``True`` if ``pt`` is contained in ``vertices``, otherwise ``False``.
Examples
--------
>>> _point_in_vertices(
... Point((1, 1)),
... [Point((0, 0)), Point((10, 0)), Point((0, 10))]
... )
True
"""
def neg_ray_intersect(p1, p2, p3) -> bool:
"""Returns whether a ray in the negative-x
direction from ``p3`` intersects the segment between.
"""
if not min(p1[1], p2[1]) <= p3[1] <= max(p1[1], p2[1]):
nr_inters = False
else:
if p1[1] > p2[1]:
vec1 = (p2[0] - p1[0], p2[1] - p1[1])
else:
vec1 = (p1[0] - p2[0], p1[1] - p2[1])
vec2 = (p3[0] - p1[0], p3[1] - p1[1])
nr_inters = vec1[0] * vec2[1] - vec2[0] * vec1[1] >= 0
return nr_inters
vert_y_set = {v[1] for v in vertices}
while pt[1] in vert_y_set:
# Perturb the location very slightly
pt = pt[0], pt[1] + -1e-14 + random.random() * 2e-14
inters = 0
for i in range(-1, len(vertices) - 1):
v1 = vertices[i]
v2 = vertices[i + 1]
if neg_ray_intersect(v1, v2, pt):
inters += 1
pt_in_poly = inters % 2 == 1
return pt_in_poly
[docs]
def point_touches_rectangle(point, rect):
"""Returns ``True`` (``1``) if the point is in the rectangle
or touches it's boundary, otherwise ``False`` (``0``).
Parameters
----------
point : {libpysal.cg.Point, tuple}
A point or point coordinates.
rect : libpysal.cg.Rectangle
A rectangle.
Returns
-------
chflag : int
``1`` if ``point`` is in (or touches
boundary of) ``rect``, otherwise ``0``.
Examples
--------
>>> rect = Rectangle(0, 0, 10, 10)
>>> a = Point((5, 5))
>>> b = Point((10, 5))
>>> c = Point((11, 11))
>>> point_touches_rectangle(a, rect)
1
>>> point_touches_rectangle(b, rect)
1
>>> point_touches_rectangle(c, rect)
0
"""
chflag = 0
if (point[0] >= rect.left and point[0] <= rect.right) and (
point[1] >= rect.lower and point[1] <= rect.upper
):
chflag = 1
return chflag
[docs]
def get_shared_segments(poly1, poly2, bool_ret=False):
"""Returns the line segments in common to both polygons.
Parameters
----------
poly1 : libpysal.cg.Polygon
A Polygon.
poly2 : libpysal.cg.Polygon
A Polygon.
bool_ret : bool
Return only a ``bool``. Default is ``False``.
Returns
-------
common : list
The shared line segments between ``poly1`` and ``poly2``.
_ret_bool : bool
Whether ``poly1`` and ``poly2`` share a
segment (``True``) or not (``False``).
Examples
--------
>>> from libpysal.cg.shapes import Polygon
>>> x = [0, 0, 1, 1]
>>> y = [0, 1, 1, 0]
>>> poly1 = Polygon(list(map(Point, zip(x, y))) )
>>> x = [a+1 for a in x]
>>> poly2 = Polygon(list(map(Point, zip(x, y))) )
>>> get_shared_segments(poly1, poly2, bool_ret=True)
True
"""
# get_rectangle_rectangle_intersection inlined for speed.
r0 = poly1.bounding_box
r1 = poly2.bounding_box
w_left = max(r0.left, r1.left)
w_lower = max(r0.lower, r1.lower)
w_right = min(r0.right, r1.right)
w_upper = min(r0.upper, r1.upper)
segments_a = set()
common = []
for part in poly1.parts + [p for p in poly1.holes if p]:
if part[0] != part[-1]: # not closed
part = part[:] + part[0:1]
a = part[0]
for b in islice(part, 1, None):
# inlining point_touches_rectangle for speed
x, y = a
# check if point a is in the bounding box intersection
if x >= w_left and x <= w_right and y >= w_lower and y <= w_upper:
x, y = b
# check if point b is in the bounding box intersection
if x >= w_left and x <= w_right and y >= w_lower and y <= w_upper:
if a > b:
segments_a.add((b, a))
else:
segments_a.add((a, b))
a = b
_ret_bool = False
for part in poly2.parts + [p for p in poly2.holes if p]:
if part[0] != part[-1]: # not closed
part = part[:] + part[0:1]
a = part[0]
for b in islice(part, 1, None):
# inlining point_touches_rectangle for speed
x, y = a
if x >= w_left and x <= w_right and y >= w_lower and y <= w_upper:
x, y = b
if x >= w_left and x <= w_right and y >= w_lower and y <= w_upper:
seg = (b, a) if a > b else (a, b)
if seg in segments_a:
common.append(LineSegment(*seg))
if bool_ret:
_ret_bool = True
return _ret_bool
a = b
if bool_ret:
if len(common) > 0:
_ret_bool = True
return _ret_bool
return common
[docs]
def distance_matrix(X, p=2.0, threshold=5e7):
r"""Calculate a distance matrix.
Parameters
----------
X : numpy.ndarray
An :math:`n \\times k` array where :math:`n` is the number
of observations and :math:`k` is the number of dimensions
(2 for :math:`x,y`).
p : float
Minkowski `p`-norm distance metric parameter where
:math:`1<=\mathtt{p}<=\infty`. ``2`` is Euclidean distance and
``1`` is Manhattan distance. Default is ``2.0``.
threshold : int
If :math:`(\mathtt{n}**2)*32 > \mathtt{threshold}` use
``scipy.spatial.distance_matrix`` instead of working in RAM,
this is roughly the amount of RAM (in bytes) that will be used.
Must be positive. Default is ``5e7``.
Returns
-------
d : numpy.ndarray
An n by :math:`m` :math:`p`-norm distance matrix.
Raises
------
TypeError
Raised when an invalid dimensional array is passed in.
Notes
-----
Needs optimization/integration with other weights in PySAL.
Examples
--------
>>> x, y = [r.flatten() for r in np.indices((3, 3))]
>>> data = np.array([x, y]).T
>>> d = distance_matrix(data)
>>> np.array(d)
array([[0. , 1. , 2. , 1. , 1.41421356,
2.23606798, 2. , 2.23606798, 2.82842712],
[1. , 0. , 1. , 1.41421356, 1. ,
1.41421356, 2.23606798, 2. , 2.23606798],
[2. , 1. , 0. , 2.23606798, 1.41421356,
1. , 2.82842712, 2.23606798, 2. ],
[1. , 1.41421356, 2.23606798, 0. , 1. ,
2. , 1. , 1.41421356, 2.23606798],
[1.41421356, 1. , 1.41421356, 1. , 0. ,
1. , 1.41421356, 1. , 1.41421356],
[2.23606798, 1.41421356, 1. , 2. , 1. ,
0. , 2.23606798, 1.41421356, 1. ],
[2. , 2.23606798, 2.82842712, 1. , 1.41421356,
2.23606798, 0. , 1. , 2. ],
[2.23606798, 2. , 2.23606798, 1.41421356, 1. ,
1.41421356, 1. , 0. , 1. ],
[2.82842712, 2.23606798, 2. , 2.23606798, 1.41421356,
1. , 2. , 1. , 0. ]])
"""
if X.ndim == 1:
X.shape = (X.shape[0], 1)
if X.ndim > 2:
msg = f"Should be 2D point coordinates: {X.ndim} dimensions present."
raise TypeError(msg)
n, k = X.shape
if (n**2) * 32 > threshold:
d = scipy.spatial.distance_matrix(X, X, p)
else:
m = np.ones((n, n))
d = np.zeros((n, n))
for col in range(k):
x = X[:, col]
x_m = x * m
dx = x_m - x_m.T
if p % 2 != 0:
dx = np.abs(dx)
dx2 = dx**p
d += dx2
d = d ** (1.0 / p)
return d