Repository: pygeoops/pygeoops (language: python)

Issue:
ENH: Add option to `centerline` to extend it to the polygon boundary
Add an option to the `centerline` function so it gets extended to the polygon boundary.

Once released, notify here: https://gis.stackexchange.com/questions/482520/drawing-a-centerline-through-a-polygon-using-python-libraries


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6f54bc3..c36f5f2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 ### Improvements
 
+- Add functions `extend_line_by_distance` and `extend_line_to_polygon` (#83)
 - Improve logging of geos warnings in `difference_all` (#77)
 - Use ruff instead of black for formatting + use mypy (#78)
 
diff --git a/docs/reference.rst b/docs/reference.rst
index b2a9555..8868520 100644
--- a/docs/reference.rst
+++ b/docs/reference.rst
@@ -21,6 +21,8 @@ methods implemented in PyGeoOps. All classes and function exposed in
    difference_all_tiled
    empty
    explode
+   extend_line_by_distance
+   extend_line_to_polygon
    make_valid
    remove_inner_rings
    simplify
diff --git a/pygeoops/__init__.py b/pygeoops/__init__.py
index 21bdb59..453f9c5 100644
--- a/pygeoops/__init__.py
+++ b/pygeoops/__init__.py
@@ -2,6 +2,7 @@
 
 from pygeoops._difference import *  # noqa: F403
 from pygeoops._centerline import *  # noqa: F403
+from pygeoops._extend_line import *  # noqa: F403
 from pygeoops._general import *  # noqa: F403
 from pygeoops._grid import *  # noqa: F403
 from pygeoops._simplify import *  # noqa: F403
diff --git a/pygeoops/_extend_line.py b/pygeoops/_extend_line.py
new file mode 100644
index 0000000..a614f57
--- /dev/null
+++ b/pygeoops/_extend_line.py
@@ -0,0 +1,227 @@
+import math
+from typing import Tuple
+import shapely
+from shapely import box, LineString, Polygon, Point
+
+
+def extend_line_by_distance(
+    line: LineString, start_distance: float, end_distance: float
+) -> LineString:
+    """
+    Extends a line with the specified distances.
+
+    Args:
+        line (LineString): The line to extend.
+        start_distance (float): The distance to extend the start of the line
+            with.
+        end_distance (float): The distance to extend the end of the line with.
+
+    Returns:
+        LineString: The extended line.
+    """
+    if start_distance == 0 and end_distance == 0:
+        return line
+
+    line_coords = shapely.get_coordinates(line)
+    line_result = line_coords.copy()
+
+    # Extend the start and end of the line
+    _, line_result[0] = _extend_segment_by_distance(
+        line_coords[1], line_coords[0], start_distance
+    )
+    _, line_result[-1] = _extend_segment_by_distance(
+        line_coords[-2], line_coords[-1], end_distance
+    )
+
+    return LineString(line_result)
+
+
+def extend_line_to_polygon(line: LineString, geometry: Polygon) -> LineString:
+    """
+    Extends a line to a given Polygon.
+
+    Args:
+        line (LineString): The line to extend.
+        geometry (Polygon): The Polygon to extend the line to.
+
+    Returns:
+        LineString: The extended line.
+    """
+    line_coords = shapely.get_coordinates(line)
+    line_result = line_coords.copy()
+    geometry_line = LineString(shapely.get_coordinates(geometry))
+    shapely.prepare(geometry_line)
+
+    # Find the closest point on geometry_line that the start of the line can extend to.
+    line_result[0] = _find_closest_extend_point(
+        line_coords[1], line_coords[0], geometry_line
+    )
+
+    # Find the closest point on geometry_line that the end of the line can extend to.
+    line_result[-1] = _find_closest_extend_point(
+        line_coords[-2], line_coords[-1], geometry_line
+    )
+
+    return LineString(line_result)
+
+
+def _find_closest_extend_point(
+    p1: Tuple[float, float],
+    p2: Tuple[float, float],
+    geometry_line: LineString,
+) -> Tuple[float, float]:
+    """
+    Find the closest point on the geometry line that can be extended to.
+
+    Args:
+        p1 (Tuple[float, float]): point 1 of the segment to extend.
+        p2 (Tuple[float, float]): point 2 of the segment to extend.
+        geometry_line (LineString): the line to extend to.
+
+    Returns:
+        Tuple[float, float]: the closest point on the geometry_line that can be
+            extended to.
+    """
+    # If p2 is already on the geometry line, return p2.
+    if geometry_line.intersects(Point(p2)):
+        return p2
+
+    # Extend the segment to the bounding box of the geometry.
+    _, p2_ext = _extend_segment_to_bbox(p1, p2, bbox=geometry_line.bounds)
+
+    # We want to find the intersection point of the extension with the geometry lines
+    # that is closest to p2.
+    end_extension = (p2, p2_ext)
+    intersections = geometry_line.intersection(LineString(end_extension))
+    intersection_points = shapely.get_coordinates(intersections)
+
+    if len(intersection_points) == 0:
+        # No intersection found, so return input point p2.
+        return p2
+    else:
+        # Find the intersection point closest to input point p2.
+        intersection_points = shapely.points(intersection_points)
+        p2_point = Point(p2)
+        points_sorted_by_distance = sorted(intersection_points, key=p2_point.distance)
+        return points_sorted_by_distance[0].coords[0]
+
+
+def _extend_segment_by_distance(
+    p1: Tuple[float, float], p2: Tuple[float, float], distance: float
+) -> Tuple[Tuple[float, float], Tuple[float, float]]:
+    """
+    Extends the segment by the distance specified.
+
+    Args:
+        p1 (Tuple[float, float]): the first point of the segment
+        p2 (Tuple[float, float]): the second point of the segment
+        distance (float): the distance to extend the segment with. If distance
+            < 0, the segment will be extended in the negative direction.
+
+    Returns:
+        Tuple[Tuple[float, float], Tuple[float, float]]: the segment extended by
+            the distance specified.
+    """
+    if distance < 0:
+        raise ValueError(f"distance must be >= 0, received: {distance}")
+
+    segment_length = math.sqrt(pow(p1[0] - p2[0], 2.0) + pow(p1[1] - p2[1], 2.0))
+    if segment_length == 0:
+        raise ValueError("lenght of input segment cannot be 0")
+    ratio = distance / segment_length
+    return _extend_segment_by_ratio(p1, p2, ratio=ratio)
+
+
+def _extend_segment_by_ratio(
+    p1: Tuple[float, float], p2: Tuple[float, float], ratio: float
+) -> Tuple[Tuple[float, float], Tuple[float, float]]:
+    """
+    Extends the segment by a ratio.
+
+    Args:
+        p1 (Tuple[float, float]): the first point of the segment
+        p2 (Tuple[float, float]): the second point of the segment
+        ratio (float): the ratio to extend the segment with. Ratio must be >= 0.
+
+    Returns:
+        Tuple[Tuple[float, float], Tuple[float, float]]: the segment extended by
+            the ratio specified.
+    """
+    if ratio < 0:
+        raise ValueError(f"ratio must be >= 0, received: {ratio}")
+
+    x_ext = p2[0] + (p2[0] - p1[0]) * ratio
+    y_ext = p2[1] + (p2[1] - p1[1]) * ratio
+    p_ext = (x_ext, y_ext)
+
+    return (p1, p_ext)
+
+
+def _extend_segment_to_bbox(
+    p1: Tuple[float, float],
+    p2: Tuple[float, float],
+    bbox: Tuple[float, float, float, float],
+) -> Tuple[Tuple[float, float], Tuple[float, float]]:
+    """
+    Extends a segment so both points are onto the boundaries of a given bounding
+    box.
+
+    Args:
+        p1 (Tuple[float, float]): The first point of the segment.
+        p2 (Tuple[float, float]): The second point of the segment.
+        bbox (Tuple[float, float, float, float]): A tuple representing the
+            boundary of the bounding box to extend to in the format
+            (minx, miny, maxx, maxy).
+
+    Returns:
+        Tuple[Tuple[float, float], Tuple[float, float]]: The extended segment on
+            the boundary of the bbox.
+    """
+    minx, miny, maxx, maxy = bbox
+    if p1[0] == p2[0]:  # vertical segment
+        if p1[1] < p2[1]:
+            # Keep segment from bottom to top
+            return ((p1[0], miny), (p1[0], maxy))
+        else:
+            # Keep segment from top to bottom
+            return ((p1[0], maxy), (p1[0], miny))
+    elif p1[1] == p2[1]:  # horizontal segment
+        if p1[0] < p2[0]:
+            # Keep segment from left to right
+            return ((minx, p1[1]), (maxx, p1[1]))
+        else:
+            # Keep segment from right to left
+            return ((maxx, p1[1]), (minx, p1[1]))
+    else:
+        # linear equation: y = k*x + m
+        k = (p2[1] - p1[1]) / (p2[0] - p1[0])
+        m = p1[1] - k * p1[0]
+        y0 = k * minx + m
+        y1 = k * maxx + m
+        x0 = (miny - m) / k
+        x1 = (maxy - m) / k
+        points_on_boundary_lines = [
+            Point(minx, y0),
+            Point(maxx, y1),
+            Point(x0, miny),
+            Point(x1, maxy),
+        ]
+        # Sort the points to the points found closest to the bbox are first
+        bbox_geom = box(minx, miny, maxx, maxy)
+        points_sorted_by_distance = sorted(
+            points_on_boundary_lines, key=bbox_geom.distance
+        )
+
+        # Make sure the point belonging to the first input param is first in result.
+        p_extended1 = next(iter(points_sorted_by_distance[0].coords))
+        p_extended2 = next(iter(points_sorted_by_distance[1].coords))
+        if p1[0] < p2[0]:
+            if p_extended1[0] < p_extended2[0]:
+                return (p_extended1, p_extended2)
+            else:
+                return (p_extended2, p_extended1)
+        else:
+            if p_extended1[0] > p_extended2[0]:
+                return (p_extended1, p_extended2)
+            else:
+                return (p_extended2, p_extended1)


Test patch:
diff --git a/tests/test_extend_line.py b/tests/test_extend_line.py
new file mode 100644
index 0000000..7c8621a
--- /dev/null
+++ b/tests/test_extend_line.py
@@ -0,0 +1,172 @@
+import math
+import pytest
+from shapely import box, LineString, Polygon
+
+import pygeoops
+from pygeoops import _extend_line
+
+
+@pytest.mark.parametrize(
+    "desc, line, start_distance, end_distance, exp_line",
+    [
+        (
+            "line with 1 diagonal segment, start_distance != end_distance",
+            LineString([(1, 1), (2, 2)]),
+            math.sqrt(2),
+            2 * math.sqrt(2),
+            LineString([(0, 0), (4, 4)]),
+        ),
+        (
+            "line with 1 diagonal segment, start_distance = end_distance",
+            LineString([(1, 1), (2, 2)]),
+            math.sqrt(2),
+            math.sqrt(2),
+            LineString([(0, 0), (3, 3)]),
+        ),
+        (
+            "line with 3 segments, start and end vertical",
+            LineString([(2, 3), (2, 2), (3, 2), (3, 1)]),
+            1,
+            2,
+            LineString([(2, 4), (2, 2), (3, 2), (3, -1)]),
+        ),
+        (
+            "0 start_distance and end_distance",
+            LineString([(1, 1), (2, 2)]),
+            0,
+            0,
+            LineString([(1, 1), (2, 2)]),
+        ),
+    ],
+)
+def test_extend_line_by_distance(desc, line, start_distance, end_distance, exp_line):
+    result = pygeoops.extend_line_by_distance(line, start_distance, end_distance)
+    assert result == exp_line
+
+
+@pytest.mark.parametrize(
+    "desc, line, poly, exp_line",
+    [
+        (
+            "each line extension intersects with one point",
+            LineString([(4, 3), (5, 5), (6, 5)]),
+            box(0, 0, 10, 10),
+            LineString([(2.5, 0), (5, 5), (10, 5)]),
+        ),
+        (
+            "one line extension intersects with multiple points",
+            LineString([(4, 3), (5, 5), (6, 5)]),
+            Polygon(
+                [(3, 0), (3, 10), (7, 10), (7, 3), (10, 10), (10, 0), (3, 0)],
+            ),
+            LineString([(3, 1), (5, 5), (7, 5)]),
+        ),
+        (
+            "one line extension intersects with a point + a line",
+            LineString([(4, 3), (5, 5), (6, 5)]),
+            Polygon(
+                [(3, 0), (3, 10), (7, 10), (7, 5), (8, 5), (10, 10), (10, 0), (3, 0)],
+            ),
+            LineString([(3, 1), (5, 5), (7, 5)]),
+        ),
+        (
+            "input line is already on the polygon",
+            LineString([(3, 1), (5, 5), (7, 5)]),
+            Polygon(
+                [(3, 0), (3, 10), (7, 10), (7, 5), (8, 5), (10, 10), (10, 0), (3, 0)],
+            ),
+            LineString([(3, 1), (5, 5), (7, 5)]),
+        ),
+        (
+            "input line is already on the polygon",
+            LineString([(4, 3), (5, 5), (6, 5)]),
+            box(4, 0, 5, 1),
+            LineString([(4, 3), (5, 5), (6, 5)]),
+        ),
+    ],
+)
+def test_extend_line_to_polygon(desc, line, poly, exp_line):
+    result = pygeoops.extend_line_to_polygon(line, poly)
+    assert result == exp_line
+
+
+@pytest.mark.parametrize(
+    "desc, p1, p2, distance, exp_p",
+    [
+        ("bl-tr, distance 0", (0, 0), (1, 1), 0, ((0, 0), (1, 1))),
+        ("bl-tr, distance 1.44", (0, 0), (1, 1), math.sqrt(2), ((0, 0), (2, 2))),
+        ("bl-tr, distance 2.88", (0, 0), (1, 1), 2 * math.sqrt(2), ((0, 0), (3, 3))),
+        ("horizontal, distance 1", (0, 0), (1, 0), 1, ((0, 0), (2, 0))),
+        ("vertical, distance 1", (0, 0), (0, 1), 1, ((0, 0), (0, 2))),
+    ],
+)
+def test_extend_segment_by_distance(desc, p1, p2, distance, exp_p):
+    result = _extend_line._extend_segment_by_distance(p1, p2, distance)
+    assert result == exp_p
+
+
+@pytest.mark.parametrize(
+    "exp_error, p1, p2, distance",
+    [
+        ("distance must be >= 0", (0, 0), (1, 1), -1),
+        ("lenght of input segment cannot be 0", (0, 0), (0, 0), 1),
+    ],
+)
+def test_extend_segment_by_distance_invalid(exp_error, p1, p2, distance):
+    with pytest.raises(ValueError, match=exp_error):
+        _ = _extend_line._extend_segment_by_distance(p1, p2, distance)
+
+
+@pytest.mark.parametrize(
+    "desc, p1, p2, ratio, exp_p",
+    [
+        ("diagonal bl-tr, ratio 0", (0, 0), (1, 1), 0, ((0, 0), (1, 1))),
+        ("diagonal bl-tr, ratio 1", (0, 0), (1, 1), 1, ((0, 0), (2, 2))),
+        ("diagonal bl-tr, ratio 0.5", (0, 0), (1, 1), 0.5, ((0, 0), (1.5, 1.5))),
+        ("diagonal tr-bl, ratio 0", (1, 1), (0, 0), 0, ((1, 1), (0, 0))),
+        ("diagonal tr-bl, ratio 1", (1, 1), (0, 0), 1, ((1, 1), (-1, -1))),
+        ("diagonal tr-bl, ratio 0.5", (1, 1), (0, 0), 0.5, ((1, 1), (-0.5, -0.5))),
+        ("diagonal br-tl, ratio 1", (1, 0), (0, 1), 1, ((1, 0), (-1, 2))),
+        ("diagonal tl-br, ratio 1", (0, 1), (1, 0), 1, ((0, 1), (2, -1))),
+        ("horizonal, ratio 0", (0, 0), (1, 0), 0, ((0, 0), (1, 0))),
+        ("horizonal, ratio 1", (0, 0), (1, 0), 1, ((0, 0), (2, 0))),
+        ("horizonal, ratio 0.5", (0, 0), (1, 0), 0.5, ((0, 0), (1.5, 0))),
+        ("vertical, ratio 0", (0, 0), (0, 1), 0, ((0, 0), (0, 1))),
+        ("vertical, ratio 1", (0, 0), (0, 1), 1, ((0, 0), (0, 2))),
+        ("vertical, ratio 0.5", (0, 0), (0, 1), 0.5, ((0, 0), (0, 1.5))),
+    ],
+)
+def test_extend_segment_by_ratio(desc, p1, p2, ratio, exp_p):
+    result = _extend_line._extend_segment_by_ratio(p1, p2, ratio)
+    assert result == exp_p
+
+
+@pytest.mark.parametrize(
+    "exp_error, p1, p2, ratio",
+    [
+        ("ratio must be >= 0", (0, 0), (1, 1), -1),
+    ],
+)
+def test_extend_segment_by_ratio_invalid(exp_error, p1, p2, ratio):
+    with pytest.raises(ValueError, match=exp_error):
+        _ = _extend_line._extend_segment_by_ratio(p1, p2, ratio)
+
+
+@pytest.mark.parametrize(
+    "desc, p1, p2, bbox, exp_line",
+    [
+        ("diag, left-right", (1, 1), (2, 2), (0, 0, 4, 4), ((0, 0), (4, 4))),
+        ("diag, left-right, switched", (1, 1), (2, 2), (4, 4, 0, 0), ((0, 0), (4, 4))),
+        ("diag, left-right, touch", (2, 2), (4, 4), (0, 0, 4, 4), ((0, 0), (4, 4))),
+        ("diag, right-left", (2, 2), (1, 1), (0, 0, 4, 4), ((4, 4), (0, 0))),
+        ("diag, right-left, touch", (4, 4), (2, 2), (0, 0, 4, 4), ((4, 4), (0, 0))),
+        ("diag, right-left, beyond", (5, 5), (2, 2), (0, 0, 4, 4), ((4, 4), (0, 0))),
+        ("horizontal, left-right", (1, 1), (2, 1), (0, 0, 4, 4), ((0, 1), (4, 1))),
+        ("horizontal, right-left", (2, 1), (1, 1), (0, 0, 4, 4), ((4, 1), (0, 1))),
+        ("vertical, bottom-top", (1, 1), (1, 2), (0, 0, 4, 4), ((1, 0), (1, 4))),
+        ("vertical, top-bottom", (1, 2), (1, 1), (0, 0, 4, 4), ((1, 4), (1, 0))),
+    ],
+)
+def test_extend_segment_to_bbox(desc, p1, p2, bbox, exp_line):
+    result = _extend_line._extend_segment_to_bbox(p1, p2, bbox)
+    assert result == exp_line