本文档主要对如何使用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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF -nomd >D:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF

Sizeis 512, 512

CoordinateSystem is `'

Subdatasets:

 SUBDATASET_1_NAME=HDF5:"D:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Height

  SUBDATASET_8_DESC=[2000x2048] //Height(16-bit integer)

 SUBDATASET_9_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://LandCover

  SUBDATASET_9_DESC=[2000x2048] //LandCover(8-bit unsigned character)

 SUBDATASET_10_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://LandSeaMask

  SUBDATASET_10_DESC=[2000x2048] //LandSeaMask(8-bit unsigned character)

  SUBDATASET_11_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Latitude

 SUBDATASET_11_DESC=[2000x2048] //Latitude (32-bit floating-point)

 SUBDATASET_12_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Longitude

 SUBDATASET_12_DESC=[2000x2048] //Longitude (32-bit floating-point)

 SUBDATASET_13_NAME=HDF5:"D:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Range

  SUBDATASET_15_DESC=[2000x2048] //Range(16-bit unsigned integer)

  SUBDATASET_16_NAME=HDF5:"D:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SensorAzimuth

  SUBDATASET_17_DESC=[2000x2048]//SensorAzimuth (16-bit integer)

 SUBDATASET_18_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SensorZenith

  SUBDATASET_18_DESC=[2000x2048] //SensorZenith(16-bit integer)

 SUBDATASET_19_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SolarAzimuth

  SUBDATASET_19_DESC=[2000x2048] //SolarAzimuth(16-bit integer)

 SUBDATASET_20_NAME=HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://SolarZenith

  SUBDATASET_20_DESC=[2000x2048] //SolarZenith(16-bit integer)

 SUBDATASET_21_NAME=HDF5:"D:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_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:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Longitude</MDI>

    <MDI key="Y_BAND">1</MDI>

    <MDI key="Y_DATASET">HDF5:"D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDF"://Latitude</MDI>

  </Metadata>

图4 修改后的VRT文件

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

3.   校正处理

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

gdalwarp.exe-geoloc D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS_sub5.vrtD:\Data\FY3_Data\FY3A_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

最新文章

  1. Redis为什么使用单进程单线程方式也这么快
  2. 月四 周2 ii
  3. 2015 Android Dev Summit(安卓开发峰会)第一天
  4. SQL2000的Enterprise Edition和Developer Edition有什么区别
  5. [转]ubuntu 14.04 系统设置不见了
  6. sgu 107 987654321 problem
  7. angularjs 遇到Error: [$injector:unpr] Unknown provider: tdpicnews-serviceProvider &lt;- tdpicnews-service &lt;- tdpic-controller 错误
  8. Linux系统安全需要注意的一些问题
  9. 第一个asp.net实例——生日邀请以及回函
  10. BZOJ_3747_[POI2015]Kinoman_线段树
  11. 记录 制作校园网登陆脚本 python编写 附源码
  12. Wannafly summer camp Day6 - D 区间权值
  13. 编译安装centos7 php7.2 mysql5.7 nginx1.9.9
  14. android selector shape 使用
  15. 5. Longest Palindromic Substring 返回最长的回文子串
  16. 让shell脚本中的echo输出带颜色
  17. DB2如何将数据库表解锁
  18. JAVA遇上HTML-----JSP 篇基本概念
  19. C#中的装箱拆箱
  20. Android 动态添加Spinner(.java文件内实现) 实现 改变spinner 内文字属性

热门文章

  1. Android文件大头10G
  2. Android 多窗口详解
  3. 实现Android5.0过渡动画兼容库
  4. tomcat启动批处理——startup.bat
  5. Linux Debugging （九） 一次生产环境下的“内存泄露”
  6. FORM中读取图片
  7. Gradle 的Daemon配置
  8. Android Demo 下拉刷新+加载更多+滑动删除
  9. UNIX网络编程——客户/服务器程序设计示范（四）
  10. (NO.00005)iOS实现炸弹人游戏(七):游戏数据的序列化表示