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.
-
Have you viewed this post. I seem to remember that I did hdf to geotiff in windows with GDAL for some MODIS data in Windows. gis.stackexchange.com/questions/58688/…– If you do not know- just GISDec 16, 2015 at 2:55
-
Yes HDF is included with the OSGEO4W version of GDAL.– If you do not know- just GISDec 16, 2015 at 3:09
3 Answers
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.
-
-
-
hi, it does not work.. I get this - ERROR 6: Unknown option name '-sd'. If I write -sds I get: ERROR 6: Too many command options Oct 18, 2017 at 11:58
-
-
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. Mar 3, 2016 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. Mar 3, 2016 at 7:22
-
Then your other installation likely doesn't support the HDF driver. check with
gdalinfo --formats
– KerstenMar 3, 2016 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. Mar 3, 2016 at 7:59
-
Thanks for your script it works. I am using gdal 2.2.2 and it seems the
GetGeoTransform()
andGetProjection()
are coorectly extracted from the source data, but not applied to the result. The output is unprojected. Thanks for your help.– MichaelOct 29, 2017 at 3:25
for extracting all the layers, I used this code
for %i IN (*.hdf) do gdal_translate -sds %i %i.tif
-
using bash, this one worked:
for i in *.hdf; do gdal_translate -sds $i $i.tif; done
– maceFeb 21, 2020 at 15:10