[R-sig-Geo] raster package. raster() function scales data when reading from HDF
Matteo Mattiuzzi
matteo.mattiuzzi at boku.ac.at
Fri May 3 12:01:55 CEST 2013
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