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 -
Here is a zoomed-in view -
Two more things I tried -
gdal_translateto 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.
- Reading the raster with R, using the
gdalUtilspackage, in which case the extent information is not even being read and I see a rectangular image such as this -
- 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 -
library(rgdal) library(raster) library(sf) library(gdalUtils) 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 %>% readGDAL %>% raster lat = sds %>% 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 -