GDAL几何校正
GDAL几何校正之Geoloc校正
1. 几何校正原理
常见的几何校正方法有几何多项式校正、有理函数模型校正、局部区域校正模型、地理查找表校正。GDAL库中可以实现的校正方法有：几何多项式校正、RPC校正（有理函数模型）、TPS(薄板函数模型)校正、Geoloc校正。
2. Geoloc校正（地理查找表校正）
2.1 直接校正
对于自带Geolocation元数据(lon/lat)的数据，可以查看波段信息（如子数据集路径等），然后直接调用gdalwarp进行校正。
比如对modis的mod02数据，我们分以下两步进行几何校正。
(1) 查看数据集路径
from osgeo import gdal
datasets = gdal.Open(r"E:\Tablefile\modis\MOD021KM.A2020350.0325.061.2020350130838.hdf")
# 可以用这个来查找EV_1KM_RefSB数据集路径
# 获取hdf中的子数据集
SubDatasets = datasets.GetSubDatasets()
# 获取子数据集的个数
SubDatasetsNum = len(datasets.GetSubDatasets())
# 输出各子数据集的信息
print("子数据集一共有{0}个: ".format(SubDatasetsNum))
for i in range(SubDatasetsNum):
print(datasets.GetSubDatasets()[i])
(2) gdalwarp.exe
# mod02 --> ref（几何校正正确，与heg/MRT结果相同）
gdalwarp.exe -geoloc -t_srs EPSG:4326 "HDF4_EOS:EOS_SWATH:"E:\Tablefile\modis\MOD021KM.A2020350.0325.061.2020350130838.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB" E:\Tablefile\MOD021KM.tif
(3) 此外可以通过gdalinfo.exe查看数据集信息
gdalinfo.exe HDF4_EOS:EOS_SWATH:"E:\Tablefile\modis\MOD021KM.A2020350.0325.061.2020350130838.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB
但同时也发现了一点问题：
# mod03 --> sala/saza/sola/soza（几何校正错误！！！猜测是lon/lat问题？）
# 同样的方法对mod03地理定位文件进行几何校正，校正明显失败，与HEG的结果不同。
gdalwarp.exe -geoloc -t_srs EPSG:4326 "HDF4_EOS:EOS_SWATH:"E:\\Tablefile\\modis\\MOD03.A2020350.0325.061.2020350090352.hdf":MODIS_Swath_Type_GEO:SensorZenith" E:\Tablefile\MOD03_salz.tif
猜测是lon/lat不配对的原因，但是查看数据集后发现mod03(2030X1354)与其对应mod02(406X271)数据的lon/lat数值范围基本相同，因此可能不是这个原因。求解？
2.2 调用VRT文件校正
对于本身不带lon/lat的数据，需要构建VRT文件进行校正。
常用于风云4A、Himawari-8等卫星数据，这两种静止轨道卫星的校正方法有两种：
-
构建lon/lat数据集，进行地理查找表校正
-
根据行偏移、列偏移量，对不同分辨率数据的行列号和经纬度进行推导，以构建经纬度和值的对应关系，从而实现校正。（速度较慢）
介绍第一种方法
(以mod02为例，lon/lat读取mod数据集中的lon/lat子数据集)
1. 查询待几何校正数据的路径 以及 lon/lat 路径(同上)
(1) 如果该卫星提供了行列号与经纬度的对应数据，如.raw(FY-4A)，那么可以用np.fromfile进行读取，生成lon/lat的栅格文件。
(2) 如果没有提供这种数据，可以根据行列号建立经纬度数据集（行列号=>经纬度）。
2. 调用gdal_translate.exe，后接需要几何校正的数据集路径，生成VRT文件
gdal_translate.exe -of VRT HDF4_EOS:EOS_SWATH:"E:\\Tablefile\\modis\\MOD021KM.A2020350.0325.061.2020350130838.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB E:/Tablefile/back.vrt
3. 修改VRT文件的X_DATASET/Y_DATASET => lon/lat，调用gdalwarp
gdalwarp.exe -geoloc -t_srs EPSG:4326 E:\Tablefile\back.vrt E:\Tablefile\warp.tif
# vrt文件内容
<MDI key="LINE_OFFSET">2</MDI>
<MDI key="LINE_STEP">5</MDI>
<MDI key="PIXEL_OFFSET">2</MDI>
<MDI key="PIXEL_STEP">5</MDI>
<MDI key="SRS">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</MDI>
<MDI key="X_BAND">1</MDI>
<MDI key="X_DATASET">HDF4_SDS:UNKNOWN:"E:\\Tablefile\\modis\\MOD03.A2020350.0325.061.2020350090352.hdf":1</MDI>
<MDI key="Y_BAND">1</MDI>
<MDI key="Y_DATASET">HDF4_SDS:UNKNOWN:"E:\\Tablefile\\modis\\MOD03.A2020350.0325.061.2020350090352.hdf":0</MDI>
3. VRT校正的调用代码
def bulidVRT(warpFile, gltFile, vrtFile, outfile, warpWord):
'''
@desc:VRT调用gdalwarp.exe几何校正
'''
datasets = gdal.Open(warpFile)
SubDatasetsNum = len(datasets.GetSubDatasets())
infoList = []
for i in range(SubDatasetsNum):
eachLine = datasets.GetSubDatasets()[i]
if warpWord in eachLine[1]: # 需要校正的数据集的路径
infoList.append(eachLine[0])
if 'Latitude' in eachLine[1]:
infoList.append(eachLine[0])
if 'Longitude' in eachLine[1]:
infoList.append(eachLine[0])
print(infoList)
datasets2 = gdal.Open(gltFile)
SubDatasetsNum2 = len(datasets2.GetSubDatasets())
infoList2 = []
for i in range(SubDatasetsNum2):
eachLine2 = datasets2.GetSubDatasets()[i]
if 'Latitude' in eachLine2[1]:
infoList2.append(eachLine2[0])
if 'Longitude' in eachLine2[1]:
infoList2.append(eachLine2[0])
print(infoList2)
cmd = 'D:/Anaconda/Lib/site-packages/osgeo/gdal_translate.exe -of VRT %s %s' % (
infoList[0], vrtFile
)
os.system(cmd)
f = open(vrtFile)
newf = open(vrtFile.replace('.vrt', 'new.vrt'), 'w')
lines = f.readlines()
for line in lines:
if ('<MDI key="X_DATASET">' in line):
newf.write(' <MDI key="X_DATASET">'+ str(infoList2[1]) +'</MDI>\n')
elif ('<MDI key="Y_DATASET">' in line):
newf.write(' <MDI key="Y_DATASET">'+ str(infoList2[0]) +'</MDI>\n')
else:
newf.write(line)
newf.close()
f.close()
cmd2 = 'D:/Anaconda/Lib/site-packages/osgeo/gdalwarp.exe -geoloc -t_srs EPSG:4326 %s %s'%(
vrtFile.replace('.vrt', 'new.vrt'), outfile
)
os.system(cmd2)
os.remove(vrtFile)
os.remove(vrtFile.replace('.vrt', '.tif.aux.xml'))
os.rename(vrtFile.replace('.vrt', 'new.vrt'), vrtFile)
4. 后记
- 思考mod03数据几何校正失败的原因
- 需要熟悉dom解析模块/re表达式，以提高文档的解析力
- 后续补充RPC校正、TPS校正的原理及调用
