Last active
January 3, 2021 04:32
-
-
Save Jeremiah-England/6377b798d6384d454583e55fc053f3cc to your computer and use it in GitHub Desktop.
Directional Distance with Shapely
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
"""A small module for finding the distance between two geometries in an arbitrary direction | |
Inspired by a question/issue in the Shapely github repository (`#1051`_). | |
Notes | |
----- | |
There is another way I can think of doing this which will probably be much faster for | |
large geometries, but less accurate. It should be possible to stretch both geometries | |
by a very large amount in the direction perpendicular to the direction you want to | |
find the distance in. Of course, you would loose the 'movement' idea (i.e. it could | |
return what should be a 'negative' distance in that direction.) But for some use cases | |
I think this method would be preferred because there is only one transform and then a | |
`distance` call. | |
Also, I should note that the method below can easily be adapted to return 'negative' | |
distances or 'absolute' distances in the given direction. | |
.. _#1051: https://github.com/Toblerity/Shapely/issues/1051/ | |
""" | |
from shapely.geometry import Point, Polygon, LineString | |
from shapely.geometry.base import BaseGeometry, BaseMultipartGeometry | |
import math | |
from typing import Tuple | |
def max_dist_upper_bound(a: BaseGeometry, b: BaseGeometry) -> float: | |
"""Calculate an upper bound for the maximum distance between two geometries using their bounding boxes | |
Parameters | |
---------- | |
a : BaseGeometry | |
first shapely geometry | |
b : BaseGeometry | |
second shapely geometry | |
Notes | |
----- | |
Perhaps a simpler way to do this is put all the geometries into a | |
collection and take the bounding box of that. The distance between | |
two opposite corners would do for an upper bound, though a bit looser | |
but that doesn't matter. | |
Returns | |
------- | |
float | |
An upper bound for the maximum distance between the two geometries using the | |
furthest corners of their bounding boxes. | |
""" | |
ea = a.envelope | |
eb = b.envelope | |
if isinstance(ea, Point): | |
if isinstance(eb, Point): | |
return ea.distance(eb) | |
else: | |
return max(ea.distance(Point(x, y)) for x, y in eb.exterior.coords) | |
else: | |
max_dist = -math.inf | |
for a_corner in ea.exterior.coords: | |
for b_corner in eb.exterior.coords: | |
dist = math.dist(a_corner, b_corner) | |
max_dist = max(max_dist, dist) | |
if max_dist < 0: | |
ValueError(f'Should never try to return a negative number: {max_dist}') | |
return max_dist | |
def directional_distance(moving_geom: BaseGeometry, stationary_geom: BaseGeometry, | |
direction: Tuple[float, float]) -> float: | |
"""Find the distance between two geometries in a given direction | |
Parameters | |
---------- | |
moving_geom : shapely geometry | |
geometry moving in the direction | |
stationary_geom : shapely geometry | |
geometry being moved toward by the other | |
direction : Tuple[float, float] | |
The direction the `moving_geom` is moving (unit vector) | |
Returns | |
------- | |
float | |
The distance `moving_geom` can move in `direction` before colliding with `stationary_geom` | |
If the geometries never collide then math.inf is returned. | |
""" | |
if not isinstance(moving_geom, BaseGeometry): | |
raise TypeError(f'The `moving_geom` argument must be a shapely geometry: {moving_geom}') | |
if not isinstance(stationary_geom, BaseGeometry): | |
raise TypeError(f'The `stationary_geom` argument must be a shapely geometry: {moving_geom}') | |
if isinstance(moving_geom, BaseMultipartGeometry) or isinstance(stationary_geom, BaseMultipartGeometry): | |
raise NotImplementedError('Cannot yet handle Multi(LineString/Polygon/Point) geometries.') | |
# return early if they already touch | |
if moving_geom.intersects(stationary_geom): | |
return 0 | |
# unitize the direction if it is not already | |
dir_length = math.hypot(*direction) | |
if dir_length != 1: | |
direction = (direction[0] / dir_length, direction[1] / dir_length) | |
ray_length = max_dist_upper_bound(moving_geom, stationary_geom) | |
ray_x = direction[0] * ray_length | |
ray_y = direction[1] * ray_length | |
# get the coordinates of both geometries as the minimum distance | |
# line will start or end at one of the vertexes | |
if isinstance(moving_geom, Polygon): | |
moving_coords = moving_geom.exterior.coords | |
else: | |
moving_coords = moving_geom.coords | |
if isinstance(stationary_geom, Polygon): | |
stationary_coords = stationary_geom.exterior.coords | |
else: | |
stationary_coords = stationary_geom.coords | |
min_dist = math.inf | |
for coord in moving_coords: | |
ray = LineString([coord, (coord[0] + ray_x, coord[1] + ray_y)]) | |
intersection = ray.intersection(stationary_geom) | |
if intersection.is_empty: | |
continue | |
ray_origin = Point(*coord) | |
dist = ray_origin.distance(intersection) | |
min_dist = min(min_dist, dist) | |
for coord in stationary_coords: | |
ray = LineString([coord, (coord[0] - ray_x, coord[1] - ray_y)]) | |
intersection = ray.intersection(moving_geom) | |
if intersection.is_empty: | |
continue | |
ray_origin = Point(*coord) | |
dist = ray_origin.distance(intersection) | |
min_dist = min(min_dist, dist) | |
return min_dist | |
if __name__ == '__main__': | |
from shapely.geometry import box | |
# run some tests | |
square1 = box(0, 0, 1, 1) | |
square2 = box(2, -2, 3, -1) | |
direction_vec = (1, -1) | |
assert directional_distance(square1, square2, direction_vec) == math.dist((1, 0), (2, -1)) | |
direction_vec = (1.5, -1) | |
assert directional_distance(square1, square2, direction_vec) == math.dist((1, 0), (2.5, -1)) | |
direction_vec = (-1, 1) | |
assert directional_distance(square1, square2, direction_vec) == math.inf | |
pt1 = Point(0, 0) | |
pt2 = Point(1, 1) | |
direction_vec = (1, 1) | |
assert directional_distance(pt1, pt2, direction_vec) == math.dist((0, 0), (1, 1)) | |
direction_vec = (1, 2) | |
assert directional_distance(pt1, pt2, direction_vec) == math.inf | |
ls1 = LineString([(0, 0), (0, 1)]) | |
ls2 = LineString([(1, 0), (1, 1)]) | |
direction_vec = (1, 1) | |
assert directional_distance(ls1, ls2, direction_vec) == math.dist((0, 0), (1, 1)) | |
direction_vec = (1, 0) | |
assert directional_distance(ls1, ls2, direction_vec) == math.dist((0, 0), (0, 1)) | |
direction_vec = (-1, -1) | |
assert directional_distance(ls1, ls2, direction_vec) == math.inf |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment