如何判断多个点位于Shapely LineString的东侧或西侧

针对包含数百个点的场景，我们可以通过点到LineString的最近线段判断的方法高效解决，核心思路是找到每个点最接近的折线线段，再通过叉积计算点相对于该线段的左右侧，映射为东西侧。

实现步骤 #

1. 对每个点，找到LineString上距离它最近的线段；
2. 利用叉积计算点相对于该线段的位置（左/右）；
3. 结合坐标系将左右侧映射为东西侧（假设东为x轴正方向，北为y轴正方向）。

代码实现 #

基础版本（适用于线段数量不多的场景） #

from shapely import geometry
from shapely.geometry import Point

def point_rel_line_side(line, point):
    # 获取点在LineString上的投影距离
    proj_dist = line.project(point)
    current_dist = 0.0
    # 遍历线段找到投影所在的区间
    for i in range(len(line.coords) - 1):
        p1 = Point(line.coords[i])
        p2 = Point(line.coords[i+1])
        seg_len = p1.distance(p2)
        if current_dist <= proj_dist <= current_dist + seg_len:
            seg = geometry.LineString([p1, p2])
            break
        current_dist += seg_len
    else:
        # 投影在最后一段线段上
        seg = geometry.LineString([line.coords[-2], line.coords[-1]])
    
    # 计算叉积：(x2-x1)*(y-y1) - (y2-y1)*(x-x1)
    x1, y1 = seg.coords[0]
    x2, y2 = seg.coords[1]
    x, y = point.coords[0]
    cross = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1)
    
    # 根据叉积符号判断东西侧（考虑浮点误差）
    if cross > 1e-9:
        return "东侧"
    elif cross < -1e-9:
        return "西侧"
    else:
        return "在线段上"

# 示例数据
xy = geometry.LineString([
    (284013.3168553732, 4929336.870078533),
    (284012.39530582004, 4929324.339915148),
    (284017.5631161644, 4929250.361344838),
    (284022.5080302361, 4929182.366627739),
    (284030.5078563042, 4929110.832004678)
])

pts = [(280800,4894570),(280900,4894620),(295750,4894620)]

# 批量判断每个点的位置
for p in pts:
    point = Point(p)
    side = point_rel_line_side(xy, point)
    print(f"点{p}位于LineString的{side}")

优化版本（适用于线段/点数量较多的场景） #

通过预计算线段累积长度，用二分查找快速定位最近线段，将时间复杂度从O(NM)优化为O(NlogM)（N为点数，M为线段数）：

from shapely import geometry
from shapely.geometry import Point
import bisect

def precompute_segment_lengths(line):
    """预计算LineString各线段的累积长度"""
    seg_end_lengths = []
    total = 0.0
    for i in range(len(line.coords)-1):
        p1 = Point(line.coords[i])
        p2 = Point(line.coords[i+1])
        length = p1.distance(p2)
        total += length
        seg_end_lengths.append(total)
    return seg_end_lengths

def point_rel_line_side_fast(line, point, seg_end_lengths):
    proj_dist = line.project(point)
    # 二分查找定位投影所在的线段
    idx = bisect.bisect_right(seg_end_lengths, proj_dist)
    if idx >= len(seg_end_lengths):
        seg = geometry.LineString([line.coords[-2], line.coords[-1]])
    else:
        seg = geometry.LineString([line.coords[idx], line.coords[idx+1]])
    
    # 计算叉积判断位置
    x1, y1 = seg.coords[0]
    x2, y2 = seg.coords[1]
    x, y = point.coords[0]
    cross = (x2 - x1) * (y - y1) - (y2 - y1) * (x - x1)
    
    if cross > 1e-9:
        return "东侧"
    elif cross < -1e-9:
        return "西侧"
    else:
        return "在线段上"

# 预计算线段累积长度
seg_end_lengths = precompute_segment_lengths(xy)

# 批量判断
for p in pts:
    point = Point(p)
    side = point_rel_line_side_fast(xy, point, seg_end_lengths)
    print(f"点{p}位于LineString的{side}")

原理说明 #

• 最近线段定位：折线是连续的，点的相对位置由最近的线段决定，通过投影距离找到对应线段；
• 叉积判断逻辑：叉积的符号反映点相对于线段的左右侧（沿线段走向），结合平面坐标系（东为x+），将左侧映射为东侧，右侧映射为西侧；
• 浮点误差处理：使用小阈值（1e-9）避免浮点计算导致的误判。

内容的提问来源于stack exchange，提问作者Italo Lopes