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)