基于Python GDAL库实现图像的几何校正详细教程

Posted 2023-04-06 倾城一少

If the implementation is hard to explain, it’s a bad idea.
If the implementation is easy to explain, it may be a good idea.
——The Zen of python

---

最近需要利用Python的GDAL库对遥感图像进行几何校正，在网上搜了搜相关资料，大部分是来自李民录老师的《GDAL源码剖析与开发指南》及其博客的C++代码，关于Python的资料较少，于是便四处参考查阅，最终实现了校正功能，现总结整理一下，如果有任何意见建议，欢迎批评指正。

一、几何校正方法

  图像校正本质是建立一种从原始图像行列号到某种投影的数学关系，即实现图像行列坐标到投影坐标的转换。不同的校正方法利用了不同的方法来表示转换关系，但本质上式相同的。常用的几何校正方法包括：几何多项式校正、有理函数模型校正、局部区域校正模型、地理查找表校正等。

  GDAL库中可以实现的校正方法就包括以上四种方法，即：1~3次的几何多项式校正、RPC（有理函数系数）校正、TPS（薄板样条）校正、GeoLoc校正。

二、转换关系的描述

  不同的校正方法需要的信息也不相同，通常我们采用地面控制点（GCPs）的方式来建立转换关系，如果是RPC校正，则需要RPC文件来提供RPC系数。有的卫星数据，例如MODIS是包含GeoLocation Arrays提供每个像素的经纬度，例如Himawari-8卫星数据则直接提供了投影和地理变换参数(仿射变换六参数，这个最简单)。

三、Python中的GDAL几何校正

  Python中的几何校正主要靠 gdal.Warp() 函数来实现的，下面说一下主要流程：

3.1 读取未经校正的图像

  主要利用 gdal.Open()：

from osgeo import gdal
from osgeo import osr

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

3.2 构造地面控制点

  构造控制点列表 gcps_list：

# 实际控制点肯定要多的多，这里只写了四个点.
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.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校正利用上述流程即可校正，对于Geoloc校正，可以参考李老师博文中的思路，分为两种情况：

• 对于自带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信息，假设现在有一幅未校正图像 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校正，我没有数据，没有测试过。。。但是经过一番搜索[,看了里面gdal的单元测试文件，有如下思路以供参考:

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'])

RPC模型

RPC模型

• 1. 简介
• 2. RPC有理函数模型
• 3. 优秀的RPC相关知识博文
• • 3.1 RPC和RPB文件相关信息
  • 3.2 Matlab实现RPC正算反算
  • 3.3 使用Gdal进行批量的影像RPC正射校正
  • 3.4 （Python）卫星RPC有理多项式模型读取与正反投影坐标计算原理与实现
  • 3.5 高分影像批处理第三回——RPC文件与几何校正
• 4. 参考

1. 简介

有理多项式系数(rational polynomial coefficients,RPC)，实质是有理函数模型（Rational Function Model-RFM）。它可建立起像点和空间坐标之间的关系，不需要内外方位元素，回避成像的几何过程，可以广泛用于线阵影像处理中。RFM将像点坐标表示为以相应地面点空间坐标为自变量的多项式的比值。

卫星成像期间卫星的姿态控制导致影像的严格几何模型（所谓严格几何模型是指基于传统共线方程的严格几何模型）形式极其复杂，要利用其提取地球空间三维信息，需要在向用户提供影像的同时把卫星详细的轨道星历、传感器成像参数、成像方式等信息一并交付，并且最终用户需要具有摄影测量的专业知识和复杂的应用处理系统。为了降低对用户专业水平的需求，扩大用户范围，同时保护卫星的核心技术参数不被泄露，RPC定位模型应运而生。

RPC是一种与传感器无关的通用型成像几何模型。RPC是传感器严格几何模型的拟合形式，这里的严格几何模型是指通过平台载荷测量的平台运行轨迹参数﹑姿态参数﹑传感器安装参数及传感器内部几何参数等构建的像一地关系几何模型。由于这些参数不可避免地存在不同性质的误差，其拟合模型RPC也就存在着相应的误差。校正 RPC误差的传统方法是对地面点通过RPC投射到像方的像点进行一个多项式纠正，使投射像点坐标与测量像点坐标相吻合，从而达到消除误差的目的。

2. RPC有理函数模型

RPC 模型利用多项式的比值建立探测器像点坐标 d(line, sample)（可理解为行列号）与其对应的地面成像点坐标 D(latitude, lontitude, height) 的关系。RPC 也是各种传感器几何模型的一种通用表达形式，RPC 的正算形式为
$$\begin{array}{rl}& F_1=L_n=\frac{NumL(U,V,W)}{\mathrm{Den}L(U,V,W)}\\ & F_2=S_n=\frac{NumS(U,V,W)}{\mathrm{DenS}(U,V,W)}\end{array}$$
式中: $\mathrm{Num}L(U,V,W)=a_1+a_{2}V+a_{3}U+a_{4}W+a_{5}VU+a_{6}VW+a_{7}UW+a_8{V}^2+a_9{U}^2+a_{10}{W}^2+a_{11}UVW+a_{12}{V}^3+a_{13}V{U}^2+a_{14}V{W}^2+a_{15}{V}^{2}U+a_{16}{U}^3+a_{17}U{W}^2+a_{18}{V}^{2}W+a_{19}{U}^{2}W+a_{20}{W}^3;$
$$\begin{array}{c}\mathrm{Den}L(U,V,W)=b_1+b_{2}V+,\cdots ,+b_{19}{U}^{2}W+b_{20}{W}^3\\ \mathrm{Num}S(U,V,W)=c_1+c_{2}V+,\cdots ,+c_{19}{U}^{2}W+c_{20}{W}^3\\ \mathrm{DenS}(U,V,W)=d_1+d_{2}V+,\cdots ,+d_{19}{U}^{2}W+d_{20}{W}^3\end{array}$$

在式中, a 1 , a 2 ,