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

[Roger Bivand]

On Fri, 16 Jan 2009, Sebastian P. Luque wrote:

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

Sorry for the delay. Having now set up a reprodicable route (OSGeo4W GDAL 
and drivers, rgdal built against OSGeo4W as in a posting a couple of days 
ago), I can see the problem, but have no resolution.

Doing gdalinfo on the data itself, I see geographical coordinates, not 
projected. GDAL does not retrieve the grid metadata correctly, and assigns 
default 1x1 grid cells starting at c(0.5, 0.5). Using gdal_translate to 
convert it to a GTiff, I see that I have to set p4s to NA to read at all, 
and the problems remain. I would be interested in knowing who or what 
wrote the HDF5 file, because it doesn't seem to be what it says it is (no 
grid metadata, different CRS from your assumption), and so the data may 
actually be in the wrong order too.

Roger

PS - I can't see that you've asked on the GDAL list, hence the reply here.

>
> ---<--------------------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,
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no