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.
If you are working with a GDAL version that supports HDF5 datasets this is how you can extract a single subdataset from it and convert it to Geotiff in Python. This example was used for MOD09 and MOD13 band extraction (hence the no_data value conversion).
from osgeo import gdal import numpy as np def hdf_subdataset_extraction(hdf_file, dst_dir, subdataset): """unpack a single subdataset from a HDF5 container and write to GeoTiff""" # open the dataset hdf_ds = gdal.Open(hdf_file, gdal.GA_ReadOnly) band_ds = gdal.Open(hdf_ds.GetSubDatasets()[subdataset], gdal.GA_ReadOnly) # read into numpy array band_array = band_ds.ReadAsArray().astype(np.int16) # convert no_data values band_array[band_array == -28672] = -32768 # build output path band_path = os.path.join(dst_dir, os.path.basename(os.path.splitext(hdf_file)) + "-sd" + str(subdataset+1) + ".tif") # write raster out_ds = gdal.GetDriverByName('GTiff').Create(band_path, band_ds.RasterXSize, band_ds.RasterYSize, 1, #Number of bands gdal.GDT_Int16, ['COMPRESS=LZW', 'TILED=YES']) out_ds.SetGeoTransform(band_ds.GetGeoTransform()) out_ds.SetProjection(band_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(band_array) out_ds.GetRasterBand(1).SetNoDataValue(-32768) out_ds = None #close dataset to write to disc
edit: This answer used GDAL version 1.9.
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.
–ot Float32 for floating point we are in Integer.
for extracting all the layers, I used this code
for %i IN (*.hdf) do gdal_translate -sds %i %i.tif