numpy.memmap 是加载超大栅格的轻量方案，通过内存映射按需读取磁盘数据；需配合 rasterio 获取元数据，严格匹配 dtype 和 shape，并始终同步 transform 与 crs 以保地理语义完整。
用 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：
立即学习“Python免费学习笔记（深入）”；
# 假设 ndvi 是 shape=(h,w) 的 float32 数组 mask = np.where(ndvi > 0.6, 1, 0) # 直接生成整型掩膜 # 或更省内存：用布尔索引原地修改 ndvi[ndvi <= 0.6] = 0 ndvi[ndvi > 0.6] = 1
注意点：
仅使用 Python 将 Claw Code 虚拟化框架与 OpenClaw 集成，提供 184 种工具和 200 多条命令，用于安全分析和代码质量...
-
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）。结果数值对，但业务解释全错。
安全做法：
- 所有邻域操作前，先用
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)
常见疏漏：
- 忘了
transform—— 输出文件无地理定位，QGIS 里显示在赤道原点 -
crs写成字符串如'EPSG:4326'，但 rasterio 严格要求CRS对象或 proj4 字符串（带+proj=） - 多波段写入时，
count和dst.write(arr, band_idx)的band_idx不匹配，导致波段错位
地理空间栅格的“高效”不在单次计算快慢，而在整个流程不丢失坐标语义——从加载、运算到落盘，每一步的 transform 和 crs 都得像影子一样跟着数据走，少一个环节，后面就得人工对齐，成本远高于前期多写两行代码。
