使用GDAL工具对FY3系列卫星数据进行校正

本文档主要对如何使用GDAL提供的工具对FY3系列卫星数据进行校正处理。FY3系列卫星提供的数据一般是以HDF5格式下发，一个典型的FY3A和FY3B的数据文件名如下：

FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF
FY3B_MERSI_GBAL_L1_20130620_0600_1000M_MS.HDF

下面分为三个部分来对FY3系列数据校正进行处理，分别是数据预处理、构造子数据集和校正三个部分，下面分别进行详述。

该文档用到的GDAL^[1]工具主要有三个，gdalinfo（用于查看数据信息），gdal_translate（用于提取子数据集和格式转换等）和gdalwarp^[2]（用于图像校正）。此外使用了QGIS软件，用来查看图像数据。

1. 数据预处理

在数据预处理过程中，主要是对原始数据进行查看，得到要处理的子数据集等信息。首先使用gdalinfo工具查看HDF文件中的子数据集信息，主要是找到要处理的子数据集和两个经纬度子数据集。使用gdalinfo命令如下：

gdalinfo.exeD:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF -nomd >D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.txt

上面的命令中，使用了-nomd参数，表示不输出元数据信息，如果不加这个参数的话，输出的元数据太多了，所以这里还是加上这个-nomd参数，不输出元数据信息。执行完上面的命令后，会将HDF数据信息输出在文件FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.txt中，内容如下所示，注意下面红色字体的三个子数据集，其中第四个子数据集是我们要进行校正处理的子数据集，后面的11、12子数据集就是用来存储经度和纬度坐标的子数据集，后面会用到这三个子数据集。

Driver:HDF5/Hierarchical Data Format Release 5
Files:D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF
Sizeis 512, 512
CoordinateSystem is `'
Subdatasets:
 SUBDATASET_1_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://BB_DN_average
  SUBDATASET_1_DESC=[20x2000] //BB_DN_average(32-bit floating-point)
 SUBDATASET_2_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EVS_Attitude_angles
  SUBDATASET_2_DESC=[200x3]//EVS_Attitude_angles (64-bit floating-point)
  SUBDATASET_3_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EVS_orb_pos
  SUBDATASET_3_DESC=[200x3] //EVS_orb_pos(64-bit floating-point)
 SUBDATASET_4_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EVS_orb_vel
  SUBDATASET_4_DESC=[200x3] //EVS_orb_vel(64-bit floating-point)
  SUBDATASET_5_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EV_1KM_RefSB
 SUBDATASET_5_DESC=[15x2000x2048] //EV_1KM_RefSB (16-bit unsignedinteger)
  SUBDATASET_6_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EV_250_Aggr.1KM_Emissive
  SUBDATASET_6_DESC=[2000x2048]//EV_250_Aggr.1KM_Emissive (16-bit unsigned integer)
 SUBDATASET_7_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EV_250_Aggr.1KM_RefSB
  SUBDATASET_7_DESC=[4x2000x2048]//EV_250_Aggr.1KM_RefSB (16-bit unsigned integer)
 SUBDATASET_8_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Height
  SUBDATASET_8_DESC=[2000x2048] //Height(16-bit integer)
 SUBDATASET_9_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://LandCover
  SUBDATASET_9_DESC=[2000x2048] //LandCover(8-bit unsigned character)
 SUBDATASET_10_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://LandSeaMask
  SUBDATASET_10_DESC=[2000x2048] //LandSeaMask(8-bit unsigned character)
  SUBDATASET_11_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Latitude
 SUBDATASET_11_DESC=[2000x2048] //Latitude (32-bit floating-point)
 SUBDATASET_12_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Longitude
 SUBDATASET_12_DESC=[2000x2048] //Longitude (32-bit floating-point)
 SUBDATASET_13_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Moon_Vector
  SUBDATASET_13_DESC=[200x3] //Moon_Vector(32-bit floating-point)
 SUBDATASET_14_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://QA_index
  SUBDATASET_14_DESC=[2000x2048] //QA_index(32-bit unsigned integer)
 SUBDATASET_15_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Range
  SUBDATASET_15_DESC=[2000x2048] //Range(16-bit unsigned integer)
  SUBDATASET_16_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SV_DN_average
  SUBDATASET_16_DESC=[20x2000] //SV_DN_average(32-bit floating-point)
 SUBDATASET_17_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SensorAzimuth
  SUBDATASET_17_DESC=[2000x2048]//SensorAzimuth (16-bit integer)
 SUBDATASET_18_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SensorZenith
  SUBDATASET_18_DESC=[2000x2048] //SensorZenith(16-bit integer)
 SUBDATASET_19_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SolarAzimuth
  SUBDATASET_19_DESC=[2000x2048] //SolarAzimuth(16-bit integer)
 SUBDATASET_20_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SolarZenith
  SUBDATASET_20_DESC=[2000x2048] //SolarZenith(16-bit integer)
 SUBDATASET_21_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Sun_Vector
  SUBDATASET_21_DESC=[200x3] //Sun_Vector(32-bit floating-point)
 SUBDATASET_22_NAME=HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://VOC_DN_average
  SUBDATASET_22_DESC=[20x2000] //VOC_DN_average(32-bit floating-point)
CornerCoordinates:
UpperLeft  (   0.0,    0.0)
LowerLeft  (   0.0,  512.0)
UpperRight (  512.0,    0.0)
LowerRight (  512.0,  512.0)
Center      ( 256.0,  256.0)

注意：使用gdalinfo等工具处理时，会提示一些列的警告信息，如图 1所示，忽略这些信息即可，不影响后续处理。

图 1 gdalinfo输出的警告信息

2. 构造子数据集

通过上面的第一步，我们找到要校正的子数据集（就是第四个子数据集）。接下来我们要把这个子数据集单独导出来用来后续处理。导出使用gdal_translate工具，导出格式使用VRT格式^[3]（该格式是一种基于XML格式的虚拟文件格式，可以使用记事本打开进行修改）。导出命令如下：

gdal_translate.exe-of VRT HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://EV_1KM_RefSBFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS_sub5.vrt

经过上面的命令处理结束后，会生成一个文件名为FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS_sub5.vrt的VRT文件，使用记事本等文本编辑软件打开该VRT文件，如图 2所示。可以使用QGIS直接打开该文件，即可显示图像。

图 2 使用记事本打开VRT文件

在VRT文件中找到</Metadata>和<GCPList>节点左右，大致为88行左右，如图 3所示。

图 3 </Metadata>和<GCPList>节点位置

接下来在这两个节点之间加入下面的节点内容，修改后的VRT文件如图 4所示。加入VRT文件的内容叫GEOLOCATION元数据，主要由九个子节点组成，下面分别进行详细说明。

1、 LINE_OFFSET：指定行偏移量

2、 LINE_STEP：指定行步长

3、 PIXEL_OFFSET：指定象元偏移量

4、 PIXEL_STEP：指定象元步长

5、 SRS：指定空间参考信息，一般都为WGS84经纬度坐标系统，使用WKT字符串格式。

6、 X_BAND：指定经度对应的波段序号

7、 X_DATASET：指定经度对应的子数据集路径，就是上面的子数据集12

8、 Y_BAND：指定经度对应的波段序号

9、 Y_DATASET：指定经度对应的子数据集路径，就是上面的子数据集11

从第一步中gdalinfo输出的信息可以看出，FY3数据中的经度和纬度数据的大小为2000x2048，与子数据集四的图像大小一致，所以上面的LINE_STEP和PIXEL_STEP均设置为1即可。对于MODIS的数据，精度和纬度的数据大小为406x271，而图像数据有可能是2030x1354，这两者差不多是5倍的关系，所以对于MODIS数据来说，所以上面的LINE_STEP和PIXEL_STEP需要设置为5。不过对于MODIS数据来说，不要这一步，具体参考本文结尾处的题外话。

 <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:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Longitude</MDI>
    <MDI key="Y_BAND">1</MDI>
    <MDI key="Y_DATASET">HDF5:"D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Latitude</MDI>
  </Metadata>

图 4 修改后的VRT文件

按照上面的步骤修改完VRT文件后保存即可，该VRT文件用于后续处理。

3. 校正处理

通过对上面的VRT文件修改之后，我们就可以对该VRT文件使用gdalwarp工具来进行校正处理，具体命令如下：

gdalwarp.exe-geoloc D:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS_sub5.vrtD:DataFY3_DataFY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS_sub5_warp.tif

执行完上面的命令后，即可将第四个子数据集校正到WGS84经纬度坐标系统下。输出的格式为GeoTiff格式，可能会丢失原来子数据集中的元数据信息（如果需要用到的话，从原始VRT文件中进行解析）。校正前的数据使用QGIS打开如图 5所示。校正后的tif数据使用QGIS打开并与一个全球的shp数据叠加显示，如图 6和图 7所示。

从图 6和图 7中看出，对于FY3的数据使用GDAL提供的工具校正应该可以达到预期的效果，但是对于精度有些偏差，如图 7中渤海湾附近的陆地和海洋边界与shp数据有一定的偏差，不过对于FY3这种低分辨率的气象卫星这种精度应该可以满足要求。

图 5 校正前数据

图 6 校正后数据

图 7 校正后放大显示

题外话：

对于MODIS数据，不需要经过第二步处理，直接通过第一步筛选要校正的子数据集，然后直接用第三步中的gdalwarp工具校正即可。MODIS数据的子数据集中直接就包含了GEOLOCATION元数据域，而FY3系列的数据，子数据集中没有包含GEOLOCATION元数据域，所以在第二步我们需要人工来添加一个GEOLOCATION元数据域，从而才能使用gdalwarp进行处理。

个人觉得这是GDAL库对于FY3系列卫星的数据解析问题，GDAL库对于MODIS数据解析时就直接构建了一个GEOLOCATION元数据域，而对于FY3系列卫星的数据没有构建，只是按照普通的HDF数据来进行解析。

最后可以参考我之前写的两篇博文，使用GDAL处理MODIS数据^[4]和一个Geoloc校正的代码^[5]。

4. 参考资料

[1] www.gdal.org

[2] http://www.gdal.org/gdalwarp.html

[3] http://www.gdal.org/gdal_vrttut.html

[4] http://blog.csdn.net/liminlu0314/article/details/8623607

[5] http://blog.csdn.net/liminlu0314/article/details/8629298