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

Florian Detsch florian.detsch at staff.uni-marburg.de
Fri Jul 18 10:21:47 CEST 2014 

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
>

More information about the R-sig-Geo mailing list