1

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.

4

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][0], 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)[0]) + "-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.

  • Hi Kersten, i tried your code but I am getting an error in the band_ds=gdal.Open(hdf_ds.GetSubDatasets()[subdataset][0], gdal.GA_ReadOnly). I am not sure what is wrong with my processing. hdf_ds = "folder\\path\\MOD.yyyydoy.hdf", band_ds=gdal.Open(hdf_ds.GetSubDatasets()["250m 16 days NDVI"][0]) and it returns 'NoneType' object has no attribute 'GetSubDatasets'. I tried to open my file using % gdalinfo file.hdf in the commandline and it returns the information from my hdf file. – user32145 Mar 3 '16 at 5:24
  • Hi Kersten, it is working if I am using gdal bundled with QGIS. But with my separate installation of GDAL it is not working. – user32145 Mar 3 '16 at 7:22
  • Then your other installation likely doesn't support the HDF driver. check with gdalinfo --formats – Kersten Mar 3 '16 at 7:36
  • Hi Kersten, my GDAL driver supports HDF4 (ros), HDF4Image (rw+),HDF5 (ros), HDF5 (image). However, I am not sure if what GDAL installation my command prompt is accessing. I have GDAL bundled with QGIS and I also have separate installation of GDAL from GIS internals. – user32145 Mar 3 '16 at 7:59
  • Thanks for your script it works. I am using gdal 2.2.2 and it seems the GetGeoTransform() and GetProjection() are coorectly extracted from the source data, but not applied to the result. The output is unprojected. Thanks for your help. – Michael Oct 29 '17 at 3:25
7

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.

2

for extracting all the layers, I used this code

for %i IN (*.hdf) do gdal_translate -sds %i %i.tif

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.