Skip to content
125 changes: 124 additions & 1 deletion geostructures/structures.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,10 +10,12 @@

from abc import ABC
import copy
from datetime import timedelta
from functools import cached_property
import math
import statistics
from typing import cast, Any, Dict, List, Optional, Tuple, Sequence, TYPE_CHECKING
import warnings

import numpy as np

Expand All @@ -23,7 +25,7 @@
_RE_LINESTRING_WKT, LineLikeMixin, PointLikeMixin, PolygonLikeMixin,
SingleShapeBase, SimpleShapeMixin
)
from geostructures.time import GEOTIME_TYPE
from geostructures.time import GEOTIME_TYPE, TimeInterval
from geostructures.coordinates import Coordinate
from geostructures.calc import (
inverse_haversine_radians,
Expand Down Expand Up @@ -1479,6 +1481,127 @@ def intersects_shape(self, shape: 'GeoShape', **kwargs) -> bool:
# because the centroid may fall in a hole
return o_edges[0][0][0] in self or s_edges[0][0][0] in shape

def split(self, distance_meters: float) -> List['GeoLineString']:
"""
Splits a GeoLineString into smaller segments based on a specified distance.

Args:
distance_meters (float): The maximum distance for each segment in meters.

Returns:
List[GeoLineString]: A list of GeoLineString objects, each of which is no longer
than the specified distance. If the total length of the
line is less than the specified distance, the original line
is returned as a single segment.

Notes:
- If the GeoLineString is time-bounded (has a datetime interval), the resulting
segments will have proportional datetime intervals based on the segment's length.
- If the specified distance is greater than the total length of the line,
a warning is issued, and the original line is returned.
"""

out = [] # List to store resulting GeoLineString segments
segments: List[Tuple[Coordinate, Coordinate]] = self.segments.copy() # Copy of all line segments
vertices = [segments[0][0]] # Initialize the first vertex from the starting point of the first segment
remaining_distance_meters = None # Remaining distance from a previous iteration
total_length_meters = sum(haversine_distance_meters(*segment) for segment in segments) # Total line length

# Handle case where the total line length is less than the specified distance
if total_length_meters <= distance_meters:
warnings.warn(
f'Total length ({total_length_meters}) is less than the distance ({distance_meters}); returning line.'
)
return [self]

properties = self._properties.copy() # Copy of properties for the GeoLineString
cumulative_length = 0 # Tracks cumulative length traversed
dt = None # Placeholder for datetime interval for each segment

start_time = end_time = total_duration_seconds = None
# If time interval is provided, calculate total duration in seconds
if self.dt and self.dt.start != self.dt.end:
start_time, end_time = self.dt.start, self.dt.end
total_duration_seconds = (end_time - start_time).total_seconds()

# Iterate through segments and split accordingly
while segments:
remaining_segment_length = haversine_distance_meters(*segments[0]) # Length of the current segment

if remaining_distance_meters is not None:
# If there's a remaining distance from the previous iteration, process it
cumulative_length += distance_meters - remaining_distance_meters
end_point = inverse_haversine_degrees(
vertices[-1],
bearing_degrees(vertices[-1], segments[0][1]),
remaining_distance_meters
)
# Calculate the time interval proportionally if applicable
if total_duration_seconds is not None:
segment_start_time = start_time + timedelta(
seconds=cumulative_length / total_length_meters * total_duration_seconds
)
dt = TimeInterval(dt.end, segment_start_time)
elif self.dt:
dt = self.dt

vertices.append(end_point) # Add the calculated endpoint to vertices
out.append(GeoLineString(vertices, properties=properties.copy(), dt=dt)) # Store the segment
vertices = [end_point] # Reset vertices for the next segment
remaining_segment_length = haversine_distance_meters(vertices[-1], segments[0][1])
remaining_distance_meters = None

while distance_meters < remaining_segment_length:
# Handle cases where the current segment is longer than the specified distance
cumulative_length += distance_meters
end_point = inverse_haversine_degrees(
vertices[-1],
bearing_degrees(vertices[-1], segments[0][1]),
distance_meters
)
# Calculate the time interval proportionally if applicable
if total_duration_seconds is not None:
segment_start_time = start_time + timedelta(
seconds=(cumulative_length - distance_meters)
/ total_length_meters * total_duration_seconds
)
segment_end_time = start_time + timedelta(
seconds=cumulative_length / total_length_meters * total_duration_seconds
)
dt = TimeInterval(segment_start_time, segment_end_time)
elif self.dt:
dt = self.dt

vertices.append(end_point) # Add the calculated endpoint to vertices
out.append(GeoLineString(vertices, properties=properties.copy(), dt=dt)) # Store the segment
vertices = [end_point] # Reset vertices for the next segment
remaining_segment_length = haversine_distance_meters(vertices[-1], segments[0][1])

# Calculate the remaining distance after processing the current segment
remaining_distance_meters = distance_meters - remaining_segment_length
cumulative_length += remaining_distance_meters
vertices.append(segments[0][1]) # Add the endpoint of the current segment

if len(segments) == 1:
break

segments.pop(0) # Remove the processed segment

# Handle the final segment if there are remaining distances
if remaining_distance_meters:
last_segment = out.pop() if out else None
if last_segment:
vertices = last_segment.vertices[:-1] + [vertices[-1]]
# Assign the time interval for the final segment if applicable
if total_duration_seconds is not None:
dt = TimeInterval(last_segment.dt.start, end_time)
elif self.dt:
dt = self.dt

out.append(GeoLineString(vertices, properties=properties.copy(), dt=dt))

return out

def to_geo_interface(self, **kwargs):
return {
**self.__geo_interface__,
Expand Down
80 changes: 80 additions & 0 deletions tests/test_structures.py
Original file line number Diff line number Diff line change
Expand Up @@ -1593,3 +1593,83 @@ def test_geopoint_from_wkt():
def test_geopoint_to_wkt(geopoint):
assert geopoint.to_wkt() == 'POINT(0.0 0.0)'


@pytest.fixture
def basic_line():
# A basic GeoLineString without time
return GeoLineString([
Coordinate(0, 0),
Coordinate(1, 1),
Coordinate(2, 2)
])


@pytest.fixture
def timed_line():
# A GeoLineString with a time interval
return GeoLineString(
[
Coordinate(0, 0),
Coordinate(1, 1),
Coordinate(2, 2)
],
dt=TimeInterval(
start=datetime(2025, 1, 1, 0, 0, 0),
end=datetime(2025, 1, 1, 2, 0, 0)
)
)


def test_split_no_split_needed(basic_line):
# Distance is greater than the total length of the line
result = basic_line.split(200000)
assert len(result) == 1
assert result[0] == basic_line


def test_split_even_segments(basic_line):
# Splitting into segments close to half the total length
half_length = 157237.40665500844 # Approximate 1/2 total length of line
result = basic_line.split(half_length)

# Check if the segments are created
assert len(result) == 2

# Validate the first segment
assert result[0].vertices[0] == basic_line.vertices[0]
assert result[0].vertices[-1].longitude != basic_line.vertices[-1].longitude # Shouldn't reach the end

# Validate the second segment starts where the first ended
assert result[1].vertices[0] == result[0].vertices[-1]
assert result[1].vertices[-1] == basic_line.vertices[-1]


def test_split_with_remainder(basic_line):
# Distance doesn't evenly divide the total length
result = basic_line.split(5000)
assert len(result) > 1
assert all(len(segment.vertices) > 1 for segment in result)


def test_split_with_time_intervals(timed_line):
# Ensure time intervals are proportional
result = timed_line.split(157249)
assert len(result) == 2
assert result[0].dt.start == timed_line.dt.start
assert result[1].dt.end == timed_line.dt.end
assert result[0].dt.end == result[1].dt.start
assert result[0].dt.end < result[1].dt.end
assert result[0].dt.start < result[1].dt.start


def test_warning_on_large_distance(basic_line):
# Warning if the split distance exceeds the total length
with pytest.warns(UserWarning):
result = basic_line.split(400000)
assert len(result) == 1


def test_split_exact_division(basic_line):
# Distance perfectly divides the line
result = basic_line.split(157249 * 2)
assert len(result) == 1