本节将从原理和代码两个方面讲解遥感图像的几何校正。
原理
首先介绍几何校正的概念：在遥感成像过程中，传感器生成的图像像元相对于地面目标物的实际位置发生了挤压、扭曲、拉伸和偏移等问题，这一现象叫做几何畸变。几何畸变会给遥感图像的定量分析、变化检测、图像融合、地图测量或更新等处理带来的很大误差，所以需要针对图像的几何畸变进行校正，即几何校正。
几何校正分为几何粗校正和几何精校正。粗校正是利用空间位置变化关系，采用计算公式和辅助参数进行的校正，叫做系统几何校正；精校正是在此基础上，使图像的几何位置符合某种地理坐标系统，与地图配准，调整亮度值，即利用地面控制点（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博客。