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校正的原理及調用