Python利⽤GDAL对图像进⾏⼏何校正
⼀、⼏何校正⽅法
在线翻译句子 图像校正本质是建⽴⼀种从原始图像⾏列号到某种投影的数学关系，即实现图像⾏列坐标到投影坐标的转换。不同的校正⽅法利⽤了不同的⽅法来表⽰转换关系，但本质上式相同的。常⽤的⼏何校正⽅法包括：⼏何多项式校正、有理函数模型校正、局部区域校正模型、地理查找表校正等。
GDAL库中可以实现的校正⽅法就包括以上四种⽅法，即：1~3次的⼏何多项式校正、RPC(有理函数系数)校正、TPS(薄板样条)校正、GeoLoc校正。
⼆、转换关系的描述
elevator 不同的校正⽅法需要的信息也不相同，通常我们采⽤地⾯控制点（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
datat = 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')
# 添加控制点
datat.SetGCPs(gcps, sr.ExportToWkt())研究生考试报名时间
4.进⾏校正
校正就直接调⽤gdal.Warp()就可以了：
# tps校正重采样：最邻近法
dst_ds = gdal.Warp(r'xxx_dst.tif', datat, 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 = Fal, width = 0, height = 0, srcSRS = None,
dstSRS = None, srcAlpha = Fal, dstAlpha = Fal, warpOptions = None,
挤压英文errorThreshold = None, warpMemoryLimit = None, creationOptions = None,
outputType = GDT_Unknown, workingType = GDT_Unknown, resampleAlg = None,
srcNodata = None, dstNodata = None, multithread = Fal, tps = Fal,
rpc = Fal, geoloc = Fal, polynomialOrder = None,
transformerOptions = None, cutlineDSName = None, cutlineLayer = None,
cutlineWhere = None, cutlineSQL = None, cutlineBlend = None,
ropToCutline = Fal, copyMetadata = True, metadataConflictValue = None,
tColorInterpretation = Fal, callback = None, callback_data = None):
# options --- 字符串数组, 字符串或者空值
# format --- 输出格式 ("GTiff", )
# outputBounds --- 结果在⽬标空间参考的边界范围(minX, minY, maxX, maxY)
# outputBoundsSRS --- 结果边界范围的空间参考, 如果在dstSRS中没有指定的话，采⽤此参数wsn是什么意思
# xRes, yRes --- 输出分辨率，即像素的⼤⼩
hook
# targetAlignedPixels --- 是否强制输出边界是输出分辨率的倍数
# width --- 输出栅格的列数
# height --- 输出栅格的⾏数
# srcSRS --- 输⼊数据集的空间参考
# dstSRS --- 输出数据集的空间参考
情人节英文# srcAlpha --- 是否将输⼊数据集的最后⼀个波段作为alpha波段
# dstAlpha --- 是否强制创建输出
# outputType --- 输出栅格的变量类型 (gdal.GDT_Byte, )
# workingType --- working type (gdal.GDT_Byte, )
# 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校正
usc# 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 --- 元数据冲突值
# tColorInterpretation --- 是否强制将输⼊栅格颜⾊表给输出栅格
# 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 ⽂件，内容如下：
口罩标价超千万<VRTDatat 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[ <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>
</VRTDatat>
其中关键的是<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('', ds, dstSRS='EPSG:4326', xRes=0.0002, yRes=0.0002,
rpc=True, transformerOptions = [r'RPC_DEM=X:\DEM.tif'])
