[R-sig-Geo] problem reading hdf image

[Agustin Lobo]

Thanks!
The following works:
r <- raster("HDF4_EOS:EOS_GRID:/media/Iomega_HDD/VOLCAN/Modis_adls_MOD09A1/MOD09A1.A2011289.h11v11.005.2011298101214.hdf:MOD_Grid_500m_Surface_Reflectance:sur_refl_b01")

And much beter, the following works as well:
modislista <- getSds(file.path(modisdir,"MOD09A1.A2011289.h11v11.005.2011298101214.hdf"))
length(modislista$SDS4gdal)
s <- stack(modislista$SDS4gdal)
names(s) <- modislista$SDS4names
inMemory(s)
[1] FALSE
show(s)
class       : RasterStack
dimensions  : 2400, 2400, 5760000, 13  (nrow, ncol, ncell, nlayers)
resolution  : 463.3127, 463.3127  (x, y)
extent      : -7783654, -6671703, -3335852, -2223901  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs
names       : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6,
layer.7, layer.8, layer.9, layer.10, layer.11, layer.12, layer.13
min values  :  -32768,  -32768,  -32768,  -32768,  -32768,  -32768,
 -16,      NA,  -32768,   -32768,   -32768,        0,        0
max values  :   32767,   32767,   32767,   32767,   32767,   32767,
7753,      NA,   32767,    32767,    32767,    65535,    65535

Note the result of inMemory(s), which is critical for all this to be
operational (thanks Robert!!).
I think we should have a wrap function (ie. readMODIS) in package
MODIS doing this. Then the user can just discard those layers in s
that would be of no interest.

Agus

On Tue, Oct 8, 2013 at 2:09 PM, Matteo Mattiuzzi
<matteo.mattiuzzi at boku.ac.at> wrote:
> Section 'Subdatasets'ie:
> SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD09A1.A2011289.h11v11.005.2011298101214.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01
>
>
> # this should work if R is in the HDF directory! if not redo gdalinfo with the file absolute path:  gdalinfo /path/to/MOD09A1.A2011289.h11v11.005.2011298101214.hdf
> r <- raster("HDF4_EOS:EOS_GRID:MOD09A1.A2011289.h11v11.005.2011298101214.hdf:MOD_Grid_500m_Surface_Reflectance:sur_refl_b01")
>
>
> Take care about the quotation marks! You can do:
> r <- raster("HDF4_EOS:EOS_GRID:MOD09A1.A2011289.h11v11.005.2011298101214.hdf:MOD_Grid_500m_Surface_Reflectance:sur_refl_b01") # removed
>
> # or
> r <- raster('HDF4_EOS:EOS_GRID:"MOD09A1.A2011289.h11v11.005.2011298101214.hdf":MOD_Grid_500m_Surface_Reflectance:sur_refl_b01') # not removed
>
>
>
> # The MODIS Package contains the function 'getSds' that extracts the SDS name, may you find something usable in there.
>
>
> I hope this works as I am writing without testing!
> Matteo
>
>

>>>>> Agustin Lobo <alobolistas at gmail.com> 10/07/13 5:42 PM >>>
>> Hi!
>>
>> I get the following error when try to read an hdf image (release 4):
>>
>>> modisima <- stack(file.path(modisdir,"MOD09A1.A2011289.h11v11.005.2011298101214.hdf"))
>> rgdal: version: 0.8-11, (SVN revision 479M)
>> Geospatial Data Abstraction Library extensions to R successfully loaded
>> Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24
>> Path to GDAL shared files: /usr/share/gdal/1.10
>> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
>> Path to PROJ.4 shared files: (autodetected)
>> Error in GDALinfo(filename, silent = silent, returnRAT = RAT,
>> returnCategoryNames = RAT) :
>>   object 'RATlist' not found
>> Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) :
>>   Cannot create a RasterLayer object from this file.
>> In addition: Warning message:
>> In dim(x) : no bands in dataset
>>
>> Also:
>>
>>> GDALinfo(file.path(modisdir,"MOD09A1.A2011289.h11v11.005.2011298101214.hdf"))
>> Error in GDALinfo(file.path(modisdir,
>> "MOD09A1.A2011289.h11v11.005.2011298101214.hdf")) :
>>   object 'scaleOffset' not found
>> In addition: Warning messages:
>> 1: In dim(x) : no bands in dataset
>> 2: In GDALinfo(file.path(modisdir,
>> "MOD09A1.A2011289.h11v11.005.2011298101214.hdf")) :
>>   GeoTransform values not available
>>
>> Nevertheless, hdf4 and 5 are listed by gdalDrivers() and, on the linux
>> terminal, I do get
>> output from gdalinfo (quite long, let me know if I should send it).
>> Also, I can display the
>> image in Qgis with no problem.
>>
>> Any idea on where the problem could be?
>> Thanks,
>> Agus
>>
>>> sessionInfo()
>> R version 3.0.2 (2013-09-25)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>>
>> locale:
>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>> LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8
>>  [5] LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8
>> LC_PAPER=es_ES.UTF-8       LC_NAME=C
>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] rgdal_0.8-11  MODIS_0.10-8  raster_2.1-49 sp_1.0-13
>>
>> loaded via a namespace (and not attached):
>> [1] grid_3.0.2      lattice_0.20-23 tools_3.0.2
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>
>
>
>
>
>