Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal
The solution using readGDAL(..., as.is = TRUE) worked out fine and all SDS layers are in a valid range now. Thank you, Lo?c and Jonathan, for the precious input. Florian
On 17.07.2014 18:49, Dutrieux, Loic wrote:
Hi Florian, Jonathan, I remember experiencing similar issues a while back with MODIS (I think NDVI product), simply because the product metadata were not correct (scaling factor) and readGDAL() uses them by default for automatic scaling on import. The trick was to import using as.is=TRUE and doing the scaling myself using values from the product description. So in your case, that would mean: # Raster import via `readGDAL` sds.rg <- readGDAL(sds, as.is=TRUE) You can also wrap that into raster if you prefer working with classes of the raster package. r0 <- raster(readGDAL(sds, as.is=TRUE)) r <- (r0 - 15000) * 0.01 Hope this helps, Best regards, Lo?c
________________________________________ From: r-sig-geo-bounces at r-project.org [r-sig-geo-bounces at r-project.org] on behalf of Jonathan Greenberg [jgrn at illinois.edu] Sent: Thursday, July 17, 2014 4:54 PM To: Florian Detsch Cc: r-sig-geo at r-project.org Subject: Re: [R-sig-Geo] Cannot properly import MODIS MOD07_L2 .hdf files into R using rgdal 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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo