测绘人新玩具:用Python和Rasterio处理海量GeoTIFF文件

测绘人新玩具：用Python和Rasterio处理海量GeoTIFF文件

做遥感和GIS的朋友，肯定被海量栅格数据折磨过。

一个项目几十上百个GeoTIFF，每个文件几百兆甚至几个G。打开慢、处理慢、导出更慢。以前用传统软件处理，等得人都快睡着了。不对，准确地说，是等到怀疑人生。

后来接触了GDAL和Rasterio，才发现处理大型影像可以这么快。

GDAL是底层库，功能强大但API比较底层。Rasterio是Python封装，用起来友好多了。读取、切片、重采样、导出，几行代码就能搞定。

读取影像数据，最基本的操作。

import rasteriowith rasterio.open('image.tif') as src:print(f'尺寸: {src.width} x {src.height}')print(f'波段数: {src.count}')print(f'坐标系: {src.crs}')    data = src.read(1)

切片处理也很方便。影像太大一次读不进来？按块读就行。

import rasterioimport numpy as npwith rasterio.open('large_image.tif') as src:    window = rasterio.windows.Window(0, 0, 1000, 1000)    data = src.read(1, window=window)    red = src.read(3, window=window)    nir = src.read(4, window=window)    ndvi = (nir - red) / (nir + red + 1e-10)

重采样也是常见需求。分辨率不统一？统一一下。

from rasterio.warp import reproject, Resamplingimport numpy as npwith rasterio.open('image.tif') as src:    dst_data = np.zeros((500, 500), dtype=src.dtypes[0])    reproject(        source=rasterio.band(src, 1),        destination=dst_data,        src_transform=src.transform,        src_crs=src.crs,        dst_transform=rasterio.transform.from_bounds(0, 0, 500, 500, 500, 500),        dst_crs=src.crs,        resampling=Resampling.bilinear    )

从原始数据到分析产出，整个工作流可以连起来。

import rasterioimport numpy as npfrom rasterio.warp import reproject, Resamplingwith rasterio.open('raw.tif') as src:    dst_data = np.zeros((1000, 1000), dtype=src.dtypes[0])    reproject(        source=rasterio.band(src, 1),        destination=dst_data,        src_transform=src.transform,        src_crs=src.crs,        dst_transform=rasterio.transform.from_bounds(0, 0, 1000, 1000, 1000, 1000),        dst_crs=src.crs,        resampling=Resampling.bilinear    )    profile = src.profile.copy()    profile.update(height=1000, width=1000, count=1)with rasterio.open('processed.tif', 'w', **profile) as dst:        dst.write(dst_data, 1)

性能这块儿，有几个小技巧。

并行处理能省不少时间。可以用Python的multiprocessing，多个文件同时处理。特别是处理几个G的影像时，分块读取避免内存用尽。还有，用rasterio的windowed读写，别一次性把整个文件塞进内存。

我上个月处理一个省域的DOM数据，几十个文件加起来上百个G。用分块加并行的方式，原来要跑一天的任务，两个小时就搞定了。

用这个工作流处理你的GeoTIFF，提升分析效率。