无人机航测数据处理太慢?Python脚本批量正射影像裁剪与重采样

引子：外业飞得爽，内业跑到哭

兄弟们，咱们干航测的，是不是都有这么个痛点？外业飞飞机那是风驰电掣，一天几百亩甚至几平方公里跟玩儿似的。可一回到内业，面对那动辄几百兆甚至几个G的正射影像（DOM），还有那一堆分幅要求，头都大了。

我就见过不少新手，甚至干了几年的人，还在用ArcGIS或者Global Mapper一张一张去裁剪，或者在那手动重采样。你要是只处理个十来张图也就罢了，要是碰到那种几百个图幅的测区，鼠标点到你手抽筋，进度条跑得比你女朋友翻脸还慢。这效率，甲方催图的时候你想死的心都有。

其实，这种重复性的体力活，本来就是代码的强项。今天老哥我就把这多年压箱底的“批量处理”思路掏出来，教你用Python写个脚本，把正射影像的裁剪和重采样搞定，让你的电脑替你加班。

工具与方法：GDAL与Rasterio的左右互搏

要在Python里搞定栅格数据，有两个库你绕不开：GDAL 和 Rasterio。

GDAL是GIS界的“祖师爷”，底层极其强大，但原生的Python接口写起来稍微有点繁琐，代码风格比较老派。而Rasterio则是后来之秀，它底层调用的还是GDAL，但是把接口封装得非常人性化，读起来就像在读普通的文件一样顺畅。

咱们今天的思路很简单：用Rasterio来读取和写入，用GDAL做底层支撑，写个循环，让电脑不知疲倦地替你干活。

脚本流程：先裁剪，后重采样，一气呵成

咱们假设你手里已经有了拼接好的正射影像（big_dom.tif）和一个存放裁剪框（Shapfile格式）的文件夹。

第一步：准备工作与裁剪

裁剪这事儿，本质上就是根据矢量范围，把栅格数据里对应的位置“抠”出来。

import osimport rasteriofrom rasterio.mask import maskimport fiona# 定义输入输出路径input_raster = 'big_dom.tif'output_folder = 'cut_results/'# 读取裁剪矢量框with fiona.open('grid_shp/grid_001.shp', "r") as shapefile:    shapes = [feature["geometry"] for feature in shapefile]# 开始裁剪with rasterio.open(input_raster) as src:    out_image, out_transform = mask(src, shapes, crop=True)    out_meta = src.meta.copy()# 更新元数据，这步很关键，不然出来的图坐标可能不对    out_meta.update({"driver": "GTiff","height": out_image.shape[1],"width": out_image.shape[2],"transform": out_transform})# 保存结果with rasterio.open(os.path.join(output_folder, 'cut_001.tif'), "w", **out_meta) as dest:        dest.write(out_image)

你看，这就完了。不用你在软件里点“导出”、选路径、起名字。你只要把这些代码封装成一个函数，配合一个for循环遍历你的分幅Shp文件夹，几百张图也就是喝杯茶的功夫。

第二步：重采样，把分辨率降下来

有时候甲方要求不一样，比如原始分辨率是2cm，但他想要10cm的数据方便查看。要是手动改，又得重新跑一遍重采样工具。来，咱们直接上代码：

from rasterio.warp import calculate_default_transform, reproject, Resamplingdef resample_raster(input_path, output_path, resampling_factor=2):with rasterio.open(input_path) as src:# 计算重采样后的变换参数，这里我们把分辨率乘以系数# 注意：这里假设是简单的倍数关系，实际项目中可根据目标分辨率精确计算        transform, width, height = calculate_default_transform(            src.crs, src.crs, src.width, src.height, *src.bounds,            resolution=(src.res[0] * resampling_factor, src.res[1] * resampling_factor)        )# 设置新文件的元数据        kwargs = src.meta.copy()        kwargs.update({'transform': transform,'width': width,'height': height        })# 执行重采样操作with rasterio.open(output_path, 'w', **kwargs) as dst:for i in range(1, src.count + 1):                reproject(                    source=rasterio.band(src, i),                    destination=rasterio.band(dst, i),                    src_transform=src.transform,                    src_crs=src.crs,                    dst_transform=transform,                    dst_crs=src.crs,                    resampling=Resampling.average)  # 这里用了平均值重采样，适合正射影

这里面有个坑我得提一句：Resampling.average（平均值）通常比 Resampling.nearest（最近邻）更适合正射影像，因为正射影像是照片，最近邻会让锯齿感很强，而平均值会让图像看起来更平滑自然。

示例数据与输出效果

我之前有个项目，大概10平方公里的测区，原始DOM大概8个G，需要裁剪成200多个标准图幅。

如果是手动操作，我在ArcGIS里设置好模型，跑一次大概要40分钟，还得时刻盯着有没有报错，电脑卡得连网页都刷不开。

用了Python脚本配合批处理，同样的电脑配置，跑完这200张图只用了不到12分钟。而且，脚本会在后台静默运行，你甚至可以切出去回个邮件，或者刷会儿手机。输出的Tif文件整整齐齐地躺在文件夹里，坐标信息、波段数一个不少，这就是自动化的魅力。

注意事项与调优：老手的避坑指南

虽然脚本好用，但有几点经验之谈得嘱咐大家：

1. 内存溢出问题：如果你的原始DOM特别大（比如几十个G），一次性读取可能会报内存错误。这时候就要用到Rasterio的“窗口读取”功能，分块处理，不要硬刚。
2. 坐标系一致性：一定要检查你的裁剪框Shp和原始影像是不是同一个坐标系！虽然Rasterio和GDAL支持动态投影，但精度多少会有点损失。作为测绘人，我建议统一转换好再做处理，保证精度万无一失。
3. 文件命名：批量处理最怕文件名冲突。写脚本的时候，记得把输出文件名设计好，比如用图幅号+分辨率命名，免得后面一锅粥。

写在最后

航测内业虽苦，但咱们不能傻苦。咱们干测绘的，其实玩的就是数据，学会用Python这把“快刀”，能帮你砍掉90%的重复劳动。

今天讲的这两个片段，拼起来就是一个完整的批量预处理流程。不要只看不练，找个小点的影像（几百兆那种），试着写个循环跑一遍。

将这份脚本应用到你的项目中，记得分享你的优化前后对比。 让咱们内业兄弟们的发际线，能多保住几根是一根。