2

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/

  • 1
    I'm just noticing that you're the author of a whole series of incrementally related questions. I suggest you install iPython and use it with the --pylab option to explore your data interactively. TAB-completion is very useful, and iPython implements Python's introspection options very well. – chryss Sep 23 '13 at 2:13
8

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).ReadAsRaster()

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

1

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!

  • The example you provided is not Python code if I am not wrong. Can you give an example of how would you use gdal_translate in Python? – multigoodverse Sep 23 '13 at 6: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.