FY-4A批量几何校正（基于GDAL）

Posted 2021-06-11 Rser_ljw

引言

​ 前次介绍了FY-4A遥感影像的几何校正方法，然而其主体基于ENVI实现，虽然可以通过ENVI_DO_IT和IDL语言实现批量化处理，但由于ENVI中GLT的生成需要经纬度查找表没有无效值，需要裁剪，所以比较麻烦。

​ 下面介绍一下基于gdal进行FY-4A几何校正的方法。仔细阅读李民录老师《gdal实现静止卫星几何校正》相关文章后，发现其主要流程都是先gdalinfo查看数据集，再gdaltranslate生成vrt，然后修改vrt添加查找表数据路径，最后gdalwarp进行几何校正。

​ 实施过程中也发现了一些问题，如圆盘数据的类型不同，有HDF/NC/Tiff，存储方式不同，读取时文件的路径也各不相同，有的文件类型如NC甚至无法生成VRT文件。综合分析后，笔者决定拆解这一问题，既然gdalwarp的核心是根据影像数组，经纬度查找表数组进行几何校正，那就可以只凭借这三个数组实现几何校正，故可以生产各自的geotiff文件，略去多余的信息。

步骤

1. 数据准备

   1.1 遥感影像nc转geotiff

   ​ FY-4A遥感影像属于netCDF数据集，在使用gdal_translate.exe生成vrt文件时，netCDF文件由于库自身的问题会导致错误，故首先需要把FY-4A数据转换成GeoTiff数据。

​ 1.2 生成lon/lat对应的geotiff数据

​ 即经纬度查找表

2. gdal_translate生成VRT文件

   有关gdal的命令如gadl_translate.exe/gdalwarp.exe可以在cmd里调用，也可以在python环境下调用gdal函数。（需要注意在cmd里调用需要声明插件的路径（proj.db/...），或者在环境变量中声明。）

gdal_translate.exe -of VRT E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif E:/Tablefile/before.vrt

gdal.Translate(\'E:/Tablefile/before.vrt\', \'E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif\', options=gdal.TranslateOptions(format=\'VRT\'))

3. 修改VRT文件，添加用于校正的lon/lat数据集

   以上操作生成VRT文件(beforer.vrt)内容如下：

<VRTDataset rasterXSize="2748" rasterYSize="2748">
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

遥感影像有几个波段，就会生成几个VRTRasterBand，部分数据（单波段、部分多波段）不会生成"TIFFTAG_XRESOLUTION"信息（情况暂时不明），需要自行添加，重复添加不会有影响。（记得在SourceFilename处添加波段路径）

本次校正修改如下：

添加经纬度查找表数据描述信息（在X_DATASET/Y_DATASET处添加lon/lat文件路径）

Modis数据由于自带lon/lat数据集，所以会自动生成这段信息；而其他不含lon/lat的卫星数据（如FY-4A、单独的geotiff）需要手动添加这段信息。

  <Metadata>
    <MDI key="TIFFTAG_XRESOLUTION">1</MDI>
    <MDI key="TIFFTAG_YRESOLUTION">1</MDI>
  </Metadata>
  <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">E:\\Tablefile\\CBHRegression\\geometryCorrect\\lookuptable\\4km\\lonDisk.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">E:\\Tablefile\\CBHRegression\\geometryCorrect\\lookuptable\\4km\\latDisk.tif</MDI>
  </Metadata>

修改后VRT文件(after.vrt)如下：

此处可以更改dataType，设置数据类型。

<VRTDataset rasterXSize="2748" rasterYSize="2748">
  <Metadata>
    <MDI key="TIFFTAG_XRESOLUTION">1</MDI>
    <MDI key="TIFFTAG_YRESOLUTION">1</MDI>
  </Metadata>
  <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">E:\\Tablefile\\CBHRegression\\geometryCorrect\\lookuptable\\4km\\lonDisk.tif</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">E:\\Tablefile\\CBHRegression\\geometryCorrect\\lookuptable\\4km\\latDisk.tif</MDI>
  </Metadata>
  <VRTRasterBand dataType="Float32" band="1">
    <NoDataValue>-999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="2">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>2</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
  <VRTRasterBand dataType="Float32" band="3">
    <NoDataValue>-999</NoDataValue>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">E:/Tablefile/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20210609130000_20210609131459_4000M_V0001.tif</SourceFilename>
      <SourceBand>3</SourceBand>
      <SourceProperties RasterXSize="2748" RasterYSize="2748" DataType="Float32" BlockXSize="2748" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
      <DstRect xOff="0" yOff="0" xSize="2748" ySize="2748" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

4. gdalwarp几何校正

最后利用gdalwarp进行几何校正即可（推荐使用函数进行）

注意指定geoloc利用查找表进行校正、利用t_srs添加投影信息，这样会自动生成含有geotrans/proj的geotiff。

gdalwarp.exe -geoloc -t_srs EPSG:4326 E:/Tablefile/after.vrt E:/Tablefile/warp.tif
# 但是这种方式有时会找不到插件路径，需要自行设置

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
gdal.Warp(r\'E:/Tablefile/warp.tif\',
                r\'E:/Tablefile/after.vrt\',
                format=\'GTiff\', geoloc=True,
                dstSRS=srs.ExportToWkt(), resampleAlg=gdal.GRIORA_NearestNeighbour)

结果

几何校正结果如下，黄色是-180_180，-9090的全球地图框线，红色是东海附近框线。由于FY-4A只拍到了地球的一侧，所以只有右半区有成像，左半边的黑色区域对应圆盘的最右侧，其经度起于最左侧地区。（此时中国正在夜间）

后记

1. 几何校正的质量和经纬度查找表的质量直接相关，错误的经纬度查找表会导致校正失败。（注意查找表的数据格式，必须是单波段的float32数组）

2. FY-4A有HDF，也有NC数据集，转换TIFF时需要注意；此外部分数据如CTH读取后数组需要.astype(np.float32)，否则在gdal写入栅格数据时gdal内部会发生数据类型转换导致错误。

3. 中国区REGC与全圆盘DISK操作方法相同，注意选择合适的经纬度查找表即可。

4. 该方法可以在已知经纬度查找表数据、遥感影像数据的情况下进行几何校正。将以上过程代码化，即可完成遥感影像的批量几何校正。