[GIS] extract scientific layers from MODIS HDF Dataeset using Python GDAL

gdalhdfmodispython

I have been trying to extract scientific data layers from MODIS HDF files but without success. Before I have been doing this type of task using Arcpy, but I want the script that I am currently working to be based only in GDAL and not in Arcpy so I can easily share it.
I used GDAL, but it seems that GDAL version in windows does not have support for HDF files (I browsed some forums about it). I tried to open the HDF file using gdal.Open() but it does not return anything, and I am pretty sure that my directory path to the HDF file is correct.
Then I tried pyHDF. So this packages allowed me to extract the scientific layers as 2D array, but I have a problem how to save it as GeoTIFF file. I am aware that for me to save a 2D array to a GeoTIFF file I need some reference, such as projection and raster properties, just like the GetGeoTransform() module in GDAL. In this step, I am not sure how to extract that same kind of information using pyHDF.

Best Answer

Here is my code I use to convert all HDFs in a folder to GeoTIFF via GDAL (OSGEO4w) in Windows when working with SST data. Remember to use the OSGEO4W version to get hdf support.

for %A in ("C:\temp\*.hdf") do gdal_translate -of GTiff -a_srs "+init=epsg:4326" -a_ullr -180 90 180 -90 -co "COMPRESS=PACKBITS" -sds "%A" "%A.tiff

for %A in ("C:\temp\*.hdf") do Loop through all files in the folder with .hdf extension.

gdal_translate -of GTiff Convert all files in the selected folder to GTiff.

-a_srs Set the WGS 1984 projection to all the Gtiffs.

-a_ullr Set the bounds of all the Gtiffs.

-co Compress all the Gtiffs.

-sds Convert all files in the hdf, subdataset.

TIP

Use –ot Float32 for floating point we are in Integer.

Related Solutions

[GIS] processing MODIS hdf files with gdal

HDF4 support is not available in GDAL by default. If you're using the GDAL binaries and python libraries from GISInternals, these do not have HDF4 support compiled in (at the time this answer was originally posted in July 2013, however HDF4 support was added to the GISInternal GDAL build in Nov. 2013). HDF4 support is compiled in the OSGeo4W GDAL binaries.

You can test if you have HDF4 support with the gdalinfo command: gdalinfo --format hdf4

You can also batch the MRT from within your python script using the subprocess module, see the "Command Line Interface" section of the MRT User Manual.

Edit: the subprocess module is easy to use, see below. See also the "Automated Batch Processing" section in the MRT User Manual.

import subprocess
def runcmd(cmd)
    proc = subprocess.Popen(cmd, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
    stdout,stderr=proc.communicate()
    exit_code=proc.wait()
    return exit_code,stdout,stderr

cmd = 'mrt command line'
exit_code,stdout,stderr = runcmd(cmd)
if exit_code:raise RuntimeError(stderr)
else:print stdout

GDAL Python – Extracting Projection Info from HDF File

This is fairly straightforward if you think of the HDF dataset as a container, where each subdataset is a raster image with its own projection.

Your error is in not opening the subdataset, as GetSubDatasets only returns the strings you need to access them.

# open the HDF container
hdf_ds = gdal.Open(hdfFile)

# this is just a string of the name of the subdataset
b3_string = hdf_ds.GetSubDatasets()[4][0]

# open the subdataset
sds_b3 = gdal.Open(hdf_ds.GetSubDatasets()[4][0])

# get the projection
proj = sds_b3.GetProjection()

Best Answer

Related Solutions

[GIS] processing MODIS hdf files with gdal

GDAL Python – Extracting Projection Info from HDF File

Related Question