Source code for libpysal.cg.standalone

"""
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 = "Should be 2D point coordinates: %s dimensions present." % X.ndim 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