5

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

UPDATED: It turns out that the gdal package on my computer is out of date so I cannt get the GCPs via gdalinfo. I updated the gdal to 1.10.0 and now everythings work fine.

  • Welcome to gis.SE. Where did you download the MODIS, and which particular image set are you using? What have you already looked at (modis-atmos.gsfc.nasa.gov/MOD04_L2/spec.html ?) What happens when you use gdalinfo? Do you see GRINGPOINTLATITUDE and GRINGPOINTLONGITUDE? – BradHards Dec 25 '13 at 5:03
  • 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. I updated the log from gdalinfo command above. And yes I see the GRINGPOINTLATITUDE and GRINGPOINTLONGITUDE but dont known what it's for? – Hung Luu Dec 25 '13 at 7:00
  • Can you please edit your question, rather than adding comments? That helps us to understand the "whole story" and helps other people. GringPointLatitude and GringPointLongitude give you the geographic bounds of the swath - see nrm.salrm.uaf.edu/~dverbyla/MODIS_Land_Products/pdf_documents/… and look for the stuff that is common. – BradHards Dec 25 '13 at 9:20
  • Thanks @BradHards, I updated the question. As you said, GringPointLatitude and GringPointLongitude give the geographic bounds of the swath, but I still cant see how I could you it for the reprojection. Can I use these points as corresponding to MODIS image corner coordinate as GCPs? – Hung Luu Dec 25 '13 at 9:36
  • I believe those are the corner points, so why not use those as part of the warp? – BradHards Dec 25 '13 at 22:16
10

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:

enter image description here


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:

enter image description here

  • Thanks @Andre Joost, your approach is absolutely right. It turns out that the gdal package on my computer is out of date so I cannt get the GCPs via gdalinfo. I updated the gdal to 1.10.0 and now everythings work fine. – Hung Luu Dec 29 '13 at 18:09
  • Testing with other sample data I found out that -tps option is required for the non-equidistant raster source. Otherwise your data gets slight offsets. – AndreJ Dec 29 '13 at 20:28
  • Hi Andre Joost, there are some files which some GCPs have value kind of (0.5,126.5) -> (-999,-999,0). Is it normal since if I gdalwarp those file, the output tif is extremely large (about 174mb) and cant be processed? Do you have any ideas? – Hung Luu Dec 30 '13 at 6:36
  • Can you give me an example file name? I would suggest to use the -refine_gcps option to filter those out. – AndreJ Dec 30 '13 at 7:04
  • example hdf file: sdrv.ms/1el4xCF gdalinfo: sdrv.ms/1el4AOR the subdataset is "mod04:Image_Optical_Depth_Land_And_Ocean" – Hung Luu Dec 30 '13 at 7:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.