查找不同线串上两点之间的长度

Find length between two points on different linestrings

我有一组铁路作为 LineString 对象和铁路停靠点作为 Point 对象。我想找到两个不同铁路段内的两个停靠站之间沿着轨道的距离。

理想情况下,算法应该是这样的:

  1. 找到 LineStrings 个点属于
  2. 找到路线的距离

数据几何:

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)

图形数据看起来像

我应该怎么做才能找到最低点和最高点之间的距离?它们位于不同的 LineString 上,它们之间还有其他 LineString

如果在同一个 LineString 上撒谎,我会简单地对每个点使用 shapely.split,然后从总长度中减去结果 line.length

但是我不知道像这种情况。

从提供的数据开始。计算过程如下:

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

输出:

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

理想情况下,如果所有 LineString 都有共同的端点,您可以使用 linemerge 将它们合并为一个 LineString 并照常计算沿轨道的距离,但是在你的情况下,线条不会相互接触:

在这种情况下,我建议寻找最接近的 LineStrings 对并将它们连接在一起,例如,使用以下代码:

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)

从这里您可以简单地使用 project 函数来获取距离:

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

只有两个注释:
1) 此距离不是以米或公里为单位。在某些步骤中,您将不得不从 lat/lon 坐标进行转换。不幸的是,我不知道该怎么做。你可以自己想办法,或者我会研究一下,稍后再编辑:)

2) 如果您有数千个分段,这可能不是最有效的方法,但对于给定的示例,它就足够了。