[R-sig-Geo] readGDAL() and HDF5 files

[Sebastian P. Luque]

Hi,

I'm trying to use HDF5 files with this structure:

---<--------------------cut here---------------start------------------->---
$ gdalinfo conc_200901011200.hdf 
Driver: HDF5/Hierarchical Data Format Release 5
Files: conc_200901011200.hdf
Size is 512, 512
Coordinate System is `'
Subdatasets:
  SUBDATASET_0_NAME=HDF5:"conc_200901011200.hdf"://Data/data[00]
  SUBDATASET_0_DESC=[760x1120] //Data/data[00] (32-bit floating-point)
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)
---<--------------------cut here---------------end--------------------->---

The file (467 kb placed in
http://members.shaw.ca/sluque/conc_200901011200.hdf) is a grid of a
single variable (ice concentration) in a well defined polar
stereographic projection, that loads ok in R, via rgdal's readGDAL():

---<--------------------cut here---------------start------------------->---
R> p4s <- "+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45"
R> ice <- readGDAL("HDF5:\"conc_200901011200.hdf\"://Data/data[00]", p4s=p4s)
HDF5:"conc_200901011200.hdf"://Data/data[00] has GDAL driver HDF5Image 
and has 760 rows and 1120 columns
R> summary(ice)
Object of class SpatialGridDataFrame
Coordinates:
  min  max
x   0 1120
y   0  760
Is projected: TRUE 
proj4string :
[+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x               0.5        1      1120
y               0.5        1       760
Data attributes:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -32800     -99     -99    -117       0     100 
R> ## We are concerned only with the 0-100 range; other values are codes with
R> ## various meanings, so we remove them
R> ice at data[ice at data[[1]] < 0 | ice at data[[1]] > 100, ] <- NA
R> summary(ice)
Object of class SpatialGridDataFrame
Coordinates:
  min  max
x   0 1120
y   0  760
Is projected: TRUE 
proj4string :
[+proj=stere +a=6378273 +b=6356889.44891 +lat_0=90 +lat_ts=70 +lon_0=-45]
Number of points: 2
Grid attributes:
  cellcentre.offset cellsize cells.dim
x               0.5        1      1120
y               0.5        1       760
Data attributes:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
     0.0      0.0      0.0     23.3     36.7    100.0 478521.0 
## but something is wrong with visualizing these data
pdf("~/tmp/ice.pdf")
spplot(ice)
dev.off()
R> sessionInfo()
R version 2.8.1 (2008-12-22) 
x86_64-pc-linux-gnu 

locale:
LC_CTYPE=en_CA.UTF-8;LC_NUMERIC=C;LC_TIME=en_CA.UTF-8;LC_COLLATE=en_CA.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_CA.UTF-8;LC_PAPER=en_CA.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_CA.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_0.5-35    sp_0.9-29       slmisc_0.7.0    lattice_0.17-20

loaded via a namespace (and not attached):
[1] grid_2.8.1
---<--------------------cut here---------------end--------------------->---

I'm not attaching the resulting file because it's too large.  However,
the plot should look like http://members.shaw.ca/sluque/full.ps I'd
appreciate any feedback on what the problem might be.  Thanks.


Cheers,

-- 
Seb