Python 进行栅格数据处理的基本步骤

栅格数据是以网格形式组织的空间数据，常用于遥感影像、数字高程模型（DEM）等。Python 中常用的栅格数据处理库是 rasterio，它基于 GDAL 提供了高效的栅格数据读写和操作功能。

栅格数据处理的基本步骤

1. 读取栅格数据：使用 rasterio.open() 读取栅格文件，获取数据、元数据（如分辨率、坐标系等）。
2. 重采样：调整栅格数据的分辨率，使其与其他数据匹配。常用方法包括最近邻法、双线性插值法等。
3. 裁剪：根据感兴趣区域（AOI）裁剪栅格数据，通常使用矢量边界或掩膜。
4. 拼接：将多个栅格数据拼接成一个更大的栅格数据集。
5. 保存结果：将处理后的栅格数据保存为文件（如 GeoTIFF）。

示例代码注释说明

import rasterio  # 导入 rasterio 库，用于处理栅格数据
from rasterio.enums import Resampling  # 导入重采样方法
from rasterio import mask, merge  # 导入裁剪和拼接功能
import numpy as np  # 导入 numpy 库，用于数组操作
from shapely.geometry import Point  # 导入 Point 类，用于定义几何对象
from tqdm import tqdm  # 导入 tqdm 库，用于显示进度条
import os  # 导入 os 库，用于文件路径操作

def process_raster_data(input_raster1, input_raster2, output_path):
    """
    处理栅格数据：重采样、裁剪、拼接
    :param input_raster1: 输入栅格文件1的路径
    :param input_raster2: 输入栅格文件2的路径
    :param output_path: 输出文件保存路径
    """
    print("开始重采样...")
    resampled_file = f"{output_path}/resampled.tif"  # 重采样后的文件路径
    with rasterio.open(input_raster1) as src:  # 打开输入栅格文件1
        target_resolution = (30, 30)  # 目标分辨率（30米 x 30米）
        
        # 计算重采样比例
        scale_factor = src.res[0] / target_resolution[0]
        new_width = int(src.width * scale_factor)  # 新宽度
        new_height = int(src.height * scale_factor)  # 新高度
        
        # 读取数据并进行重采样
        data = src.read(
            out_shape=(src.count, new_height, new_width),  # 输出形状
            resampling=Resampling.bilinear  # 使用双线性插值法
        )
        
        # 计算新的变换矩阵
        transform = src.transform * src.transform.scale(
            (src.width / new_width),
            (src.height / new_height)
        )
        
        # 保存重采样后的数据
        with rasterio.open(
            resampled_file,
            'w',
            driver='GTiff',  # 文件格式为 GeoTIFF
            height=new_height,
            width=new_width,
            count=src.count,  # 波段数
            dtype=src.dtypes[0],  # 数据类型
            crs=src.crs,  # 坐标系
            transform=transform,  # 变换矩阵
        ) as dst:
            dst.write(data)  # 写入数据
    
    print("开始裁剪...")
    clipped_file = f"{output_path}/clipped.tif"  # 裁剪后的文件路径
    with rasterio.open(resampled_file) as src:  # 打开重采样后的文件
        height, width = src.height, src.width  # 获取高度和宽度
        center_y, center_x = height // 2, width // 2  # 计算中心点
        Y, X = np.ogrid[:height, :width]  # 创建网格
        radius = min(height, width) // 8  # 计算裁剪半径
        
        # 创建圆形掩膜
        mask = (X - center_x)**2 + (Y - center_y)**2 <= radius**2
        
        # 读取数据并应用掩膜
        data = src.read()
        for i in range(src.count):
            data[i][~mask] = src.nodata if src.nodata else 0  # 将掩膜外的值设为无数据
        
        # 保存裁剪后的数据
        with rasterio.open(
            clipped_file,
            'w',
            driver='GTiff',
            height=height,
            width=width,
            count=src.count,
            dtype=src.dtypes[0],
            crs=src.crs,
            transform=src.transform,
            nodata=src.nodata  # 无数据值
        ) as dst:
            dst.write(data)
    
    print("开始拼接...")
    mosaic_file = f"{output_path}/mosaic.tif"  # 拼接后的文件路径
    with rasterio.open(clipped_file) as src1, \  # 打开裁剪后的文件
         rasterio.open(input_raster2) as src2:  # 打开输入栅格文件2
        if src1.count != src2.count:  # 检查波段数是否匹配
            print(f"警告：输入文件的波段数不匹配（{src1.count} vs {src2.count}）")
        
        # 拼接两个栅格文件
        mosaic, out_transform = merge.merge([src1, src2])
        
        # 保存拼接后的数据
        with rasterio.open(
            mosaic_file,
            'w',
            driver='GTiff',
            height=mosaic.shape[1],
            width=mosaic.shape[2],
            count=mosaic.shape[0],
            dtype=src1.dtypes[0],
            crs=src1.crs,
            transform=out_transform,
            nodata=src1.nodata
        ) as dst:
            dst.write(mosaic)

if __name__ == "__main__":
    input_raster1 = r'<栅格文件1路径>'
    input_raster2 = r'<栅格文件2路径>'
    output_path = r'<输出文件保存路径>'
    
    os.makedirs(output_path, exist_ok=True)  # 创建输出目录（如果不存在）
    
    try:
        process_raster_data(input_raster1, input_raster2, output_path)  # 处理栅格数据
        print("处理完成！")
    except Exception as e:
        print(f"处理过程中出现错误: {str(e)}")  # 捕获并打印错误信息