[R-sig-Geo] raster package. raster() function scales data when reading from HDF
Dutrieux, Loic
loic.dutrieux at wur.nl
Fri May 3 15:19:45 CEST 2013
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
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
More information about the R-sig-Geo
mailing list