本節將從原理和代碼兩個方面講解遙感圖像的幾何校正。
原理
首先介紹幾何校正的概念：在遙感成像過程中，傳感器生成的圖像像元相對於地面目標物的實際位置發生了擠壓、扭曲、拉伸和偏移等問題，這一現象叫做幾何畸變。幾何畸變會給遙感圖像的定量分析、變化檢測、圖像融合、地圖測量或更新等處理帶來的很大誤差，所以需要針對圖像的幾何畸變進行校正，即幾何校正。
幾何校正分為幾何粗校正和幾何精校正。粗校正是利用空間位置變化關系，采用計算公式和輔助參數進行的校正，叫做系統幾何校正；精校正是在此基礎上，使圖像的幾何位置符合某種地理坐標系統，與地圖配准，調整亮度值，即利用地面控制點（GCP）做的幾何精校正。
幾何校正步驟：1.空間位置的變換（像元坐標）2.像元灰度值的重新計算，即重采樣。
1. 坐標變換
坐標變換分為直接法和間接法。
1）直接法：從原始圖像陣列出發，依次計算每個像元在輸出圖像中的坐標。直接法輸出的像元值大小不會發生變化，但輸出圖像中的像元分布不均勻。
2）間接法：從輸出圖像陣列出發，依次計算每個像元在原始圖像中的位置，然后計算原始圖像在該位置的像元值，再將計算的像元值賦予輸出圖像像元。此方法保證校正后的圖像的像元在空間上均勻分布，但需要進行灰度重采樣。該方法是最常用的幾何校正方法。
由上圖可見，直接法直接以原始圖像的坐標為基准點，坐標偏移到校正后的圖像，坐標的位置有很多出現在了像元的中間位置，所以直接輸出像元值大小導致像元分布不均勻。
而對於間接法。以輸出圖像的坐標為基准點，已經定義在了格點的位置上，此時反算出該點在原始圖像上對應的圖像坐標，坐標多數落在像元的中間位置。這里采用最鄰近法、雙線性內插和三次卷積法來計算該點的灰度值，達成重采樣的目的。
2. 重采樣
圖像數據經過坐標變換之后，像元中心的位置發生改變，其在原始圖像的位置不一定是整數行\列，需要根據輸出圖像各像元在原始圖像中對應的位置，對原始圖像重采樣，建立新的柵格矩陣。
主要有最近鄰法、雙線性內插法、三次卷積法。原理這里不再敘述。其中最近鄰法沒有改變圖像光譜信息，能保證定量分析的准確性，但是圖像不夠平滑。后兩者雖然數據連續平滑了，但是光譜信息遭到了破壞，比較適合於制圖，不適用於定量分析。
代碼實現
在熟悉了幾何校正的原理后，可以動手實現這一環節。這里主要介紹三種方法，具體如下。
1.首先介紹最基礎也是最底層的方法。(這里是分步切片得到偏移后的坐標，然后重采樣）
from osgeo import gdal def get_indices(source_ds, target_width, target_height): source_geotransform = source_ds.GetGeoTransform() source_width = source_ds.GetGeoTransform[1] source_height = source_ds.GetGeoTransform[5] dx = target_width/source_width dy = target_height/source_height target_x = np.arange(dx/2, source_ds.RasterXSize, dx) target_y = np.arange(dy/2, source_ds.RasterYSize, dy) return np.meshgrid(target_x, target_y) def bilinear(in_data, x, y): x-=0.5 y-=0.5 x0=np.floor(x).astype(int) x1=x0+1 y0=np.floor(y).astype(int) y1=y0+1 ul=in_data[y0,x0]*(y1-y)*(x1-x) ur=in_data[y0,x1]*(y1-y)*(x-x0) ll=in_data[y1,x0]*(y-y0)*(x1-x) lr=in_data[y1,x1]*(y-y0)*(x-x0) return ul+ur+ll+lr if __name__ =="__main__": in_fn=r'' out_fn=r'' cell_size=(0.02,-0.02) in_ds = gdal.Open(in_fn) x, y=get_indices(in_ds,*cell_size) outdata=bilinear(in_ds.ReadAsArray(), x, y) driver=gdal.GetDriverByName('GTiff') rows,cols=outdata.shape out_ds=driver.Create(out_fn,cols,rows,1,gdal.GDT_ Int32) gt=list(in_ds.GetGeoTransform()) gt[1]=cell_size[0] gt[5]=cell_size[1] out_ds.SetGeoTransform(gt) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(outdata) out_band.FlushCache() out_band.ComputeStatistics(False)
2.第二種方法是ReadAsArray()函數，設置buf_xsize, buf_ysize這個緩沖區大小，就可以按這個大小重采樣（最鄰近法）。同時需要設置新的gt[1]\gt[5]以及行列數，這里代碼省略，有興趣查閱之前的gdal筆記。
3.第三種方法是gdalwarp或gdal.ReprojectImage。這里有重采樣的部分，具體如下。
from osgeo import gdal, gdalconst # 1.設置分辨率投影 # inputfile = r'E:\桌面文件保存路徑\GIS\test\quickbird_test.dat' # outputfile = r'E:\桌面文件保存路徑\GIS\test\resampling_test.tif' # # inputds = gdal.Open(inputfile, gdal.GA_ReadOnly) # 待重采樣圖像 # band1 = inputds.GetRasterBand(1) # xsize = inputds.RasterXSize # ysize = inputds.RasterYSize # nbands = inputds.RasterCount # inputProj = inputds.GetProjection() # inputTrans = inputds.GetGeoTransform() # # # new_x = int(xsize * (inputTrans[1] / 10)) # 計算新行列數 # new_y = int(ysize * (inputTrans[5] / -10)) # 注意抵消負號 # # print(inputTrans[1]) # print(inputTrans[5]) # # inputTrans = list(inputTrans) # inputTrans[1] = 10 # 分辨率 # inputTrans[5] = -10 # 注意設置為-，否則無法SetGeoTransform # # print(inputTrans) # # driver = gdal.GetDriverByName('GTiff') # out_ds = driver.Create(outputfile, new_x, new_y, nbands, band1.DataType) # out_ds.SetGeoTransform(inputTrans) # out_ds.SetProjection(inputProj) # # gdal.ReprojectImage(inputds, out_ds, inputProj, out_ds.GetProjection(), gdalconst.GRA_Bilinear) # # del out_ds # 2.從圖1重采樣到圖2 inputfile = r'E:\桌面文件保存路徑\IDL\imgs\TMSPOT\TM-30m.img' #高光譜 referencefile = r'E:\桌面文件保存路徑\IDL\imgs\TMSPOT\bldr_sp.img' #高分辨率 outputfile = r'E:\test.tif' #重采樣輸出圖像 inputras = gdal.Open(inputfile, gdal.GA_ReadOnly) #待重采樣圖像，高光譜低分辨率 inputProj = inputras.GetProjection() inputTrans = inputras.GetGeoTransform() nbands = inputras.RasterCount reference = gdal.Open(referencefile, gdal.GA_ReadOnly) #目標圖像，高分辨率全色波段 referenceProj = reference.GetProjection() referenceTrans = reference.GetGeoTransform() bandreference = reference.GetRasterBand(1) x = reference.RasterXSize y = reference.RasterYSize driver = gdal.GetDriverByName('GTiff') output = driver.Create(outputfile, x, y, nbands, bandreference.DataType) output.SetGeoTransform(referenceTrans) output.SetProjection(referenceProj) gdal.ReprojectImage(inputras, output, inputProj, referenceProj, gdalconst.GRA_Bilinear) del output
4.此外完整的幾何校正和重采樣還可以用opencv的一些api。
from matplotlib import pyplot as plt import cv2 import numpy as np img = cv2.imread('aier.jpg') rows,cols = img.shape[:2] pts1 = np.float32([[50,50],[200,50],[50,200]]) #對應校正點 pts2 = np.float32([[10,100],[200,20],[100,250]]) M = cv2.getAffineTransform(pts1,pts2) //計算仿射變換矩陣 #第三個參數：變換后的圖像大小 res = cv2.warpAffine(img,M,(rows,cols)) //根據矩陣和圖像重采樣 plt.subplot(121) plt.imshow(img[:,:,::-1]) plt.subplot(122) plt.imshow(res[:,:,::-1]) #畫出圖片，但不顯示 plt.show() #顯示
需要注明，本文中原理部分引用了朱文泉老師著書的內容，部分代碼摘自csdn博客。