raster package. raster() function scales data when reading from HDF
Thanks Matteo,
I dug the problem a bit, and found out where it may come from. There is a as.is argument in the getRasterData function (rgdal) that if set to FALSE applies a scaling factor fetched from the metadata to the data. The argument default to FALSE. When set to TRUE, no scaling is applied.
# Get sds from data previously downloaded
sds = getSds('~/MODIS_test/MOD13A2.A2003065.h10v10.005.2008301134056.hdf')
# Read the data first using defaults
a = readGDAL(sds$SDS4gdal[1])
summary(a)
Data attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-19790000 20790000 43840000 43680000 63740000 99860000 1043545
# With as.is=TRUE
a = readGDAL(sds$SDS4gdal[1], as.is=TRUE)
summary(a)
Data attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
-1979 2079 4384 4368 6374 9986 1043545
I'm not convinced that automatically scaling data is such a good default setting.
Unfortunately the raster() function does not seem to take the as.is argument and therefore uses the default settings of getRasterData(). Any way we could get more control of that?
Cheers,
Lo?c
From: Matteo Mattiuzzi [matteo.mattiuzzi at boku.ac.at]
Sent: Friday, May 03, 2013 12:01 PM
To: R-sig-Geo at stat.math.ethz.ch; Dutrieux, Loic
Subject: Re: [R-sig-Geo] raster package. raster() function scales data when reading from HDF
Sent: Friday, May 03, 2013 12:01 PM
To: R-sig-Geo at stat.math.ethz.ch; Dutrieux, Loic
Subject: Re: [R-sig-Geo] raster package. raster() function scales data when reading from HDF
Dear Loic, I don't know why this happens, but I'm not so sure if it depends directly on your GDAL installation. By running your code I get the same results, but converting the file to geotiff (using internally a direct call to gdal GDAL 1.9.2, released 2012/10/08) the values seam to be correct (except NA that should be -3000 and not 0 but that's not a raster package related problem!). I hope this helps a little bit further... Matteo # downloading and converting > library(MODIS) > runGdal(product="MOD13A2",begin="2003065",end="2003065",tileH=10,tileV=10,SDSstring=1,outDirPath="~/MODIS_test",job="test") pixelSize = asIn resamplingType = near outProj = +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs Output Directory = ~/MODIS_test/test Getting file from: LPDAAC ############################ --2013-05-03 11:33:46-- http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2003.03.06/MOD13A2.A2003065.h10v10.005.2008301134056.hdf Resolving e4ftl01.cr.usgs.gov (e4ftl01.cr.usgs.gov)... 152.61.133.5, 2001:49c8:4000:121d::5 Connecting to e4ftl01.cr.usgs.gov (e4ftl01.cr.usgs.gov)|152.61.133.5|:80... connected. HTTP request sent, awaiting response... 200 OK Length: 7314533 (7.0M) [application/x-hdf] Saving to: `gdalinfo /home/matteo/MODIS_test/test/MOD13A2.A2003065.1_km_16_days_NDVI.tif' 100%[===========================================================================================================================================================================>] 7,314,533 1.77M/s in 5.3s 2013-05-03 11:33:52 (1.31 MB/s) - `/home/matteo/MODIS_ARC/MODIS/MOD13A2.005/2003.03.06/MOD13A2.A2003065.h10v10.005.2008301134056.hdf' saved [7314533/7314533] Downloaded by the first try! Creating output file that is 1200P x 1200L. Processing input file HDF4_EOS:EOS_GRID:/home/matteo/MODIS_ARC/MODIS/MOD13A2.005/2003.03.06/MOD13A2.A2003065.h10v10.005.2008301134056.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI. Using internal nodata values (eg. -3000) for image HDF4_EOS:EOS_GRID:/home/matteo/MODIS_ARC/MODIS/MOD13A2.005/2003.03.06/MOD13A2.A2003065.h10v10.005.2008301134056.hdf:MODIS_Grid_16DAY_1km_VI:1 km 16 days NDVI. 0...10...20...30...40...50...60...70...80...90...100 - done. > a <- raster('~/MODIS_test/test/MOD13A2.A2003065.1_km_16_days_NDVI.tif') > summary(a) MOD13A2.A2003065.1_km_16_days_NDVI Min. -1760 1st Qu. 0 Median 0 3rd Qu. 865 Max. 9965 NA's 0 > sessionInfo() R version 2.15.3 (2013-03-01) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] XML_3.96-1.1 RCurl_1.95-4.1 rgeos_0.2-16 rgdal_0.8-8 MODIS_0.8-08 [6] bitops_1.0-5 raster_2.1-25 sp_1.0-9 loaded via a namespace (and not attached): [1] grid_2.15.3 lattice_0.20-15 tools_2.15.3 >>> "Dutrieux, Loic" <loic.dutrieux at wur.nl> 05/03/13 10:55 AM >>> Dear all, I've noticed some unexpected behavior while loading HDF4 data with the raster() function. See the example below. # In the terminal mkdir ~/MODIS_test cd ~/MODIS_test/ wget ftp://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2003.03.06/MOD13A2.A2003065.h10v10.005.2008301134056.hdf # 7MB # R library(raster) library(MODIS) sds = getSds('~/MODIS_test/MOD13A2.A2003065.h10v10.005.2008301134056.hdf') a = raster(sds$SDS4gdal[1]) summary(a) Returns tMin. -19580000 1st Qu. 20930000 Median 43820000 3rd Qu. 63740000 Max. 99650000 NA's 0 Which is outside of the valid range by a factor of 10000. My first guess is that it has to do with the version of GDAL intalled; this started after upgrading from gdal1.7.3. to gdal1.9.2. Any idea? > sessionInfo() R version 2.15.2 (2012-10-26) Platform: x86_64-pc-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] MODIS_0.6-4 bitops_1.0-5 rgdal_0.8-6 raster_2.1-25 sp_1.0-8 loaded via a namespace (and not attached): [1] grid_2.15.2 lattice_0.20-15 tools_2.15.2 Thanks in advance Best regards, -- Lo?c Dutrieux Laboratory of Geo-Information Science and Remote Sensing Wageningen University The Netherlands _______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo