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校正。[1]
二、转换关系的描述
不同的校正方法需要的信息也不相同,通常我们采用地面控制点(GCPs)的方式来建立转换关系,如果是RPC校正,则需要RPC文件来提供RPC系数。有的卫星数据,例如MODIS是包含GeoLocation Arrays提供每个像素的经纬度,例如Himawari-8卫星数据则直接提供了投影和地理变换参数(仿射变换六参数[2],这个最简单)。
三、Python中的GDAL几何校正
Python中的几何校正主要靠 gdal.Warp()
函数来实现的,下面说一下主要流程:
1.读取未经校正的图像
主要利用 gdal.Open()
:
1 2 3 4
| from osgeo import gdal from osgeo import osr
dataset = gdal.Open(r'xxx.tif', gdal.GA_Update)
|
2.构造地面控制点
构造控制点列表 gcps_list
:
1 2 3 4 5
| 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()
的说明[3]:
1 2 3 4
| gdal.GCP([x], [y], [z], [pixel], [line], [info], [id])
|
3.添加地面控制点到图像
在添加之前需要指定一个空间参考,这里以WGS84基准的地理坐标系(经纬度)为例:
1 2 3 4
| sr = osr.SpatialReference() sr.SetWellKnownGeogCS('WGS84')
dataset.SetGCPs(gcps, sr.ExportToWkt())
|
4.进行校正
校正就直接调用gdal.Warp()
就可以了:
1 2 3 4
| 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()
的参数说明[4]:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
| gdal.Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)
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):
|
多项式校正和TPS校正利用上述流程即可校正,对于Geoloc校正,可以参考李老师博文[5] [6]中的思路,分为两种情况:
- 对于自带GeoLocation元数据段的MODIS数据,利用
gdal.Info()
查看波段信息,直接利用gdal.Warp()
设置geoloc=True
后,对目标波段进行校正即可:
1 2 3 4
| 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 文件,内容如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
| <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校正,我没有数据,没有测试过。。。但是经过一番搜索[8] ,看了里面gdal的单元测试文件,有如下思路以供参考:
1 2 3 4 5 6 7 8 9
| 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'])
|
四、总结
人生苦短,我用Python。
参考资料
[1]:李民录. GDAL源码剖析与开发指南[M]. 人民邮电出版社, 2014.
[2]:https://blog.csdn.net/ivan_ljf/article/details/9226463 关于GDAL计算图像坐标的几个问题
[3]:Garrad C. Geoprocessing with Python[M]. 2016.
[4]:gdal.py Python gdal 包中定义代码
[5]:https://blog.csdn.net/liminlu0314/article/details/9181959 使用GDAL工具对FY3系列卫星数据进行校正
[6]:https://blog.csdn.net/liminlu0314/article/details/8623607 使用GDAL对HDF数据进行校正
[7]:https://stackoverflow.com/questions/50472624/reference-lat-lon-vrt-files-in-data-vrt-for-gdaltranslaste Reference Lat/Lon vrt files in data vrt for gdal translaste
[8]:https://svn.osgeo.org/gdal/ 这里面有gdal的单元测试文件