如何从rasterio rasterize生成的NDArray快速提取沿线像素值？

一、Numpy 高效索引替代循环 #

你要的替代低效循环的方法就是Numpy高级索引，直接用数组对数组进行索引，完全是向量化操作，效率比Python循环高几个数量级：

v = a[rows, cols]

针对你的示例，运行后直接得到：

array([71, 82], dtype=uint8)

这种方式底层基于C实现，没有Python循环的额外开销，处理数十亿点的场景也能大幅提速。

二、针对视线判断场景的优化 #

你的核心需求是判断线是否与建筑物相交（即线上是否存在建筑物像素），不需要提取所有像素值，因此可以进一步优化：
假设你栅格化时用1标记建筑物，0标记空白区域，那么只需判断线上是否存在1即可：

# 检查线上是否有建筑物像素
has_obstacle = numpy.any(a[rows, cols] == 1)

numpy.any()会短路求值（找到第一个满足条件的元素就停止），比提取所有值再判断更快。

三、大规模数据处理的进阶建议 #

面对单城市超10万组工业区对、数百万条线的场景，仅靠索引优化还不够，可结合以下方案：

• 批量处理：将多条线的rows和cols数组合并成大数组，一次性提取所有像素值，再按线的边界拆分结果，减少重复的IO和索引计算。
• 先做范围过滤：先计算每条线的外接矩形，如果矩形完全不与建筑物栅格的有效区域重叠，直接标记为"无遮挡"，跳过像素检查。
• 窗口读取栅格：用rasterio的窗口（Window）功能，只读取线对应的栅格区域，避免加载整个大栅格到内存，降低内存占用。
• 并行计算：工业区对之间相互独立，可使用多进程/多线程框架（如concurrent.futures）并行处理不同的工业区对，充分利用CPU资源。

关于Rasterio原生方法 #

Rasterio没有直接提取线上所有像素值的原生方法，但通过rasterio.transform.rowcol()转换坐标得到索引后，结合Numpy的向量化索引，就是最高效的实现方式，不需要额外的Rasterio API。

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