专栏名称: 锐多宝
遥感技术教程、资讯与前沿论文
目录
相关文章推荐
疯狂区块链  ·  如何才能不返贫? ·  昨天  
疯狂区块链  ·  谁还在相信人性? ·  3 天前  
51好读  ›  专栏  ›  锐多宝

批量对齐栅格数据代码

锐多宝  · 公众号  ·  · 2024-11-15 23:22

正文

批量对齐栅格数据代码

最近经常遇到把一个栅格数据对齐到另外一个栅格数据。

经常会遇到同一个区域的几个栅格数据的分辨率、行列号和投影都不一致。

那可以用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






请到「今天看啥」查看全文