Python+GDAL几何校正任意自带经纬度数据的遥感影像

前言

李民录老师在他的博客中使用C++ GDAL的gdaltranslate.exe、gdalwarp.exe等工具对FY3A_MERSI数据进行了几何校正，思路是构建一个需要校正波段数据的VRT虚拟栅格文件，然后将FY数据集自带的地理位置数据写入VRT，然后利用warp进行几何校正。本篇博客主要介绍如何使用Python + GDAL完成这一实现，理论上自身带有地理位置数据的影像数据集都可用此方法进行几何校正，如FY-3系列、MODIS、FY-4A， Himawari-8数据。

构建虚拟数据集

通过查看官方手册，可知translate API主要是提供一个数据格式转换的功能。
TranslateOptions

在Python中通过以下代码便可将指定数据集或者子数据构造成一个VRT虚拟栅格文件。

 vrtFile = gdal.Translate(vrtDir, 
                         subDataset,
                         format='vrt')

根据上面提到的李民录老师的博客中对GEOLOCATION元数据的描述，向VRT文件中写入以下。关于各个key的说明，李老师博客中有详细解释。特别需要注意的是X_DATASET和Y_DATASET分别指定了经度和纬度数据的存储位置。

<Metadata domain="GEOLOCATION">
 <MDI key="LINE_OFFSET">1</MDI>
 <MDI key="LINE_STEP">1</MDI>
 <MDI key="PIXEL_OFFSET">1</MDI>
 <MDI key="PIXEL_STEP">1</MDI>
 <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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"]],AUTHORITY["EPSG","4326"]]</MDI>
 <MDI key="X_BAND">1</MDI> 
 <MDI key="X_DATASET">HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Longitude</MDI>
 <MDI key="Y_BAND">1</MDI>
 <MDI key="Y_DATASET">HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Latitude</MDI>
</Metadata>

几何校正

几何校正可利用warp实现，可参考之前的博客Python GDAL学习笔记(二)，校正的对象为上面构建的VRT虚拟栅格数据集。

完整代码

这里根据一景FY-3D MERSI-2 250m分辨率数据25波段为例，需要注意的该数据的波段数据和地理定位数据分别在两个HDF5文件中，所以定义函数的时候，有两个输入文件的路径。

# -*- coding: utf-8 -*-
import gdal, osr
import os 

def geoMERSI2(file, geoFile, afterGeoFile):
    '''
    
    Parameters
    ----------
    file : 文件绝对路径
        需要校正的文件.
    geoFile : 文件绝对路径
        地理坐标存在的文件.
    afterGeoFile : tif文件绝对路径
        经过校正后的tif格式 Albers投影 250m空间分辨率.
    Returns
    -------
    None.

    '''
    #os.chdir(r'F:\Py_project\gdalDealMERSI2')  
    dataset = gdal.Open(file) 
    # geoDataset = gdal.Open(geoFile)    
    #print(gdal.Info(dataset)) # 查看元数据
    #print(gdal.Info(geoDataset))
    #subdataset = dataset.GetSubDatasets() # 查看子数据集
    #print(subdataset)
    '''
    查看子数据集位置：
       SUBDATASET_6_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_0250M_MS.hdf"://Data/EV_250_Emissive_b24
       SUBDATASET_6_DESC=[8000x8192] //Data/EV_250_Emissive_b24 (16-bit unsigned integer) 
    '''
    '''
    地理校正数据的位置：
      SUBDATASET_1_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Latitude
      SUBDATASET_1_DESC=[8000x8192] //Latitude (32-bit floating-point)
      SUBDATASET_2_NAME=HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Longitude
      SUBDATASET_2_DESC=[8000x8192] //Longitude (32-bit floating-point)  
    '''
    vrtDir = os.path.splitext(file)[0] + '.vrt'
    subDataset = dataset.GetSubDatasets()[6][0] # B25
    vrtFile = gdal.Translate(vrtDir, 
                         subDataset,
                         format='vrt')
    '''
    需要写入描述GEOLOCATION元数据域的信息：
        <Metadata domain="GEOLOCATION">
         <MDI key="LINE_OFFSET">1</MDI>
         <MDI key="LINE_STEP">1</MDI>
         <MDI key="PIXEL_OFFSET">1</MDI>
         <MDI key="PIXEL_STEP">1</MDI>
         <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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"]],AUTHORITY["EPSG","4326"]]</MDI>
         <MDI key="X_BAND">1</MDI>
         <MDI key="X_DATASET">HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Longitude</MDI>
         <MDI key="Y_BAND">1</MDI>
         <MDI key="Y_DATASET">HDF5:"F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf"://Latitude</MDI>
        </Metadata> 
    '''
    srs = osr.SpatialReference()
    srs.ImportFromProj4('+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
    
    lines = []
    with open(vrtDir, 'r') as f:
        for line in f:
            lines.append(line) 
    lines.insert(1,'<Metadata domain="GEOLOCATION">\n')
    lines.insert(2,' <MDI key="LINE_OFFSET">1</MDI>\n')
    lines.insert(3, ' <MDI key="LINE_STEP">1</MDI>\n')
    lines.insert(4, ' <MDI key="PIXEL_OFFSET">1</MDI>\n')
    lines.insert(5, ' <MDI key="PIXEL_STEP">1</MDI>\n')
    lines.insert(6, ' <MDI key="SRS">GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",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"]],AUTHORITY["EPSG","4326"]]</MDI>\n')
    lines.insert(7, ' <MDI key="X_BAND">1</MDI>')
    lines.insert(8, ' <MDI key="X_DATASET">HDF5:"{}"://Longitude</MDI>\n'.format(geoFile))
    lines.insert(9, ' <MDI key="Y_BAND">1</MDI>\n')
    lines.insert(10, ' <MDI key="Y_DATASET">HDF5:"{}"://Latitude</MDI>\n'.format(geoFile))
    lines.insert(11, '</Metadata>\n')
    with open(vrtDir, 'w') as f:
        for line in lines:
            f.writelines(line)
    '''
    geoData = gdal.Warp('F:\Py_project\gdalDealMERSI2\Warp_B24_WGS.tif',  
              'F:\Py_project\gdalDealMERSI2\MERSI2_0403_B24_vrt.vrt', 
              format='GTiff', geoloc=True, dstSRS="EPSG:4326",
              xRes=0.25, yRes=0.25)
    '''
    geoData = gdal.Warp(afterGeoFile, vrtDir, 
              format='GTiff', geoloc=True, dstSRS=srs,
              resampleAlg=gdal.GRIORA_Bilinear, xRes=250, yRes=250)
    
    if geoData == None:
        print('deal failure!')
    del geoData
    print('{} finish Geo\n'.format(file)) 
    
if __name__ == '__main__':
    file = r'F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_0250M_MS.hdf'
    geoFile = r'F:\FY3D_MERSI_data\FY3D_MERSI_GBAL_L1_20190403_0445_GEOQK_MS.hdf'
    afterGeoFile = r'F:\Py_project\gdalDealMERSI2\Warp_B25_albers.tif'
    geoMERSI2(file, geoFile, afterGeoFile)

校正结果如下，没有去除bowties，但速度比ENVI的GLT校正快很多。