Use gdal tool to correct fy3 series satellite data

This document describes how to use gdal tools to correct fy3 series satellite data. The data provided by the fy3 series satellites is generally delivered in hdf5 format. A typical data file name of fy3a and fy3b is as follows:

FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDFFY3B_MERSI_GBAL_L1_20130620_0600_1000M_MS.HDF

The following three sections are used to process the fy3 series data correction, which are data preprocessing, data subdataset construction, and data correction.

The gdal [1] tools used in this document mainly include gdalinfo (used to view data information) and gdal_translate (used to extract subdatasets and format conversion) and gdalwarp [2] (for Image Correction ). In addition, qgis is used to view image data.

1. data preprocessing

In the data pre-processing process, you can view the raw data and obtain information about the sub-dataset to be processed. First, use the gdalinfo tool to view the sub-dataset information in the HDF file, mainly to find the sub-dataset to be processed and two sub-datasets of longitude and latitude. Use the gdalinfo command as follows:

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

In the preceding command, the-nomd parameter is used to indicate that no metadata information is output. If this parameter is not added, too much metadata is output. Therefore, the-nomd parameter is added here, no metadata is output. After the complete command is executed, HDF data information will be output in the fy3a_mersi_gbal_l0000201003200000150_1000m_ms.txt file. The content is as follows. Note the three child datasets in red text below. The Fourth Child dataset is the Child dataset to be corrected, the following 11 and 12 subdatasets are subdatasets used to store longitude and latitude coordinates. These subdatasets are used later.

Driver:HDF5/Hierarchical Data Format Release 5Files:D:\Data\FY3_Data\FY3A_MERSI_GBAL_L1_20100321_0150_1000M_MS.HDFSizeis 512, 512CoordinateSystem 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)

Note:When you use gdalinfo or other tools for processing, warning information of some columns will be displayed, as shown in 1. Ignore this information without affecting subsequent processing.

Figure 1 warning output by gdalinfo

2. Construct a subdataset

Through the first step above, we can find the child dataset to be corrected (that is, the fourth child dataset ). Next, we will export this subdataset for subsequent processing. Gdal_translate is used for export. The export format is VRT [3] (This format is a virtual file format based on XML format, which can be opened in Notepad for modification ). The Export command is as follows:

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

After the preceding command processing is complete, a VRT file named fy3a_mersi_gbal_l0000201003200000150_1000m_ms_sub5.vrt is generated, and the VRT file is opened by using notepad and other text editing software, as shown in figure 2. You can use qgis to open the file directly to display the image.

Figure 2 use NotePad to open a VRT File

Locate the </metadata> and <gcplist> nodes in the VRT file, roughly 88 rows, as shown in 3.

Figure 3 </metadata> and <gcplist> node locations

Add the following content between the two nodes, as shown in VRT file 4 after modification. The content added to the VRT file is called the geolocation metadata, which consists of nine subnodes. The details are as follows.

1. line_offset: Specifies the row offset.

2. line_step: Specifies the row step size.

3. pixel_offset: Specify the object offset.

4. pixel_step: Specify the pixel step size.

5. SrS: Specify the space reference information. Generally, it is a WGS84 coordinate system, which uses the WKT string format.

6. x_band: Specify the band number corresponding to the longitude

7. x_dataset: Specify the sub-dataset path corresponding to the longitude, that is, the sub-dataset 12 above.

8. y_band: Specify the band number corresponding to the longitude.

9. y_dataset: Specify the sub-dataset path corresponding to the longitude, that is, the sub-dataset 11 above.

The output information of gdalinfo in step 1 shows that the longitude and latitude data in fy3 is 2000X2048, which is consistent with the image size in subdataset 4, therefore, set line_step and pixel_step to 1. The data size of the precision and latitude of the image data is 406x271, and the image data may be 2030x1354, which is almost five times the relationship between the two, therefore, the preceding line_step and pixel_step must be set to 5 for the data of the patients who need to perform the operation. However, this step is not required for the data of the MNS. For more information, see the topic at the end of this article.

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

Figure 4 modified VRT File

Follow the steps above to modify the VRT file and save it. The VRT file will be used for subsequent processing.

3. Correction

After modifying the VRT file, we can use the gdalwarp tool to correct the VRT file. The specific command is as follows:

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

After executing the preceding command, you can correct the fourth subdataset to the WGS84 Coordinate System of longitude and latitude. If the output format is geotiff, metadata in the original subdataset may be lost (if necessary, it is parsed from the original VRT file ). The data before correction is opened using qgis 5. The corrected TIF data is opened using qgis and displayed in combination with a global SHP data, as shown in Figure 6 and Figure 7.

As shown in figure 6 and Figure 7, the use of the tool provided by gdal for fy3 data should be able to achieve the expected results, but there are some deviations in accuracy, 7. There is a certain deviation between the land and ocean boundaries near the Bohai Bay and SHP data. However, the precision of the low-resolution meteorological satellite fy3 should be able to meet the requirements.

Figure 5 pre-correction data

Figure 6 corrected data

Figure 7 zoom in after correction

Digress:

If you do not need to process the data in step 2, you can directly filter the child data to be corrected in step 1, and then use the gdalwarp tool in step 3 to correct the data. The sub-data set of the motile data directly contains the geolocation metadata domain, while the sub-data of the fy3 series does not contain the geolocation metadata domain, therefore, in the second step, we need to manually add a geolocation metadata field to use gdalwarp for processing.

I personally think this is a problem for the gdal database to parse the data of the fy3 series satellites. The gdal database directly constructs a geolocation meta-data domain for the parsing of the modemdata, but not the data of the fy3 series satellites, it is only parsed according to common HDF data.

At last, you can refer to the two blog posts I wrote earlier, using gdal to process the modem' data [4] and a geoloc Correction Code [5].

4. References

[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