Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal

Question

I've been trying to accurately import and process some .hdf files from the MODIS Atmospheric Profile product (MOD07_L2) in R. There's something going wrong during data import. The below code can be reproduced using one example .hdf file (MOD07_L2.A2013001.0835.006.2013001192145.hdf) that can be downloaded from Dropbox.

library(rgdal)

# Extraction of metadata via `GDALinfo`    
filename <- "MOD07_L2.A2013001.0835.006.2013001192145.hdf"
gdalinfo <- GDALinfo(filename, returnScaleOffset = FALSE)
metadata <- attr(gdalinfo, "subdsmdata")

# Extraction of SDS string for parameter 'Skin_Temperature' (formerly 'Surface_Temperature')    
sds <- metadata[grep("Skin_Temperature", metadata)[1]]
sds <- sapply(strsplit(sds, "="), "[[", 2)

# Raster import via `readGDAL`   
sds.rg <- readGDAL(sds)

So far, so good, but here comes the confusing part:

> summary(sds.rg$band1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 -14870  -14850  -14850  -14840  -14840  -14820   53529 

Considering the fact that Skin_Temperature has an official valid range from 150 to 350 K (see MOD07_L2:Format & Content), the mean would inherit a value of

> (-14840 - (-15000)) * 0.01
[1] 1.6

after considering the corresponding add_offset (-15000) and scale_factor (0.01). Note that we're still talking about Kelvin, not °C. Extracting SDS No. 8, i.e. Skin_Temperature, using

library(gdalUtils)
gdal_translate(filename, dst_dataset = "tmp.tif", sd_index = 8)

and opening the resulting file called "tmp.tif" in QGIS resulted in seemingly reliable values centered around 15000, i.e. roughly 27 °C. However, importing "tmp.tif" back into R using raster again resulted in values comparable to the ones shown above:

> summary(raster("tmp.tif"))
              tmp
Min.    -14867.31
1st Qu. -14848.13
Median  -14845.89
3rd Qu. -14840.53
Max.    -14819.93
NA's         0.00

I've been searching the internet and stumbled across similar problems related to rgdal. However, when I tried to cast toUnSigned on band 1 of my previously generated 'SpatialGridDataFrame', I received the following error message:

> toUnSigned(sds.rg$band1, 16)
Error in toUnSigned(sds.rg$band1, 16) : band not integer

Apparently, the data imported into R is not even of type integer (what it is supposed to be), but numeric:

> sds.rg$band1[1:5]
[1]        NA        NA -14839.40 -14840.25 -14839.26

Is there an apparent mistake in my code, or is there any point I miss when importing the .hdf and .tif files using rgdal?

Answer 1

The answer is surprisingly simple.

sgr_lst <- readGDAL(sds, as.is = TRUE)

solves the issue. The only thing left to do then is transform the resulting SpatialGridDataFrame to a proper Raster* object using raster(). For this purpose, it is necessary to retrieve the west, east, south, and north bounding coordinates (see ?extent: xmin, xmax, ymin, ymax) from the metadata e.g. via

meta <- attr(gdalinfo, "mdata")

## string patterns
crd_str <- paste0(c("WEST", "EAST", "SOUTH", "NORTH"), "BOUNDINGCOORDINATE")
## search for patterns in metadata
crd_id <- sapply(crd_str, function(i) grep(i, meta))
## extract information
crd <- meta[crd_id]
## create 'extent' object
crd <- as.numeric(sapply(strsplit(crd, "="), "[[", 2))
ext <- extent(crd)

The thus determined extent (together with the EPSG code) can then be passed on to the finally created RasterLayer which, after applying the scale factor and offset, definitely looks fine now.

## rasterize, apply offset and scale factor
rst_lst <- (raster(sgr_lst) - -1.5e+04) * 1.0e-02
## set extent and coordinate reference system
extent(rst_lst) <- crd
projection(rst_lst) <- "+init=epsg:4326"

Answer 2

I have not managed to do it in R either.. and I have also spent countless hours. What I do now is this:

1. Use the ModisReprojectionTool to extract the layers and the
   subsets I need as binary files
2. read the binary files and if necessary convert them to "raster" objects. I mainly use them as matrix and in the end transform them into raster-objects to write them as TIFF

the code to run MRT:

# Reproject from HDF to plain binary with MRT
setwd('/Volumes/Archive1/MOD11C3V5/')
DIR <- getwd()
# Run the Modis Reprojection Tool once with one HDF-file and save the desired parameterization in a .prm-file
# Use this .prm-file here:
ReprojectionParamter <- 'Pamir0.05_binray.prm'

# Getting the file list that you want to process
FileList <- list.files()
FileListHDF <- FileList[which(regexpr(pattern='hdf$',FileList)>0)] # only the hdf-files but not the hdf.xml files

### Setting environmental variables for MRT_DATA_DIR
Sys.setenv(MRT_DATA_DIR='/Volumes/DATA/ModisReprojectionTool/data')
Sys.setenv(MRTDATADIR='/Volumes/DATA/ModisReprojectionTool/data')

for (i in FileListHDF){
  system(command=paste('/Volumes/DATA/ModisReprojectionTool/bin/resample -p ',ReprojectionParamter,' -i ',DIR,'/',i,' -o ',DIR,'/binary_small/',i,'.hdr',sep=''),wait=T,)
}

Then you have to set the "what" argument of readBin to int, numeric, etc... I always have the same extent of my files (defined in the "prm"-file) I get the extent and resolution directly from there:

ReprojectionParamter <- scan('Pamir0.05_binary2.prm', nmax=90,what='character')
SpatExtent.minLon <- as.numeric(ReprojectionParamter[31])
SpatExtent.minLat <- as.numeric(ReprojectionParamter[36])
SpatExtent.maxLon <- as.numeric(ReprojectionParamter[37])
SpatExtent.maxLat <- as.numeric(ReprojectionParamter[30])
SpatExtent.RES <- as.numeric(ReprojectionParamter[72])
# readBinary function
UInt8_LST <- function(f, ...) readBin(f, what = "integer", signed = FALSE, endian = "little", size = 2, ...)
# Finally read the data; Nlat,Nlon etc you can calculate easily from the info in the prm-file; 0.02 was my scale factor
NC.LST_night <- matrix(UInt8_LST(f=LST_night.filename, n= Nlon*Nlat)),nrow=Nlat,ncol=Nlon,byrow=T)*0.02-273.15

Maybe it helps.

Answer 3

Just want to say that fdetsch's answer may not be correct because the extent created from lonmin, lonmax, latmin, latmax may not match with the actual grid boundary. The image is usually tilted.

The dropbox link is invalid so I download the same date file (updated version) from here

Below is what I plotted with fdetsch's code:

Below is my code and plot:

hdfname="MOD07_L2.A2013001.0835.061.2017296134818.hdf"
GDALinfo(hdfname, returnScaleOffset = FALSE) 
#note that the 29th and 30th subdatasets are lat and lon
one <- getSds(hdfname)    
lat=raster(readGDAL(one$SDS4gdal[29], as.is = TRUE))
lon=raster(readGDAL(one$SDS4gdal[30], as.is = TRUE))
Skin_Temperature <- brick(readGDAL(one$SDS4gdal[8], as.is = TRUE))
Skin_Temperature <- (Skin_Temperature - (-15000))*0.009999999776482582 #It seems that the scaling factor changed after data updates
color <- rev(terrain.colors(100))[as.numeric(cut(Skin_Temperature[],breaks = 100))]
plot(lon[],lat[],pch=20,col = color,xlab="lon",ylab="lat") #my plot