"""
Computational geometry code for PySAL: Python Spatial Analysis Library.
"""
__author__ = "Sergio J. Rey, Xinyue Ye, Charles Schmidt, Andrew Winslow"
__credits__ = "Copyright (c) 2005-2011 Sergio J. Rey"
# ruff: noqa: B028, F403, F405
import copy
import math
import warnings
from .rtree import *
from .shapes import *
from .standalone import *
__all__ = ["Grid", "BruteForcePointLocator", "PointLocator", "PolygonLocator"]
dep_msg = "is deprecated and will be removed in a future version of libpysal"
[docs]
class Grid:
"""
Representation of a binning data structure.
"""
[docs]
def __init__(self, bounds, resolution):
"""
Returns a grid with specified properties.
__init__(Rectangle, number) -> Grid
Parameters
----------
bounds : the area for the grid to encompass
resolution : the diameter of each bin
Examples
--------
TODO: complete this doctest
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
"""
warnings.warn("Grid " + dep_msg, FutureWarning)
if resolution == 0:
raise Exception("Cannot create grid with resolution 0")
self.res = resolution
self.hash = {}
self.x_range = (bounds.left, bounds.right)
self.y_range = (bounds.lower, bounds.upper)
try:
self.i_range = int(
math.ceil((self.x_range[1] - self.x_range[0]) / self.res)
)
self.j_range = int(
math.ceil((self.y_range[1] - self.y_range[0]) / self.res)
)
except Exception as e:
raise Exception(
"Invalid arguments for Grid(): ("
+ str(x_range)
+ ", "
+ str(y_range)
+ ", "
+ str(res)
+ ")"
) from e
[docs]
def in_grid(self, loc):
"""
Returns whether a 2-tuple location _loc_ lies inside the grid bounds.
Test tag: <tc>#is#Grid.in_grid</tc>
"""
return (
self.x_range[0] <= loc[0] <= self.x_range[1]
and self.y_range[0] <= loc[1] <= self.y_range[1]
)
def __grid_loc(self, loc):
i = min(self.i_range, max(int((loc[0] - self.x_range[0]) / self.res), 0))
j = min(self.j_range, max(int((loc[1] - self.y_range[0]) / self.res), 0))
return (i, j)
[docs]
def add(self, item, pt):
"""
Adds an item to the grid at a specified location.
add(x, Point) -> x
Parameters
----------
item : the item to insert into the grid
pt : the location to insert the item at
Examples
--------
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
>>> g.add('A', Point((4.2, 8.7)))
'A'
"""
if not self.in_grid(pt):
raise Exception(
"Attempt to insert item at location outside grid bounds: " + str(pt)
)
grid_loc = self.__grid_loc(pt)
if grid_loc in self.hash:
self.hash[grid_loc].append((pt, item))
else:
self.hash[grid_loc] = [(pt, item)]
return item
[docs]
def remove(self, item, pt):
"""
Removes an item from the grid at a specified location.
remove(x, Point) -> x
Parameters
----------
item : the item to remove from the grid
pt : the location the item was added at
Examples
--------
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
>>> g.add('A', Point((4.2, 8.7)))
'A'
>>> g.remove('A', Point((4.2, 8.7)))
'A'
"""
if not self.in_grid(pt):
raise Exception(
"Attempt to remove item at location outside grid bounds: " + str(pt)
)
grid_loc = self.__grid_loc(pt)
self.hash[grid_loc].remove((pt, item))
if self.hash[grid_loc] == []:
del self.hash[grid_loc]
return item
[docs]
def bounds(self, bounds):
"""
Returns a list of items found in the grid within the bounds specified.
bounds(Rectangle) -> x list
Parameters
----------
item : the item to remove from the grid
pt : the location the item was added at
Examples
--------
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
>>> g.add('A', Point((1.0, 1.0)))
'A'
>>> g.add('B', Point((4.0, 4.0)))
'B'
>>> g.bounds(Rectangle(0, 0, 3, 3))
['A']
>>> g.bounds(Rectangle(2, 2, 5, 5))
['B']
>>> sorted(g.bounds(Rectangle(0, 0, 5, 5)))
['A', 'B']
"""
x_range = (bounds.left, bounds.right)
y_range = (bounds.lower, bounds.upper)
items = []
lower_left = self.__grid_loc((x_range[0], y_range[0]))
upper_right = self.__grid_loc((x_range[1], y_range[1]))
for i in range(lower_left[0], upper_right[0] + 1):
for j in range(lower_left[1], upper_right[1] + 1):
if (i, j) in self.hash:
items.extend(
[
item[1]
for item in [
item
for item in self.hash[(i, j)]
if x_range[0] <= item[0][0] <= x_range[1]
and y_range[0] <= item[0][1] <= y_range[1]
]
]
)
return items
[docs]
def proximity(self, pt, r):
"""
Returns a list of items found in the grid
within a specified distance of a point.
proximity(Point, number) -> x list
Parameters
----------
pt : the location to search around
r : the distance to search around the point
Examples
--------
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
>>> g.add('A', Point((1.0, 1.0)))
'A'
>>> g.add('B', Point((4.0, 4.0)))
'B'
>>> g.proximity(Point((2.0, 1.0)), 2)
['A']
>>> g.proximity(Point((6.0, 5.0)), 3.0)
['B']
>>> sorted(g.proximity(Point((4.0, 1.0)), 4.0))
['A', 'B']
"""
items = []
lower_left = self.__grid_loc((pt[0] - r, pt[1] - r))
upper_right = self.__grid_loc((pt[0] + r, pt[1] + r))
for i in range(lower_left[0], upper_right[0] + 1):
for j in range(lower_left[1], upper_right[1] + 1):
if (i, j) in self.hash:
items.extend(
[
item[1]
for item in [
item
for item in self.hash[(i, j)]
if get_points_dist(pt, item[0]) <= r
]
]
)
return items
[docs]
def nearest(self, pt):
"""
Returns the nearest item to a point.
nearest(Point) -> x
Parameters
----------
pt : the location to search near
Examples
--------
>>> g = Grid(Rectangle(0, 0, 10, 10), 1)
>>> g.add('A', Point((1.0, 1.0)))
'A'
>>> g.add('B', Point((4.0, 4.0)))
'B'
>>> g.nearest(Point((2.0, 1.0)))
'A'
>>> g.nearest(Point((7.0, 5.0)))
'B'
"""
search_size = self.res
while self.proximity(pt, search_size) == [] and (
get_points_dist((self.x_range[0], self.y_range[0]), pt) > search_size
or get_points_dist((self.x_range[1], self.y_range[0]), pt) > search_size
or get_points_dist((self.x_range[0], self.y_range[1]), pt) > search_size
or get_points_dist((self.x_range[1], self.y_range[1]), pt) > search_size
):
search_size = 2 * search_size
items = []
lower_left = self.__grid_loc((pt[0] - search_size, pt[1] - search_size))
upper_right = self.__grid_loc((pt[0] + search_size, pt[1] + search_size))
for i in range(lower_left[0], upper_right[0] + 1):
for j in range(lower_left[1], upper_right[1] + 1):
if (i, j) in self.hash:
items.extend(
[
(get_points_dist(pt, item[0]), item[1])
for item in self.hash[(i, j)]
]
)
if items == []:
return None
return min(items)[1]
class BruteForcePointLocator:
"""
A class which does naive linear search on a set of Point objects.
"""
def __init__(self, points):
"""
Creates a naive index of the points specified.
__init__(Point list) -> BruteForcePointLocator
Parameters
----------
points : a list of points to index (Point list)
Examples
--------
>>> pl = BruteForcePointLocator([Point((0, 0)), Point((5, 0)), Point((0, 10))])
"""
warnings.warn("BruteForcePointLocator " + dep_msg, FutureWarning)
self._points = points
def nearest(self, query_point):
"""
Returns the nearest point indexed to a query point.
nearest(Point) -> Point
Parameters
----------
query_point : a point to find the nearest indexed point to
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = BruteForcePointLocator(points)
>>> n = pl.nearest(Point((1, 1)))
>>> str(n)
'(0.0, 0.0)'
"""
return min(self._points, key=lambda p: get_points_dist(p, query_point))
def region(self, region_rect):
"""
Returns the indexed points located inside a rectangular query region.
region(Rectangle) -> Point list
Parameters
----------
region_rect : the rectangular range to find indexed points in
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = BruteForcePointLocator(points)
>>> pts = pl.region(Rectangle(-1, -1, 10, 10))
>>> len(pts)
3
"""
return [
p
for p in self._points
if get_rectangle_point_intersect(region_rect, p) is not None
]
def proximity(self, origin, r):
"""
Returns the indexed points located within some distance of an origin point.
proximity(Point, number) -> Point list
Parameters
----------
origin : the point to find indexed points near
r : the maximum distance to find indexed point from the origin point
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = BruteForcePointLocator(points)
>>> neighs = pl.proximity(Point((1, 0)), 2)
>>> len(neighs)
1
>>> p = neighs[0]
>>> isinstance(p, Point)
True
>>> str(p)
'(0.0, 0.0)'
"""
return [p for p in self._points if get_points_dist(p, origin) <= r]
[docs]
class PointLocator:
"""
An abstract representation of a point indexing data structure.
"""
[docs]
def __init__(self, points):
"""
Returns a point locator object.
__init__(Point list) -> PointLocator
Parameters
----------
points : a list of points to index
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = PointLocator(points)
"""
warnings.warn("PointLocator " + dep_msg, FutureWarning)
self._locator = BruteForcePointLocator(points)
[docs]
def nearest(self, query_point):
"""
Returns the nearest point indexed to a query point.
nearest(Point) -> Point
Parameters
----------
query_point : a point to find the nearest indexed point to
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = PointLocator(points)
>>> n = pl.nearest(Point((1, 1)))
>>> str(n)
'(0.0, 0.0)'
"""
return self._locator.nearest(query_point)
[docs]
def region(self, region_rect):
"""
Returns the indexed points located inside a rectangular query region.
region(Rectangle) -> Point list
Parameters
----------
region_rect : the rectangular range to find indexed points in
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = PointLocator(points)
>>> pts = pl.region(Rectangle(-1, -1, 10, 10))
>>> len(pts)
3
"""
return self._locator.region(region_rect)
overlapping = region
[docs]
def polygon(self, polygon):
"""
Returns the indexed points located inside a polygon
"""
# get points in polygon bounding box
# for points in bounding box, check for inclusion in polygon
[docs]
def proximity(self, origin, r):
"""
Returns the indexed points located within some distance of an origin point.
proximity(Point, number) -> Point list
Parameters
----------
origin : the point to find indexed points near
r : the maximum distance to find indexed point from the origin point
Examples
--------
>>> points = [Point((0, 0)), Point((1, 6)), Point((5.4, 1.4))]
>>> pl = PointLocator(points)
>>> len(pl.proximity(Point((1, 0)), 2))
1
"""
return self._locator.proximity(origin, r)
[docs]
class PolygonLocator:
"""
An abstract representation of a polygon indexing data structure.
"""
[docs]
def __init__(self, polygons):
"""
Returns a polygon locator object.
__init__(Polygon list) -> PolygonLocator
Parameters
----------
polygons : a list of polygons to index
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> pl = PolygonLocator([p1, p2])
>>> isinstance(pl, PolygonLocator)
True
"""
warnings.warn("PolygonLocator " + dep_msg, FutureWarning)
self._locator = polygons
# create and rtree
self._rtree = RTree()
for polygon in polygons:
x = polygon.bounding_box.left
y = polygon.bounding_box.lower
_x = polygon.bounding_box.right
_y = polygon.bounding_box.upper
self._rtree.insert(polygon, Rect(x, y, _x, _y))
[docs]
def inside(self, query_rectangle):
"""
Returns polygons that are inside query_rectangle
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> p3 = Polygon([Point((7, 1)), Point((8, 7)), Point((9, 1))])
>>> pl = PolygonLocator([p1, p2, p3])
>>> qr = Rectangle(0, 0, 5, 5)
>>> res = pl.inside( qr )
>>> len(res)
1
>>> qr = Rectangle(3, 7, 5, 8)
>>> res = pl.inside( qr )
>>> len(res)
0
>>> qr = Rectangle(10, 10, 12, 12)
>>> res = pl.inside( qr )
>>> len(res)
0
>>> qr = Rectangle(0, 0, 12, 12)
>>> res = pl.inside( qr )
>>> len(res)
3
Notes
-----
inside means the intersection of the query rectangle and a
polygon is not empty and is equal to the area of the polygon
"""
left = query_rectangle.left
right = query_rectangle.right
upper = query_rectangle.upper
lower = query_rectangle.lower
# rtree rect
qr = Rect(left, lower, right, upper)
# bb overlaps
res = [r.leaf_obj() for r in self._rtree.query_rect(qr) if r.is_leaf()]
qp = Polygon(
[
Point((left, lower)),
Point((right, lower)),
Point((right, upper)),
Point((left, upper)),
]
)
ip = []
gppi = get_polygon_point_intersect
for poly in res:
lower = poly.bounding_box.lower
right = poly.bounding_box.right
upper = poly.bounding_box.upper
left = poly.bounding_box.left
p1 = Point((left, lower))
p2 = Point((right, upper))
if gppi(qp, p1) and gppi(qp, p2):
ip.append(poly)
return ip
[docs]
def overlapping(self, query_rectangle):
"""
Returns list of polygons that overlap query_rectangle
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> p3 = Polygon([Point((7, 1)), Point((8, 7)), Point((9, 1))])
>>> pl = PolygonLocator([p1, p2, p3])
>>> qr = Rectangle(0, 0, 5, 5)
>>> res = pl.overlapping( qr )
>>> len(res)
2
>>> qr = Rectangle(3, 7, 5, 8)
>>> res = pl.overlapping( qr )
>>> len(res)
1
>>> qr = Rectangle(10, 10, 12, 12)
>>> res = pl.overlapping( qr )
>>> len(res)
0
>>> qr = Rectangle(0, 0, 12, 12)
>>> res = pl.overlapping( qr )
>>> len(res)
3
>>> qr = Rectangle(8, 3, 9, 4)
>>> p1 = Polygon([Point((2, 1)), Point((2, 3)), Point((4, 3)), Point((4,1))])
>>> p2 = Polygon([Point((7, 1)), Point((7, 5)), Point((10, 5)), Point((10, 1))])
>>> pl = PolygonLocator([p1, p2])
>>> res = pl.overlapping(qr)
>>> len(res)
1
Notes
-----
overlapping means the intersection of the query rectangle and a
polygon is not empty and is no larger than the area of the polygon
"""
left = query_rectangle.left
right = query_rectangle.right
upper = query_rectangle.upper
lower = query_rectangle.lower
# rtree rect
qr = Rect(left, lower, right, upper)
# bb overlaps
res = [r.leaf_obj() for r in self._rtree.query_rect(qr) if r.is_leaf()]
# have to check for polygon overlap using segment intersection
# add polys whose bb contains at least one of the corners of the query
# rectangle
sw = (left, lower)
se = (right, lower)
ne = (right, upper)
nw = (left, upper)
pnts = [sw, se, ne, nw]
cs = []
for pnt in pnts:
c = [r.leaf_obj() for r in self._rtree.query_point(pnt) if r.is_leaf()]
cs.extend(c)
cs = list(set(cs))
overlapping = []
# first find polygons with at least one vertex inside query rectangle
remaining = copy.copy(res)
for polygon in res:
vertices = polygon.vertices
for vertex in vertices:
xb = vertex[0] >= left
xb *= vertex[0] < right
yb = vertex[1] >= lower
yb *= vertex[1] < upper
if xb * yb:
overlapping.append(polygon)
remaining.remove(polygon)
break
# for remaining polys in bb overlap check if vertex chains intersect
# segments of the query rectangle
left_edge = LineSegment(Point((left, lower)), Point((left, upper)))
right_edge = LineSegment(Point((right, lower)), Point((right, upper)))
lower_edge = LineSegment(Point((left, lower)), Point((right, lower)))
upper_edge = LineSegment(Point((left, upper)), Point((right, upper)))
for polygon in remaining:
vertices = copy.copy(polygon.vertices)
if vertices[-1] != vertices[0]:
vertices.append(vertices[0]) # put on closed cartographic form
nv = len(vertices)
for i in range(nv - 1):
head = vertices[i]
tail = vertices[i + 1]
edge = LineSegment(head, tail)
li = get_segments_intersect(edge, left_edge)
if (
li
or get_segments_intersect(edge, right_edge)
or get_segments_intersect(edge, lower_edge)
or get_segments_intersect(edge, upper_edge)
):
overlapping.append(polygon)
break
# check remaining for explicit containment of the bounding rectangle
# cs has candidates for this check
sw = Point(sw)
se = Point(se)
ne = Point(ne)
nw = Point(nw)
for polygon in cs:
if (
get_polygon_point_intersect(polygon, sw)
or get_polygon_point_intersect(polygon, se)
or get_polygon_point_intersect(polygon, ne)
or get_polygon_point_intersect(polygon, nw)
):
overlapping.append(polygon)
break
return list(set(overlapping))
[docs]
def nearest(self, query_point, rule="vertex"):
"""
Returns the nearest polygon indexed to a query point based on
various rules.
nearest(Polygon) -> Polygon
Parameters
----------
query_point : a point to find the nearest indexed polygon to
rule : representative point for polygon in nearest query.
vertex -- measures distance between vertices and query_point
centroid -- measures distance between centroid and
query_point
edge -- measures the distance between edges and query_point
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> pl = PolygonLocator([p1, p2])
>>> try: n = pl.nearest(Point((-1, 1)))
... except NotImplementedError: print("future test: str(min(n.vertices())) == (0.0, 1.0)")
future test: str(min(n.vertices())) == (0.0, 1.0)
""" # noqa: E501
raise NotImplementedError
[docs]
def region(self, region_rect):
"""
Returns the indexed polygons located inside a rectangular query region.
region(Rectangle) -> Polygon list
Parameters
----------
region_rect : the rectangular range to find indexed polygons in
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> pl = PolygonLocator([p1, p2])
>>> n = pl.region(Rectangle(0, 0, 4, 10))
>>> len(n)
2
"""
n = self._locator
for polygon in n:
points = polygon.vertices
pl = BruteForcePointLocator(points)
pts = pl.region(region_rect)
if len(pts) == 0:
n.remove(polygon)
return n
[docs]
def contains_point(self, point):
"""
Returns polygons that contain point
Parameters
----------
point: point (x,y)
Returns
-------
list of polygons containing point
Examples
--------
>>> p1 = Polygon([Point((0,0)), Point((6,0)), Point((4,4))])
>>> p2 = Polygon([Point((1,2)), Point((4,0)), Point((4,4))])
>>> p1.contains_point((2,2))
1
>>> p2.contains_point((2,2))
1
>>> pl = PolygonLocator([p1, p2])
>>> len(pl.contains_point((2,2)))
2
>>> p2.contains_point((1,1))
0
>>> p1.contains_point((1,1))
1
>>> len(pl.contains_point((1,1)))
1
>>> p1.centroid
(3.3333333333333335, 1.3333333333333333)
>>> pl.contains_point((1,1))[0].centroid
(3.3333333333333335, 1.3333333333333333)
"""
# bbounding box containment
res = [r.leaf_obj() for r in self._rtree.query_point(point) if r.is_leaf()]
# explicit containment check for candidate polygons needed
return [poly for poly in res if poly.contains_point(point)]
[docs]
def proximity(self, origin, r, rule="vertex"):
"""
Returns the indexed polygons located within some distance of an
origin point based on various rules.
proximity(Polygon, number) -> Polygon list
Parameters
----------
origin : the point to find indexed polygons near
r : the maximum distance to find indexed polygon from the origin point
rule : representative point for polygon in nearest query.
vertex -- measures distance between vertices and query_point
centroid -- measures distance between centroid and
query_point
edge -- measures the distance between edges and query_point
Examples
--------
>>> p1 = Polygon([Point((0, 1)), Point((4, 5)), Point((5, 1))])
>>> p2 = Polygon([Point((3, 9)), Point((6, 7)), Point((1, 1))])
>>> pl = PolygonLocator([p1, p2])
>>> try:
... len(pl.proximity(Point((0, 0)), 2))
... except NotImplementedError:
... print("future test: len(pl.proximity(Point((0, 0)), 2)) == 2")
future test: len(pl.proximity(Point((0, 0)), 2)) == 2
"""
raise NotImplementedError