Directional Distance with Shapely

directional_distance.py

"""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