python活用gdal库进行批量重投影、重采样、裁剪

自然地理的数据经常保存为nc文件，因为nc文件包含经纬度、时间三个轴的基本信息。另一种保存方式是保存为tif文件，但tif文件只有经纬度信息，丢失了时间信息，折衷方案是对每个tif文件进行包含时间信息的命名，例如著名的GIMMS NDVI数据"PKU_GIMMS_NDVI_V1.2_19820101.tif"。

这会导致一个问题，就是一个日期的数据都保存为一个tif，如果要处理1982-2020年总共39年的NDVI数据，就需要循环很多次（39*24次）。arcgis的模型构建器或arcgis自带的python功能可以解决这个问题，反正循环就对了，但是我对arcgis的低效处理速度深恶痛绝，因而结合其他人的博客探索出了基于python的gdal库进行批量重投影、重采样、裁剪的高效方案。

具体来说，我不希望一个日期一个tif，那我就把所有数据组合成一个tif文件（参考我的前一篇博客：Python进行nc文件、tif文件、ndarray数据的相互转换-优快云博客）。我希望能快速重采样、裁剪、重投影，那我就用gdal库的warp功能。最终效果就是一次性对所有数据进行重采样/裁剪/重投影。下面是我分别对GRACE卫星数据进行裁剪、重采样的最后结果，运行速度10s以内。

1. 批量重投影

下面代码是针对tif文件的重投影，你可能会问哪里批处理了，好问题。只要我把所有数据保存为一个tif，再运行下面的代码就是批处理了

def gdal_reprojection(input_tif_path, output_tif_path, reference_tif_path=None):

    input_tif = gdal.Open(input_tif_path)

    if reference_tif_path==None:
        dstSRS = 'EPSG:4326'
    else:
        reference_tif = gdal.Open(reference_tif_path)
        dstSRS = reference_tif.GetProjection()

    gdal.Warp(output_tif_path,
              input_tif,
              dstSRS=dstSRS,  # 目标投影 (WGS 84)
              resampleAlg=gdal.GRA_NearestNeighbour,  # 重采样算法 (可选择其他，如 GRA_NearestNeighbour)
              dstNodata=-9999)  # 设置 NoData 值（可选）

    print("重投影完成，输出文件路径：", output_tif_path)

    return

2. 批量重采样

input_file_path是多个tif合并后的单个tif文件输入路径，reference_file_path是最后输出的参考tif文件，会修改原文件的投影方式、地理转换参数、经纬度范围（就是arcgis的重采样功能），重采样的模式可修改输入参数method。提供两种重采样的方案，gdal.ReprojectImage()和gdal.Warp()，两种都可以，亲测结果一致，用户只需使用的话，准备好输入文件就行，不用深究两种方案的差异。

def gdal_resample(input_file_path, reference_file_path, output_root=None, output_file_name=None, method=gdal.gdalconst.GRA_Bilinear):

    # 获取待重采样影像信息
    input_ras_file = gdal.Open(input_file_path, gdal.GA_ReadOnly)
    input_proj = input_ras_file.GetProjection()
    input_RasterCount = input_ras_file.RasterCount
    input_nodata = input_ras_file.GetRasterBand(1).GetNoDataValue()

    # 获取参考影像信息
    reference_file = gdal.Open(reference_file_path, gdal.GA_ReadOnly)
    reference_file_proj = reference_file.GetProjection()
    reference_file_trans = reference_file.GetGeoTransform()
    Width = reference_file.RasterXSize
    Height = reference_file.RasterYSize
    reference_lons = np.array([reference_file_trans[0] + i * reference_file_trans[1] for i in range(Width)])  # 经度
    reference_lats = np.array([reference_file_trans[3] + i * reference_file_trans[5] for i in range(Height)])  # 纬度
    
    # 创建输出文件路径
    if output_root is None:
        output_root = os.path.split(input_file_path)[0]
    input_file_name = os.path.basename(input_file_path)
    if output_file_name is None:
        output_file_name = os.path.splitext(input_file_name)[0] + "_resample.tif"
    output_file_path = os.path.join(output_root, output_file_name)

    ##### 重采样 方案1
    driver = gdal.GetDriverByName('GTiff')
    ds_output = driver.Create(output_file_path, Width, Height, input_RasterCount, gdal.GDT_Float32)
    ds_output.SetGeoTransform(reference_file_trans)
    ds_output.SetProjection(reference_file_proj)
    ds_output.GetRasterBand(1).SetNoDataValue(input_nodata)
    gdal.ReprojectImage(input_ras_file, ds_output, input_proj, reference_file_proj, method, 0.0, 0.0, )
    # 上一行的输入变量含义：输入数据集、输出数据集、输入投影、参考投影、重采样方法(最邻近内插\双线性内插\三次卷积等)、回调函数

    # gdalconst.GRA_NearestNeighbour: 'near',
    # gdalconst.GRA_Bilinear: 'bilinear',
    # gdalconst.GRA_Cubic: 'cubic',
    # gdalconst.GRA_CubicSpline: 'cubicspline',
    # gdalconst.GRA_Lanczos: 'lanczos',
    # gdalconst.GRA_Average: 'average',
    # gdalconst.GRA_RMS: 'rms',
    # gdalconst.GRA_Mode: 'mode',
    # gdalconst.GRA_Max: 'max',
    # gdalconst.GRA_Min: 'min',
    # gdalconst.GRA_Med: 'med',
    # gdalconst.GRA_Q1: 'q1',
    # gdalconst.GRA_Q3: 'q3',
    # gdalconst.GRA_Sum: 'sum',

    ##### 重采样 方案2
    # options = gdal.WarpOptions(
    #     outputBounds=(reference_lons.min(), reference_lats.min()+reference_file_trans[5], reference_lons.max()+reference_file_trans[1], reference_lats.max()),
    #     xRes=reference_file_trans[1],
    #     yRes=reference_file_trans[5],
    #     srcSRS=input_proj,
    #     dstSRS=reference_file_proj,
    #     resampleAlg=method,
    #     format='GTiff',
    #     dstNodata=input_nodata,
    #     srcNodata=input_nodata
    #     )
    # gdal.Warp(output_file_path, input_ras_file, options=options)

    print(output_file_path, ' 重采样完成')
    return

3. 批量裁剪

其实第2部分的代码也提供了裁剪的功能，与之不同的是，第2部分的裁剪是根据参考tif文件的网格信息进行裁剪的，第3部分的代码是根据shapefile文件裁剪的，分别对应于arcgis中使用矢量文件和栅格文件裁剪的功能。

def gdal_warp(input_file_path, shp_file_path=None, output_root=None, output_file_name=None):

    # 获取输出影像信息
    input_ras_file = gdal.Open(input_file_path, gdal.GA_ReadOnly)
    input_nodata = input_ras_file.GetRasterBand(1).GetNoDataValue()

    # 创建输出文件路径
    if output_root is None:
        output_root = os.path.split(input_file_path)[0]
    input_file_name = os.path.basename(input_file_path)
    if output_file_name is None:
        output_file_name = os.path.splitext(input_file_name)[0] + "_warp.tif"
    output_file_path = os.path.join(output_root, output_file_name)

    # 输入边界shp文件
    if shp_file_path is None:
        shp_file_path = r'F:\Arcgis\Boundary\Shp boundary\result.shp'

    # 裁剪
    gdal.Warp(output_file_path, input_ras_file,
                format='GTiff',
                cutlineDSName=shp_file_path,
                cropToCutline=True,
                srcNodata=input_nodata,
                dstNodata=input_nodata,)

    print(output_file_path, ' 裁剪完成')
    return