使用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

原文地址:https://www.cnblogs.com/xiaowangba/p/6313948.html