Find length between two points on different linestrings

Question

I have a set of railways as LineString objects and railway halts as Point objects. I would like to find the distance along the track between two halts that are within two different rail segments.

Ideally, the algorithm should be like:

find to which LineStrings points belong
find the distance of the route

The data geometry:

LINESTRING (39.6578321 46.9014975, 39.6570969 46.9001473, 39.6559263 46.8989562, 39.6547141 46.8979608, 39.6334559 46.8814632, 39.631423 46.8798763)
LINESTRING (39.5993901 46.8638073, 39.5908007 46.8593408, 39.5779403 46.8525964, 39.5768336 46.8520311, 39.5751525 46.8511401, 39.5746447 46.8508623, 39.5738707 46.8504168, 39.5730826 46.849873, 39.5724067 46.8493054, 39.5719514 46.8488622, 39.5714961 46.848352, 39.571218 46.8480119, 39.5707563 46.8471837, 39.5705042 46.8467756, 39.5704733 46.8467123, 39.5701367 46.8460454, 39.5699599 46.8454946, 39.5697689 46.8447907, 39.5696772 46.844289, 39.5695147 46.8433865, 39.5690113 46.8372209, 39.5684811 46.830653, 39.5678923 46.8235858, 39.567313 46.8164827, 39.5664922 46.8059292, 39.5658968 46.7983574, 39.5648453 46.7849665, 39.5644591 46.7808009, 39.5641158 46.7759223, 39.5616257 46.7440695, 39.5614031 46.7412215, 39.5612064 46.7386438, 39.5609742 46.7360506, 39.5604696 46.7299663, 39.5604382 46.7295303, 39.5603741 46.728642, 39.5603751 46.7277074, 39.5604071 46.7263769, 39.5604906 46.7251804, 39.5734779 46.676293, 39.573848 46.6748994, 39.5790219 46.6550449)
LINESTRING (39.6312641 46.8797614, 39.6304511 46.8792256, 39.625757 46.8768124, 39.6199837 46.8739702)
LINESTRING (39.6198448 46.8739025, 39.614424 46.8712966, 39.6034935 46.8657577, 39.5995486 46.8638936)
POINT (39.6547141 46.8979608)
POINT (39.5609742 46.7360506)

Graphically data looks like

What should I do to find the distance between the bottom point and the top point? They lie on different LineStrings and there are other LineStrings between them.

If the lied on the same LineString I would simply use shapely.split with each point and then subtract resulted line.length from the total length.

But in a case like this one I have no idea.

score 1 · Accepted Answer · answered Jan 23 '20 at 14:08

Ideally, if all the LineStrings had common endpoints, you could just use linemerge to join them in one LineString and find the distance along the track as usual, but in your case, the lines do not touch each other:

In this case, I propose to look for the closest LineStrings pairs and join them together, for example, by using the following code:

from shapely.geometry import LineString
from shapely.wkt import loads


def merge_nontouching(lines):
    """Makes a LineString out of a list of LineString's by joining the closest pairs"""
    line, *candidates = lines
    while True:
        if not candidates:
            return line
        closest_candidate = min(candidates, key=line.distance)
        line = join_two(line, closest_candidate)
        candidates.pop(candidates.index(closest_candidate))


def join_two(line, other):
    """Joins two nontouching LineString's in one"""
    start, end = line.boundary
    if start.distance(other) > end.distance(other):
        return LineString([*line.coords, *other.coords])
    return LineString([*line.coords[::-1], *other.coords])


l1 = loads('LINESTRING (39.6578321 46.9014975, 39.6570969 46.9001473, 39.6559263 46.8989562, 39.6547141 46.8979608, 39.6334559 46.8814632, 39.631423 46.8798763)')
l2 = loads('LINESTRING (39.5993901 46.8638073, 39.5908007 46.8593408, 39.5779403 46.8525964, 39.5768336 46.8520311, 39.5751525 46.8511401, 39.5746447 46.8508623, 39.5738707 46.8504168, 39.5730826 46.849873, 39.5724067 46.8493054, 39.5719514 46.8488622, 39.5714961 46.848352, 39.571218 46.8480119, 39.5707563 46.8471837, 39.5705042 46.8467756, 39.5704733 46.8467123, 39.5701367 46.8460454, 39.5699599 46.8454946, 39.5697689 46.8447907, 39.5696772 46.844289, 39.5695147 46.8433865, 39.5690113 46.8372209, 39.5684811 46.830653, 39.5678923 46.8235858, 39.567313 46.8164827, 39.5664922 46.8059292, 39.5658968 46.7983574, 39.5648453 46.7849665, 39.5644591 46.7808009, 39.5641158 46.7759223, 39.5616257 46.7440695, 39.5614031 46.7412215, 39.5612064 46.7386438, 39.5609742 46.7360506, 39.5604696 46.7299663, 39.5604382 46.7295303, 39.5603741 46.728642, 39.5603751 46.7277074, 39.5604071 46.7263769, 39.5604906 46.7251804, 39.5734779 46.676293, 39.573848 46.6748994, 39.5790219 46.6550449)')
l3 = loads('LINESTRING (39.6312641 46.8797614, 39.6304511 46.8792256, 39.625757 46.8768124, 39.6199837 46.8739702)')
l4 = loads('LINESTRING (39.6198448 46.8739025, 39.614424 46.8712966, 39.6034935 46.8657577, 39.5995486 46.8638936)')
merged = merge_nontouching([l1, l2, l3, l4])
print(merged.wkt)
# LINESTRING (39.6578321 46.9014975, 39.6570969 46.9001473, 39.6559263 46.8989562, 39.6547141 46.8979608, 39.6334559 46.8814632, 39.631423 46.8798763, 39.6312641 46.8797614, 39.6304511 46.8792256, 39.625757 46.8768124, 39.6199837 46.8739702, 39.6198448 46.8739025, 39.614424 46.8712966, 39.6034935 46.8657577, 39.5995486 46.8638936, 39.5993901 46.8638073, 39.5908007 46.8593408, 39.5779403 46.8525964, 39.5768336 46.8520311, 39.5751525 46.8511401, 39.5746447 46.8508623, 39.5738707 46.8504168, 39.5730826 46.849873, 39.5724067 46.8493054, 39.5719514 46.8488622, 39.5714961 46.848352, 39.571218 46.8480119, 39.5707563 46.8471837, 39.5705042 46.8467756, 39.5704733 46.8467123, 39.5701367 46.8460454, 39.5699599 46.8454946, 39.5697689 46.8447907, 39.5696772 46.844289, 39.5695147 46.8433865, 39.5690113 46.8372209, 39.5684811 46.830653, 39.5678923 46.8235858, 39.567313 46.8164827, 39.5664922 46.8059292, 39.5658968 46.7983574, 39.5648453 46.7849665, 39.5644591 46.7808009, 39.5641158 46.7759223, 39.5616257 46.7440695, 39.5614031 46.7412215, 39.5612064 46.7386438, 39.5609742 46.7360506, 39.5604696 46.7299663, 39.5604382 46.7295303, 39.5603741 46.728642, 39.5603751 46.7277074, 39.5604071 46.7263769, 39.5604906 46.7251804, 39.5734779 46.676293, 39.573848 46.6748994, 39.5790219 46.6550449)

From here you can simply use the project function to get the distance:

p1 = loads('POINT (39.6547141 46.8979608)')
p2 = loads('POINT (39.5609742 46.7360506)')
distance_along_track = abs(merged.project(p1) - merged.project(p2))
print(distance_along_track)
# 0.21041206903426987

Just two notes:
1) This distance is not in meters or kilometers. On some step, you will have to do the conversion from the lat/lon coordinates. Unfortunately, I'm not aware of how to do it. You can figure it out by yourself, or I'll look into it and edit it in later :)

2) This may be not the most efficient approach if you have thousands of segments, but for the given example it will suffice.

`shapely.ops.transform` could be applied to transform from one coordinate reference system to another [docs](https://shapely.readthedocs.io/en/latest/manual.html#shapely.ops.transform) . I transformed `'EPSG:4326'` to `'EPSG:32633'` and it got me almost the same result as I measured using OSM. The question is: is this transformation correct and If so which EPSG fits better? — EfYak, Jan 23 '20 at 16:29
@EfYak Unfortunately, I'm not familiar with geographical coordinates systems transformations. You might try looking or asking a question on [GIS Stack Exchange](https://gis.stackexchange.com/). — Georgy, Jan 23 '20 at 16:44
@EfYak Nothing fancy, really. Just ordinary interactive Matplotlib plot to get close-ups and MS Paint to put all the plots together :) — Georgy, Jan 24 '20 at 08:57

score -1 · Answer 2 · answered Jan 23 '20 at 13:08

Starting from the provided data. The computation goes like this:

import shapely.wkt as wkt
import pyproj

wkt1 = 'LINESTRING (39.6578321 46.9014975, 39.6570969 46.9001473, 39.6559263 46.8989562, 39.6547141 46.8979608, 39.6334559 46.8814632, 39.631423 46.8798763)'
wkt2 = 'LINESTRING (39.5993901 46.8638073, 39.5908007 46.8593408, 39.5779403 46.8525964, 39.5768336 46.8520311, 39.5751525 46.8511401, 39.5746447 46.8508623, 39.5738707 46.8504168, 39.5730826 46.849873, 39.5724067 46.8493054, 39.5719514 46.8488622, 39.5714961 46.848352, 39.571218 46.8480119, 39.5707563 46.8471837, 39.5705042 46.8467756, 39.5704733 46.8467123, 39.5701367 46.8460454, 39.5699599 46.8454946, 39.5697689 46.8447907, 39.5696772 46.844289, 39.5695147 46.8433865, 39.5690113 46.8372209, 39.5684811 46.830653, 39.5678923 46.8235858, 39.567313 46.8164827, 39.5664922 46.8059292, 39.5658968 46.7983574, 39.5648453 46.7849665, 39.5644591 46.7808009, 39.5641158 46.7759223, 39.5616257 46.7440695, 39.5614031 46.7412215, 39.5612064 46.7386438, 39.5609742 46.7360506, 39.5604696 46.7299663, 39.5604382 46.7295303, 39.5603741 46.728642, 39.5603751 46.7277074, 39.5604071 46.7263769, 39.5604906 46.7251804, 39.5734779 46.676293, 39.573848 46.6748994, 39.5790219 46.6550449)'

# Create line objects from WKT code
line_1 = wkt.loads( wkt1 )
line_2 = wkt.loads( wkt2 )

# Get coordinates from specific vertices of lines
lon0, lat0 = line_1.coords[3]
lon1, lat1 = line_2.coords[32]

# Prep geodetic ref for computation
geod = pyproj.Geod(ellps='WGS84')

# Compute geodetic inverse problem
forward, back, dist = geod.inv(lon0, lat0, lon1, lat1)
print('Forward bearing: {}\nBackward bearing: {}\nDistance (m): {}'.format(forward, back, dist))

The output:

Forward bearing: -158.29035951647245
Backward bearing: 21.64128791894413
Distance (m): 19368.65144050848

I think OP wanted to calculate the distance *along the track*, not just the distance between two points. — Georgy, Jan 23 '20 at 13:36
@swatchai, thank you for you answer but it is not quite what I was looking for. Surely, I did not state the problem transparently. I meant the distance of the route between to points along lines not the actual geodetic distance. Sorry for that. — EfYak, Jan 23 '20 at 13:58
Also, I prefer to use, for example, relationship function within() to check if a point belongs to a line and then take coords rather than manually perform this operation as You did in your example. — EfYak, Jan 23 '20 at 14:05
@EfYak I only show you the relevant steps, excluding the details. Did you provide us `the minimum runnable code` for the answerers? For the computation of distance I show you, you can adapt to iterate through all the vertices along the LINESTRING that runs between the 2 points along the tracks to get the 'length' you want. — swatchai, Jan 23 '20 at 15:30

Find length between two points on different linestrings

2 Answers2