优化Python代码：Numpy数组过滤性能提升方案咨询

首先得提一句：你的原代码逻辑存在准确性问题——当前的筛选方式是分别检查经度是否接近线段的某个端点经度、纬度是否接近线段的某个端点纬度，这和“点到LineString的实际距离小于阈值”不是同一个条件。比如一个点的经度接近线段端点的经度，但纬度相差很大，实际距离远超过阈值，但你的代码会误判它符合条件。所以我们先修正逻辑，再优化性能。

针对140万规模的坐标数组，这里有两种高效的解决方案：

方案1：向量化计算点到线段的精确距离（推荐，精确且高效） #

直接用numpy向量化计算每个点到线段的最短距离，时间复杂度为O(n)，比原方法的O(n*m)快很多，同时结果完全准确。

from shapely.geometry import LineString
import numpy as np
import time

def point_to_segment_distance(points, seg_start, seg_end):
    """
    计算批量点到线段的最短距离
    :param points: (N, 2) 数组，每行是一个点的[经度, 纬度]
    :param seg_start: (2,) 数组，线段起点的[经度, 纬度]
    :param seg_end: (2,) 数组，线段终点的[经度, 纬度]
    :return: (N,) 数组，每个点到线段的距离
    """
    # 计算线段向量和点到起点的向量
    seg_vec = seg_end - seg_start
    point_vec = points - seg_start
    
    # 计算点在线段上的投影比例，超出线段范围则取端点
    t = np.dot(point_vec, seg_vec) / np.dot(seg_vec, seg_vec)
    t = np.clip(t, 0.0, 1.0)
    
    # 计算最近点坐标并求距离
    closest_points = seg_start + t[:, np.newaxis] * seg_vec
    return np.linalg.norm(points - closest_points, axis=1)

# 初始化数据（模拟大规模场景）
HLat22 = np.random.rand(1413917) * 100
HLong22 = np.random.rand(1413917) * (-200) + 100
polygon2 = LineString([(-86.79191, 32.47155), (-86.786797, 32.47005)])
thresh = 0.005

start_time = time.time()
# 构造点数组（经度在前，纬度在后，对应LineString的坐标顺序）
points = np.column_stack((HLong22, HLat22))
seg_start = np.array(polygon2.coords[0])
seg_end = np.array(polygon2.coords[1])

# 计算所有点到线段的距离并生成布尔索引
distances = point_to_segment_distance(points, seg_start, seg_end)
indices = distances < thresh

print(f"耗时：{time.time() - start_time:.4f} 秒")
print(f"符合条件的点数量：{np.sum(indices)}")

这个方法一步完成所有计算，没有中间过滤步骤，对于140万数据，耗时通常在0.1秒以内，远快于原方法。

方案2：使用KDTree快速查询（适合复杂LineString） #

如果你的LineString是由多段线段组成的复杂路径，可以先在路径上采样足够多的点，然后用KDTree查询每个点到采样点的最近距离，近似判断是否在阈值内。这种方法适合复杂路径，速度也很可观：

from shapely.geometry import LineString
from scipy.spatial import KDTree
import numpy as np
import time

def sample_linestring(linestring, num_samples=200):
    """
    在LineString上采样均匀分布的点
    :param linestring: 目标LineString对象
    :param num_samples: 采样点数量
    :return: (num_samples, 2) 数组，采样点的[经度, 纬度]
    """
    # 计算线段总长度的归一化距离
    coords = np.array(linestring.coords)
    segment_lengths = np.linalg.norm(coords[1:] - coords[:-1], axis=1)
    cumulative_lengths = np.cumsum(segment_lengths)
    normalized_distances = np.insert(cumulative_lengths, 0, 0) / cumulative_lengths[-1]
    
    # 生成采样点
    t_values = np.linspace(0, 1, num_samples)
    sampled_points = np.array([
        linestring.interpolate(t, normalized=True) for t in t_values
    ])
    return np.column_stack((sampled_points.x, sampled_points.y))

# 初始化数据（模拟大规模场景）
HLat22 = np.random.rand(1413917) * 100
HLong22 = np.random.rand(1413917) * (-200) + 100
polygon2 = LineString([(-86.79191, 32.47155), (-86.786797, 32.47005)])
thresh = 0.005

start_time = time.time()
# 采样线段上的点并构建KDTree
sampled_points = sample_linestring(polygon2, num_samples=200)
kdtree = KDTree(sampled_points)

# 查询所有点的最近距离
points = np.column_stack((HLong22, HLat22))
distances, _ = kdtree.query(points, k=1)
indices = distances < thresh

print(f"耗时：{time.time() - start_time:.4f} 秒")
print(f"符合条件的点数量：{np.sum(indices)}")

注意：采样点数量越多，结果越精确，但耗时会略有增加，一般200个采样点对于单线段来说完全足够。

为什么原代码的np.in1d这么慢？ #

原代码中np.in1d(HLong22, The_X) & np.in1d(HLat22, The_Y)的时间复杂度是O(N*M)（N是140万，M是1.5万），这相当于210亿次比较，自然很慢。而上面的向量化方法是O(N)，KDTree方法是O(N log M)，都远快于原方法。

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