I have an issue with georeferencing of HDF files for the MYD04 product from MODIS (Dark Target Aerosol).

A sample HDF file is given in the following link -


Dragging and dropping the above file into QGIS, then choosing one of the layers, places it more or less in the right location. However, a closer look suggests the geo-referencing is not accurate. Seems like the image is a little too stretched. For example, note how the Gulf of Aqaba and the Nile river are misplaced -

enter image description here

Here is a zoomed-in view -

enter image description here

Two more things I tried -

  1. using gdal_translate to convert from HDF to TIF with -

gdal_translate -of GTiff \'HDF4_EOS:EOS_SWATH:MYD04_L2.A2002187.1050.006.2014195211826.hdf:mod04:Deep_Blue_Aerosol_Optical_Depth_550_Land\' modis.tif

This works, but the resulting TIF image has exactly the same misplacement issue as when uploading the HDF into QGIS.

  1. Reading the raster with R, using the gdalUtils package, in which case the extent information is not even being read and I see a rectangular image such as this -

enter image description here

  1. I tried looking for something useful in the image metadata, attached here -


I see that it contains GCPs, which I have no idea how to use in a better way than what QGIS does.

My conclusion from all of this is that the HDF image is not reoreferenced via a rectangular extent (xmin, ymin, xmax, ymax), but rather through some set of GCPs that require resampling. However I have no idea what is the cause for the geo-referencing problem - do QGIS and GDAL read the metadata in a wrong way? Are the GCPs inaccurate?

Following the valuable suggestions by AndreJ (see his comment), I used the 'lon' and 'lat' reference grid embedded in layers 71 & 72 of the HDF to create a point layer. Using R code as follows -


filename = "MYD04_L2.A2002187.1050.006.2014195211826.hdf"

# Subdataset
sds = get_subdatasets(filename)
s = sds[grepl("Deep_Blue_Aerosol_Optical_Depth_550_Land$", sds)]

# Values
r = s %>% readGDAL %>% raster

# Lon-Lat
lon = sds[71] %>% readGDAL %>% raster
lat = sds[72] %>% readGDAL %>% raster

# Stack
r = stack(lon, lat, r)
names(r) = c("lon", "lat", "value")

# To points
dat = as.data.frame(r)
dat = st_as_sf(dat, coords = c("lon", "lat"), crs = 4326)
st_write(dat, "test.geojson")

Here is what the layer looks like, now correctly positioned -

enter image description here

  • 2
    gdal_translate does not use the GCP, nor the lat and lon bands in sds 0 and 1. You have to use gdalwarp instead. See gis.stackexchange.com/questions/81361/… for a full workflow.
    – AndreJ
    May 16 '18 at 18:31
  • Great, this solved the problem! Thank you very much May 17 '18 at 8:23