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

Roger Bivand Roger.Bivand at nhh.no
Sun Jan 25 21:35:12 CET 2009


On Sun, 25 Jan 2009, Roger Bivand wrote:

> On Sat, 24 Jan 2009, Sebastian P. Luque wrote:
>
> Hi,
>
> If you have the cell centre offset for the lower left cell, then:
>
> is.na(ice$band1) <- ice$band1 < 0
> vice <- ice$band1
> mice <- matrix(vice, ncol=760, byrow=TRUE)
> mice1 <- t(mice[1120:1,])
> grd <- GridTopology(c(-4941217, -4941217), c(10000,10000), c(760, 1120))

With:

grd <- GridTopology(c(-((760/2)*10000) + 5000, -((1120/2)*10000) + 5000),
   c(10000,10000), c(760, 1120))

that is assuming that the input grid is centred at the pole, and that the 
cell sizes are 10km by 10km, we get much closer, but the whole grid is 
shifted South. Once we establish the metadata by reverse engineering, 
setting up your workflow ought to be feasible.

Roger


> SG <- SpatialGrid(grd, proj4string=CRS("+init=epsg:3411"))
> im <- list(x=sort(coordinatevalues(getGridTopology(SG))$s1),
> y=sort(coordinatevalues(getGridTopology(SG))$s2), z=mice1)
> image(im, asp=1)
> SGDF <- image2Grid(im, p4="+init=epsg:3411")
> image(SGDF, axes=TRUE)
>
>
> gets you there, but here is crucially dependent on knowing the offset (here 
> taken as -90E, 31.08831N (value from netCDF description)). This is wrong, as:
>
> SPixDF <- as(SGDF, "SpatialPixelsDataFrame")
> SPDF <- as(SPixDF, "SpatialPointsDataFrame")
> SPDF_ll <- spTransform(SPDF, CRS("+proj=longlat"))
>
> then:
>
> library(maptools)
> xx <- GE_SpatialGrid(SPDF_ll)
> png("ice.png", width=xx$width, height=xx$height, bg="transparent")
> par(mar=c(0,0,0,0), xaxs="i", yaxs="i")
> plot(SPDF_ll, cex=0.3, col="orange", xlim=xx$xlim, ylim=xx$ylim,
>  setParUsrBB=TRUE)
> dev.off()
> kmlOverlay(xx, "ice.kml", "ice.png")
>
> shows - it puts Greenland well West and a little North of its position when 
> viewed in GE.
>
> Getting there ...
>
> Roger
>

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




More information about the R-sig-Geo mailing list