批量对齐栅格数据代码
最近经常遇到把一个栅格数据对齐到另外一个栅格数据。
经常会遇到同一个区域的几个栅格数据的分辨率、行列号和投影都不一致。
那可以用gdal和rasterio来实现。代码如下：
import rasterio
from rasterio.warp import reproject, Resampling
import os
def align_raster(input_path, reference_path, output_path):
"""
将输入栅格对齐到参考栅格
"""
try:
# 读取参考栅格获取目标参数
with rasterio.open(reference_path) as src2:
dst_crs = src2.crs
dst_transform = src2.transform
dst_width = src2.width
dst_height = src2.height
# 对输入栅格进行重投影和重采样
with rasterio.open(input_path) as src1:
# 准备输出元数据
meta = src1.meta.copy()
meta.update({
'crs': dst_crs,
'transform': dst_transform,
'width': dst_width,
'height': dst_height
})
# 创建输出文件并进行重投影
with rasterio.open(output_path, 'w', **meta) as dst1:
for i in range(1, src1.count + 1):
reproject(
source=rasterio.band(src1, i),
destination=rasterio.band(dst1, i),
src_transform=src1.transform,
src_crs=src1.crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=Resampling.bilinear
)
print(f"已保存对齐后的栅格: {output_path}")
return True
except Exception as e:
print(f"处理文件 {input_path} 时发生错误: {str(e)}")
return False
def main():
# 定义输入输出路径
input_dir = r'你的文件夹路径'
output_dir = r'输出路径'
reference_path = r'参考栅格的tif'
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
# 获取所有tif文件
tif_files = [f for f in os.listdir(input_dir) if f.endswith('.tif')]
total_files = len(tif_files)
print(f"共找到 {total_files} 个TIF文件需要处理")
# 处理每个文件
successful = 0
failed = 0
for i, tif_file in enumerate(tif_files, 1):
input_path = os.path.join(input_dir, tif_file)
output_path = os.path.join(output_dir, tif_file)
print(f"\n处理第 {i}/{total_files} 个文件: {tif_file}")
if align_raster(input_path, reference_path, output_path):
successful += 1
else:
failed += 1
# 打印处理结果统计
print("\n处理完成:")
print(f"成功处理: {successful} 个文件")
print(f"处理失败: {failed} 个文件")
print(f"总计文件: {total_files} 个")
# 验证最后一个处理的文件
if successful > 0:
print("\n验证最后处理的文件:")
with rasterio.open(output_path) as aligned1, rasterio.open(reference_path) as src2:
print(f"参考栅格大小: {src2.width} x {src2.height}")
print(f"对齐后栅格大小: {aligned1.width} x {aligned1.height}")
print(f"参考栅格变换矩阵: {src2.transform}")
print(f"对齐后变换矩阵: {aligned1.transform}")
if (aligned1.transform == src2.transform and
aligned1.width == src2.width and
aligned1.height == src2.height):
print("栅格已成功对齐到参考栅格的网格和分辨率")
else:
print("对齐结果有误，请检查")
if __name__ == "__main__":
main()
环境：python 3.10 ，安装gdal和rasterio包。
原理是：
读取参考栅格，获取目标参数： 坐标系统(CRS) 变换矩阵(transform) 影像宽度和高度 对每个输入栅格文件： 读取源文件 更新元数据(metadata) 使用双线性重采样方法进行重投影 保存对齐后的文件
注：该代码由chatgpt o1一次性生成，生成后就可以直接运行，没有报错。感叹一下，太强了。
