[R-sig-Geo] Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal

[Florian Detsch]

Dear list,

I've been trying to import and process some .hdf files from the MODIS 
Atmospheric Profile product (MOD07_L2) in R for several days now. Still, 
there's something going wrong during data import. The below code can be 
reproduced using an exemplary .hdf file 
(MOD07_L2.A2013001.0835.006.2013001192145.hdf) that can be downloaded 
from Dropbox via the following link: 
https://www.dropbox.com/s/nz6p2011od0bmn4/MOD07_L2.A2013001.0835.006.2013001192145.hdf.

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 e.g. 
http://modis-atmos.gsfc.nasa.gov/MOD07_L2/format.html), 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', from the .hdf file 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 a similar problem 
on R-sig-Geo related to rgdal and signed/unsigned raster data import 
(http://grokbase.com/t/r/r-sig-geo/098mmnbt41/rgdal-and-unsigned-integers). 
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 as shown by the first 5 cell values 
of band 1:

 > 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`? I would be 
extremely grateful for any kind of help.

Cheers,
Florian

-- 
Florian Detsch, M.Sc. Physical Geography
Environmental Informatics
Faculty of Geography
Philipp University of Marburg
Deutschhausstraße 12
35037 Marburg, Germany
Tel. +49 6421 2825323