coordinate system - Reprojecting MODIS Swath data to WGS84 using GDAL

Im new to GIS and now working with MODIS04L2 Aerosol Optical Thickness image. I need to reproject the image to other coordinate system such as WGS84 for later work but now Im stuck with it. I have try two ways:

1. Write my own code using python with gdal lib: I follow the tutorial which using GCPs point as paramater for gdalwarp command. The problem is I can't get the GCPs point. The dataset.GetGCPs() function return a zero-length array of GCPs.


2. Using programs like HEG tool, MRT tool but both of it have issues. When I open MODIS file (with extension .hdf) with MRT tool I get the error "Error in Module: ReadHdrFile".

I download MODIS from ladsweb.nascom.nasa.gov/data/search.html with parameter: Satellite/Instrument:Terra, Group: Terra Atmosphere Level 2 Product, Products: MOD04L2_Leve 2 Aerrosol and Collection: MODIS Collection 5.1. I want to extract subdataset "mod04:Image_Optical_Depth_Land_And_Ocean" from hdf MODIS file to geotiff which reprojected to WGS84.

The sample hdf file http://sdrv.ms/JpVQNh The log from gdalinfo:http://sdrv.ms/JpTQok

Anyone here has experiment with this process?

---

It turns out that the gdal package on my computer is out of date so I cannot get the GCPs via gdalinfo. I updated the gdal to 1.10.0 and now everything is working fine.

Answer

You need the subdataset full name from the query on the file:

gdalinfo MOD04_L2.A2003001.0005.051.2010313005421.hdf >2003.txt

With the subdataset name, you get the GCP coordinates in pixel and lon/lat:

gdalinfo HDF4_EOS:EOS_SWATH:"MOD04_L2.A2003001.0005.051.2010313005421.hdf":mod04:Image_Optical_Depth_Land_And_Ocean >>2003a.txt

With the following content:

Driver: HDF4Image/HDF4 Dataset
Files: MOD04_L2.A2003001.0005.051.2010313005421.hdf

Size is 135, 204
Coordinate System is `'
GCP[  0]: Id=, Info=
          (0.5,0.5) -> (146.491271972656,2.82822823524475,0)
GCP[  1]: Id=, Info=
          (12.5,0.5) -> (150.015686035156,2.33101844787598,0)
GCP[  2]: Id=, Info=
          (24.5,0.5) -> (152.117935180664,2.02903413772583,0)
....

Note that the source file raster is not equidistant in WGS84, so the 4 corner coordinates are not sufficient, and thin plate spline interpolation is needed.

Additionally, the subdatasets

SUBDATASET_70_NAME=HDF4_SDS:UNKNOWN:"MOD04_L2.A2003001.0005.051.2010313005421.hdf":0
SUBDATASET_70_DESC=[204x135] Longitude (32-bit floating-point)
SUBDATASET_71_NAME=HDF4_SDS:UNKNOWN:"MOD04_L2.A2003001.0005.051.2010313005421.hdf":1
SUBDATASET_71_DESC=[204x135] Latitude (32-bit floating-point)

contain longitudes and latitudes as raster values.

Apart from that, you can use gdalwarp to extract the subdataset using those GCP points:

gdalwarp -of GTIFF -tps -t_srs EPSG:4326 HDF4_EOS:EOS_SWATH:"MOD04_L2.A2003001.0005.051.2010313005421.hdf":mod04:Image_Optical_Depth_Land_And_Ocean 2003.tif

This file is already projected in WGS84.

---

EDIT

The second sample file you provide has a complete line of garbage at line 126.5. There is no way of automatically removing it, but you can create a VRT using gdal_translate, remove the false GCP entries manually from that, then gdalwarp to WGS84:

gdal_translate -of VRT HDF4_EOS:EOS_SWATH:"MOD04_L2.A2011001.0030.051.2011001135339.hdf":mod04:Image_Optical_Depth_Land_And_Ocean 2011.vrt
(edit 2011.vrt)
gdalwarp -of GTIFF -tps 2011.vrt 2011-3.tif

The result still contains a blank line, but the rest is usable:

---

EDIT 2

Another approach, that prevents GCP interpolation issues, is to export the data as XYZ file, and do the same with the longitude and latitude subdatasets:

gdal_translate -of XYZ HDF4_EOS:EOS_SWATH:"MOD04_L2.A2011001.0030.051.2011001135339.hdf":mod04:Image_Optical_Depth_Land_And_Ocean 2011.xyz
gdal_translate -of XYZ HDF4_SDS:UNKNOWN:"MOD04_L2.A2011001.0030.051.2011001135339.hdf":0 2011-lon.xyz
gdal_translate -of XYZ HDF4_SDS:UNKNOWN:"MOD04_L2.A2011001.0030.051.2011001135339.hdf":1 2011-lat.xyz

These text files can be put together in a spreadsheet, and imported as delimited text layers using the lon and lat columns. The garbage coordinates get placed outside the globe, and can easily be cut with a +/- 180/90° cutline polygon. I took this as reference to find out how to use gdalwarp parameters correctly: