GDAL重投影重采样像元配准对齐

研究通常会涉及到多源数据，需要进行基于像元的运算，在此之前需要对数据进行地理配准、空间配准、重采样等操作。那么当不同来源，不同分辨率的数据重采样为同一空间分辨率之后，各个像元不一一对应，有偏移该怎么办呢？

在ArcGIS进行重采样操作时（resample 或者project raster）可以通过设置Environment --> Processing Extent --> Snap Raster 为参考栅格数据，解决这一问题。详见我的这一篇博客和知乎文章

但面对大批量数据的时候，我们希望通过编程解决这一问题，下面分享gdal中对这一问题的解决思路。

# -*- coding: utf-8 -*-

"""

Created on Wed Sep 23 23:32:25 2020

以栅格A参考，对栅格B进行重投影操作，输出结果和栅格A像元大小一致，相互对齐

https://gis.stackexchange.com/questions/234022/resampling-a-raster-from-python-without-using-gdalwarp

https://gis.stackexchange.com/questions/268119/how-can-2-geotiff-rasters-be-aligned-and-forced-to-have-same-resolution-in-pytho

@author: pan

"""

from osgeo import gdal

from osgeo import gdalconst

def pixel_geo_register(infile,outfile,reffile,methods):

    '''

    infile:输入文件

    outfile:输出文件

    reffile:参考文件

    methods:重采样方法

            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_Mode:mode

    '''

    # 打开tif文件

    in_ds  = gdal.Open(infile, gdalconst.GA_ReadOnly) # 输入文件

    ref_ds = gdal.Open(reffile, gdalconst.GA_ReadOnly) # 参考文件

    # 参考文件与输入文件的的地理仿射变换参数与投影信息

    in_trans = in_ds.GetGeoTransform()

    in_proj = in_ds.GetProjection()

    ref_trans = ref_ds.GetGeoTransform()

    ref_proj = ref_ds.GetProjection()

    # 参考文件的波段参考信息

    band_ref = ref_ds.GetRasterBand(1)

    # 输入文件的行列数

    x = ref_ds.RasterXSize

    y = ref_ds.RasterYSize

    # 创建输出文件

    driver= gdal.GetDriverByName('GTiff')

    output = driver.Create(outfile, x, y, 1, gdalconst.GDT_UInt16)

    # 设置输出文件地理仿射变换参数与投影

    output.SetGeoTransform(ref_trans)

    output.SetProjection(ref_proj)

    # 重投影，插值方法为双线性内插法

    gdal.ReprojectImage(in_ds, output, in_proj, ref_proj, methods)

    # 关闭数据集与driver

    in_ds = None

    ref_ds = None

    driver  = None

    output = None

if __name__ == '__main__':

    infile=r''

    outfile=r''

    reffile=r''

    methods=gdalconst.GRA_Bilinear

    pixel_geo_register(infile,outfile,reffile,methods)

GDAL重投影重采样像元配准对齐的相关教程结束。