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

[Jonathan Greenberg]

Hi Florian:

I took a look at your input, and I SUSPECT what you are seeing is a
problem with the documentation, but I agree you are seeing some odd
things.  Using only gdalUtils:

filename <- "MOD07_L2.A2013001.0835.006.2013001192145.hdf"
gdalinfo(modis_file,sd=8)
# Confirmed that your offset is -15000, scale is: 0.00999999977648258 or ~ 0.01

modis_trans <- gdal_translate(filename, dst_dataset = "tmp.tif",
sd_index = 8,output_Raster=TRUE)
summary(modis_trans)
# Confirmed the same values you got above
plot(modis_trans) # Looks good

modis_trans_conv <- (modis_trans+15000)*0.00999999977648258
summary(modis_trans_conv)
# Yep:
Min.        1.32690
1st Qu.     1.51870
Median      1.54110
3rd Qu.     1.59475
Max.        1.80070
NA's    53529.00000

I also confirmed the 15000 values in QGIS, which is a bit odd...

I'm not convinced this is a file translation issue, or at least one
the R community may be able to help with, since that subdataset is
documented as being a 16-bit integer (not unsigned) (QGIS also
confirmed this).  Also, gdalUtils should be using whatever GDAL you
have installed -- I've got the latest and greatest on my machine and
can replicate this issue, so if the translation is wrong, this is an
issue with the main GDAL developers (gdalUtils is simply wrapping
those functions).

So, I'd recommend a few email barrages:
1) Contact the MODIS folks and CONFIRM that once the gains/offsets are
applied, the output should be Kelvin, and not Celsius (which it
appears to be -- 1.6 deg C = 274.75 K, in the range you were looking
for).  Basically: they may have documented this wrong.
2) Contact the GDAL people and show them the gdal_translate and
gdalinfo results (run gdalUtils "gdal_translate" with the verbose=TRUE
flag and copy the "GDAL command being used:" line, which is the "pure"
gdal_translate command).
3) Contact the QGIS folks and let them know what you are seeing :)

--j


On Thu, Jul 17, 2014 at 2:27 AM, Florian Detsch
<florian.detsch at staff.uni-marburg.de> wrote:
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



-- 
Jonathan A. Greenberg, PhD
Assistant Professor
Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
259 Computing Applications Building, MC-150
605 East Springfield Avenue
Champaign, IL  61820-6371
Phone: 217-300-1924
http://www.geog.illinois.edu/~jgrn/
AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007