ST_交叉口行为

Question

我有一个线串，它可以追踪其自身路径的一部分（一条往返轨迹）。当我尝试将此线串与多边形相交以尝试获取给定多边形内行驶的距离时，相交操作仅保留旅程的重复部分的一条腿。有没有办法保留两条轨迹，以便准确测量行驶的距离。

就 MRE 而言，我想要

select 
    ST_Length(ST_Intersection(
        ST_GeomFromText('POLYGON((-1 -1, 1 -1, 1 1, -1 1, -1 -1))'),
        ST_GeomFromText('LINESTRING(0 -1, 0 0, 1 0, 0 0, 0 1)')))

给出结果 4 而不是 3。

我似乎能够做到这一点的唯一方法是将线串分成其组成部分并将每个部分与多边形相交。然而，这需要更多的计算时间才能完成，因为轨迹可能很长，有很多段，并且有大约 10,000 个多边形要相交。

我实际上是在 Python 中使用 shapely 来完成此操作，但我想我应该先与 PostGIS 检查一下，看看这种行为是否符合 OGC 标准（确实如此）。

Answer 1

我通过执行以下操作在 geopandas 中有效地实现了这一点

import geopandas as gpd
from shapely import LineString, Polygon

line = LineString([[0, -2],[0, 0],[2, 0],[0, 0],[0, 2]])

polygons = gpd.GeoDataFrame(geometry=[
    Polygon([[-1, -1], [-1, -3], [1, -3], [1, -1], [-1, -1]]),
    Polygon([[-1, 1], [-1, -1], [1, -1], [1, 1], [-1, 1]]),
    Polygon([[-1, 3], [-1, 1], [1, 1], [1, 3], [-1, 3]])
]).rename_geometry('polygon')
polygons['code'] = ['A', 'B', 'C']
polygons.set_index('code', inplace=True)

##################################################################

points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(*line.xy)).rename_geometry('point')
points['next_point'] = points.shift(-1)['point']
segments = gpd.GeoDataFrame(geometry = points['point'].shortest_line(points['next_point']).iloc[:-1:])

s = (segments
    .sjoin(polygons)
    .set_index('code')
    .intersection(polygons, align = True))

(s[s.notna()]
    .length
    .groupby(level=0).sum())