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

【问题标题】：Find length between two points on different linestrings查找不同线串上两点之间的长度
【发布时间】：2020-01-24 00:21:13
【问题描述】：

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

理想情况下，算法应该是这样的：

查找LineStrings 点属于哪个
查找路线的距离

数据几何：

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)

图形数据看起来像

我应该怎么做才能找到底点和顶点之间的距离？它们位于不同的LineStrings 上，它们之间还有其他LineStrings。

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

但在这种情况下，我不知道。

【问题讨论】：

标签： python distance shapely

【解决方案1】：

理想情况下，如果所有LineStrings 都有共同的端点，您可以使用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) 这个距离不是以米或公里为单位的。在某些步骤中，您必须从纬度/经度坐标进行转换。不幸的是，我不知道该怎么做。你可以自己弄清楚，或者我会研究它并在稍后编辑它:)

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

【讨论】：

shapely.ops.transform 可用于从一个坐标参考系转换到另一个 docs 。我将'EPSG:4326' 转换为'EPSG:32633'，它得到的结果与我使用OSM 测量的结果几乎相同。问题是：这种转换是否正确？如果正确，哪种 EPSG 更适合？
@EfYak 不幸的是，我不熟悉地理坐标系统转换。您可以尝试在GIS Stack Exchange 上查找或提问。
请问你是用什么来绘制这些图表的？
@EfYak 没什么特别的，真的。只是普通的交互式 Matplotlib 绘图来获得特写镜头和 MS Paint 将所有绘图放在一起:)

【解决方案2】：

从提供的数据开始。计算如下：

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

【讨论】：

我认为 OP 想要计算沿轨道的距离，而不仅仅是两点之间的距离。
@swatchai，谢谢你的回答，但这不是我想要的。当然，我没有透明地说明问题。我的意思是沿线的点之间的路线距离，而不是实际的大地距离。很抱歉。
另外，我更喜欢使用，例如，关系函数 inside() 来检查一个点是否属于一条线，然后获取坐标，而不是像您在示例中那样手动执行此操作。
@EfYak 我只向您展示相关步骤，不包括细节。您是否为回答者提供了the minimum runnable code？对于我向您展示的距离计算，您可以适应迭代沿 LINESTRING 沿轨道的 2 个点之间运行的所有顶点，以获得您想要的“长度”。