Python 利用GDAL對圖像進行幾何校正

　　原文鏈接：https://blog.csdn.net/qq_27045589/article/details/81062586

一、幾何校正方法

　　圖像校正本質是建立一種從原始圖像行列號到某種投影的數學關系，即實現圖像行列坐標到投影坐標的轉換。不同的校正方法利用了不同的方法來表示轉換關系，但本質上式相同的。常用的幾何校正方法包括：幾何多項式校正、有理函數模型校正、局部區域校正模型、地理查找表校正等。
  GDAL庫中可以實現的校正方法就包括以上四種方法，即：1~3次的幾何多項式校正、RPC(有理函數系數)校正、TPS(薄板樣條)校正、GeoLoc校正。

二、轉換關系的描述

　　不同的校正方法需要的信息也不相同，通常我們采用地面控制點（GCPs）的方式來建立轉換關系，如果是RPC校正，則需要RPC文件來提供RPC系數。有的衛星數據，例如MODIS是包含GeoLocation Arrays提供每個像素的經緯度，例如Himawari-8衛星數據則直接提供了投影和地理變換參數(仿射變換六參數。

三、Python中的GDAL幾何校正

　　Python中的幾何校正主要靠 gdal.Warp() 函數來實現的，下面說一下主要流程：

1.讀取未經校正的圖像

　　利用 gdal.Open()：

from osgeo import gdal
from osgeo import osr

dataset = gdal.Open(r'xxx.tif', gdal.GA_Update)

2.構造地面控制點

# 實際控制點肯定要多的多，這里只寫了四個點.做成人機交互更好
gcps_list = [gdal.GCP(-111.931075, 41.745836, 0, 1078, 648),
             gdal.GCP(-111.901655, 41.749269, 0, 3531, 295),
             gdal.GCP(-111.899180, 41.739882, 0, 3722, 1334),
             gdal.GCP(-111.930510, 41.728719, 0, 1102, 2548)]
    

　　附控制點構造函數gdal.GCP()的說明

gdal.GCP([x], [y], [z], [pixel], [line], [info], [id])
# x、y、z是控制點對應的投影坐標，默認為0;
# pixel、line是控制點在圖像上的列、行位置，默認為0;
# info、id是用於說明控制點的描述和點號的可選字符串，默認為空.

3.添加地面控制點到圖像

　　在添加之前需要指定一個空間參考，這里以WGS84基准的地理坐標系(經緯度)為例：

sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('WGS84')
# 添加控制點
dataset.SetGCPs(gcps, sr.ExportToWkt())

4.進行校正

　　校正就直接調用gdal.Warp()就可以了：

# tps校正 重采樣:最鄰近法
dst_ds = gdal.Warp(r'xxx_dst.tif', dataset, format='GTiff', tps=True, 
                   xRes=0.05, yRes=0.05, dstNodata=65535, srcNodata=65535, 
                   resampleAlg=gdal.GRIORA_NearestNeighbour， outputType=gdal.GDT_Int32)

　　附gdal.Warp()的參數說明：

gdal.Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)

# destNameOrDestDS --- 輸出數據集路徑或對象
# srcDSOrSrcDSTab --- 數據集對象或文件名or數據集對象或文件名的數組
# 關鍵字參數是gdal.WarpOptions()的返回值，或者直接定義gdal.WarpOptions()

gdal.WarpOptions(options = [], format = 'GTiff', outputBounds = None, 
                 outputBoundsSRS = one, xRes = None, yRes = None, 
                 targetAlignedPixels = False, width = 0, height = 0, srcSRS = None,
                 dstSRS = None, srcAlpha = False, dstAlpha = False, warpOptions = None,
                 errorThreshold = None, warpMemoryLimit = None, creationOptions = None,
                 outputType = GDT_Unknown, workingType = GDT_Unknown, resampleAlg = None,
                 srcNodata = None, dstNodata = None, multithread = False, tps = False, 
                 rpc = False, geoloc = False, polynomialOrder = None, 
                 transformerOptions = None, cutlineDSName = None, cutlineLayer = None,
                 cutlineWhere = None, cutlineSQL = None, cutlineBlend = None, 
                 ropToCutline = False, copyMetadata = True, metadataConflictValue = None,
                 setColorInterpretation = False, callback = None, callback_data = None):
# options --- 字符串數組, 字符串或者空值
# format --- 輸出格式 ("GTiff", etc...)
# outputBounds --- 結果在目標空間參考的邊界范圍(minX, minY, maxX, maxY)
# outputBoundsSRS --- 結果邊界范圍的空間參考, 如果在dstSRS中沒有指定的話，采用此參數
# xRes, yRes --- 輸出分辨率，即像素的大小
# targetAlignedPixels --- 是否強制輸出邊界是輸出分辨率的倍數
# width --- 輸出柵格的列數
# height --- 輸出柵格的行數
# srcSRS --- 輸入數據集的空間參考
# dstSRS --- 輸出數據集的空間參考
# srcAlpha --- 是否將輸入數據集的最后一個波段作為alpha波段
# dstAlpha --- 是否強制創建輸出
# outputType --- 輸出柵格的變量類型 (gdal.GDT_Byte, etc...)
# workingType --- working type (gdal.GDT_Byte, etc...)
# warpOptions --- list of warping options
# errorThreshold --- 近似轉換的誤差閾值(誤差像素數目)
# warpMemoryLimit --- 工作內存限制 Bytes
# resampleAlg --- 重采樣方法
# creationOptions --- list of creation options
# srcNodata --- 輸入柵格的無效值
# dstNodata --- 輸出柵格的無效值
# multithread --- 是否多線程和I/O操作
# tps --- 是否使用Thin Plate Spline校正方法
# rpc --- 是否使用RPC校正
# geoloc --- 是否使用地理查找表校正
# polynomialOrder --- 幾何多項式校正次數
# transformerOptions --- list of transformer options
# cutlineDSName --- cutline數據集名稱
# cutlineLayer --- cutline圖層名稱
# cutlineWhere --- cutline WHERE 子句
# cutlineSQL --- cutline SQL語句
# cutlineBlend --- cutline blend distance in pixels
# cropToCutline --- 是否使用切割線范圍作為輸出邊界
# copyMetadata --- 是否復制元數據
# metadataConflictValue --- 元數據沖突值
# setColorInterpretation --- 是否強制將輸入柵格顏色表給輸出柵格
# callback --- 回調函數
# callback_data --- 用於回調的用戶數據

# 參數很多，有些也沒有試驗過，有翻譯不對的地方也請批評指正。

　　多項式校正和TPS校正經過上述步驟即可實現，分為兩種情況：

• 對於自帶GeoLocation元數據段的MODIS數據，利用gdal.Info()查看波段信息，直接利用gdal.Warp()設置geoloc=True后,對目標波段進行校正即可：

ds = gdal.Warp(r'X:\ModisNearest.tif',
               r'HDF4_EOS:EOS_SWATH:"X:\MOD021KM.A2018130.0220.061.2018130132414\MOD021KM.A2018130.0220.061.2018130132414.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB',
               width=2030, height=1354, format='GTiff', geoloc=True,
               dstSRS=sr.ExportToWkt(), resampleAlg=gdal.GRIORA_NearestNeighbour)

• 對於其他類型數據，則需要構造VRT文件，然后指定geoloc信息 [7]，假設現在有一幅未校正圖像 XXX.tif，還有 longitude.tif，latitude.tif 兩個經緯度文件(寫成一個文件也可以，只不過需要改 X_BAND 和 Y_BAND 的值)，於是我們構造一個 xxx.vrt 文件，內容如下：

<VRTDataset rasterXSize="60" rasterYSize="39">
  <Metadata domain="GEOLOCATION">
    <MDI key="SRS">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]</MDI>
    <MDI key="X_DATASET">X:\longitude.tif</MDI>
    <MDI key="X_BAND">1</MDI>
    <MDI key="PIXEL_OFFSET">0</MDI>
    <MDI key="PIXEL_STEP">1</MDI>
    <MDI key="Y_DATASET">X:/latitude.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="LINE_OFFSET">0</MDI>
    <MDI key="LINE_STEP">1</MDI>
  </Metadata>
  <VRTRasterBand dataType="Int16" band="1">
    <ColorInterp>Gray</ColorInterp>
    <NoDataValue>65535</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">X:/XXX.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="100" ySize="100"/>
      <DstRect xOff="0" yOff="0" xSize="100" ySize="100"/>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

　　其中關鍵的是<Metadata>段中的9個key，其中X_DATASET和Y_DATASET指定了行列對應的經緯度波段，其他標簽含義從略，可看參考資料。VRT文件后半部分的<SourceFilename>指定了需要校正的文件。
　　有了VRT文件，我們就可以進行校正了，輸入改為vrt文件路徑，geoloc=True用Warp()校正。

RPC校正

rpc = [ "HEIGHT OFF=l09", "LINE NUM COEFF=-0. 001245683 -0. 09427649 -1. 006342 "
        "-1. 954469e-05 0. 001033926 2.020534e-08 -3.845472e-07 一0.002075817 "
        "0.0005520694 0 -4.642442e-06 -3.271793e-06 2. 705977e-05 -7.634384e-07 "
        "-2.132832e-05 -3.248862e-05 -8.17894e-06 -3.678094e-07 2.002032e-06 "
        "3.693162e-08", "LONG OFF=7.1477"
        "SAMP DEN COEFF=l " ......]
ds.SetMetadata(rpc,'RPC')
dst_ds = gdal.Warp('output.xxx', ds, dstSRS='EPSG:4326', xRes=0.0002, yRes=0.0002,
                    rpc=True, transformerOptions = [r'RPC_DEM=X:\DEM.tif'])