How to extract subdataset from HDF raster?

Question

I am trying to extract a subdataset from HDF raster. Later on, I intend to make operations such as mosaicking and extracting raster averages from the extracted subdataset. The following code is supposed to extract the subdataset, but it seems to me that all it does is extracting the subdataset NAME and not the object itself:

import gdal
from gdalconst import *
filename= "Image.hdf"
dataset = gdal.Open(filename,GA_ReadOnly)
subdataset = dataset.GetSubDatasets()[1]
print(type(dataset))
print(type(subdataset))

Now, if you run this, the type of the dataset and subdataset objects will be printed out:

<class 'osgeo.gdal.Dataset'>
<type 'tuple'>

So, the extracted subdataset is just a tuple. How can I use it to do mosaicking? Am I missing something here? P.S. HDF rasters can be easily downloaded from this FTP site: ftp://ladsweb.nascom.nasa.gov/allData/5/MOD11A1/2012/193/

Answer 1

What you're seeing is just the API that GDAL provides to access an HDF4 dataset and its subdatasets. The output of dataset.GetSubDatasets() is a list of tuples. Here's an example from one of the MOD11A1 files:

In [1]: dataset.GetSubDatasets()
Out[1]: [
  ('HDF4_EOS:EOS_GRID:"MOD11A1.A2012193.h00v08.005.2012196013548.hdf.hdf":MODIS_Grid_Daily_1km_LST:LST_Day_1km',
   '[1200x1200] LST_Day_1km MODIS_Grid_Daily_1km_LST (16-bit unsigned integer)'),  
  ('HDF4_EOS:EOS_GRID:"MOD11A1.A2012193.h00v08.005.2012196013548.hdf.hdf":MODIS_Grid_Daily_1km_LST:QC_Day',    
  '[1200x1200] QC_Day MODIS_Grid_Daily_1km_LST (8-bit unsigned integer)'),  ...]

Each tuple has two strings. The first is the subdataset name, the second a descriptor. dataset.GetSubDatasets()[1] is therefore the tuple describing the second subdataset (indices start at 0), and the subdataset name can be retrieved as dataset.GetSubDatasets()[1][0]. To address the second subdataset (the QC flags for the LST dataset), you do:

subdatasets = dataset.GetSubDatasets()
mysubdataset_name = subdatasets[1][0]
mydata = gdal.Open(mysubdataset_name, gdal.GA_ReadOnly).ReadAsArray()

This returns a 1200x1200 numpy array, which you can mosaic in any way you prefer.

Answer 2

To extract a sub dataset you can use gdal_translate like this:

gdal_translate src_ds dst_ds

check the help page for some examples on how to use it but generally the syntax for hdf4 files is the following:

gdaltranslate  HDF4_EOS:EOS_GRID:"/path/to/file/MOD11A1.A2012193.h00v08.005.2012196013548.hdf":MODIS_Grid_Daily_1km_LST:LST_Day_1km out.tiff

If you want to create just the mosaic, you can use gdal_merge which :

will automatically mosaic a set of images. All the images must be in the same coordinate system and have a matching number of bands, but they may be overlapping, and at different resolutions. In areas of overlap, the last image will be copied over earlier ones.

And since gdal_merge.py is a python script, you can open it with your favorite text processor, study it and then use it as refernce!