如何在Python中通过NumPy高效处理地理空间栅格-Python教程-PHP中文网

numpy.memmap 是加载超大栅格的轻量方案，通过内存映射按需读取磁盘数据；需配合 rasterio 获取元数据，严格匹配 dtype 和 shape，并始终同步 transform 与 crs 以保地理语义完整。

如何在python中通过numpy高效处理地理空间栅格

用 `numpy.memmap` 加载超大栅格不崩内存

地理栅格动辄几GB，直接 np.load 或 rasterio.open().read() 全部读入内存，Python 很快 OOM。关键不是“怎么读”，而是“按需读哪块”。numpy.memmap 是最轻量的解法——它不把数据全载入 RAM，而是让 NumPy 把磁盘文件当数组“映射”进来，访问时才触发页加载。

实操要点：

先用 rasterio 读元数据（宽、高、dtype、transform），别碰 .read()
用 np.memmap(filename, dtype='float32', mode='r', shape=(height, width)) 创建映射对象；注意 dtype 必须和原始栅格一致，否则数值错乱
切片操作如 arr[1000:1500, 2000:2500] 会自动触发对应磁盘区域读取，速度取决于硬盘随机 I/O 性能
避免对 memmap 对象调 np.mean() 这类全量聚合——它会强制加载全部数据；改用分块迭代或 np.average() 配合 weights 参数做加权局部统计

用 `np.where` 和掩膜数组替代 Python 循环做像元级条件运算

写 for i in range(h): for j in range(w): if arr[i,j] > 100: ... 处理百万级像元？慢百倍且易出错。NumPy 的向量化逻辑才是地理栅格的正确打开方式。

典型场景：提取 NDVI > 0.6 的植被像元，并设为 1，其余为 0：

# 假设 ndvi 是 shape=(h,w) 的 float32 数组
mask = np.where(ndvi > 0.6, 1, 0)  # 直接生成整型掩膜
# 或更省内存：用布尔索引原地修改
ndvi[ndvi  0.6] = 1

注意点：

np.where 返回新数组，内存翻倍；高频操作优先选布尔索引原地修改
涉及多个波段组合（如计算 NDVI = (nir-red)/(nir+red)），务必用 np.divide(nir-red, nir+red, out=out_arr, where=(nir+red)!=0) 避免除零警告和 NaN 扩散
如果掩膜本身是 GeoTIFF（比如土地利用分类图），用 rasterio.mask.mask() 先裁剪再运算，别手动套坐标循环

用 `scipy.ndimage` 做地理邻域分析时小心坐标系错位

地理栅格的“上-下-左-右”对应的是行、列索引，但 scipy.ndimage 的滤波核（如 uniform_filter）默认按数组轴方向滑动——这没问题；真正容易翻车的是空间尺度混淆：你算的是 3×3 像元均值，但实际地理范围可能是 30m×30m（Landsat）或 10m×10m（Sentinel-2）。结果数值对，但业务解释全错。

Claw Code Suite (Python-Only Edition)

仅使用 Python 将 Claw Code 虚拟化框架与 OpenClaw 集成，提供 184 种工具和 200 多条命令，用于安全分析和代码质量...

下载

安全做法：

所有邻域操作前，先用 rasterio 拿到 transform，确认 transform.a（x 方向像素尺寸）和 transform.e（y 方向像素尺寸，通常为负）
若需固定地理尺度（如 1km×1km 窗口），先换算成像元数：kernel_size = int(1000 / abs(transform.a))，再传给 uniform_filter(arr, size=kernel_size)
避免用 gaussian_filter 做平滑——它假设各向同性，而地理坐标系中 x/y 尺寸常不同，会导致畸变；改用 uniform_filter 或自定义核

保存处理结果时必须同步写入 `transform` 和 `crs`

用 np.save 存数组？np.savetxt 导出 ASCII？这些操作丢掉所有地理参考信息，后续无法叠加、无法投影、GIS 软件打不开。真正的“保存”是写回 GeoTIFF。

核心步骤就三行：

with rasterio.open(
    'output.tif', 'w',
    driver='GTiff',
    height=arr.shape[0],
    width=arr.shape[1],
    count=1,
    dtype=arr.dtype,
    crs='+proj=epsg:4326',  # 或 rasterio.crs.CRS.from_epsg(4326)
    transform=original_transform  # 必须！从输入文件读出来
) as dst:
    dst.write(arr, 1)

常见疏漏：