[R-sig-Geo] Spatial3dArray - coordinates method
Torleif Markussen Lunde
torleif.lunde at cih.uib.no
Thu Nov 19 23:34:37 CET 2009
Hi Mike
At the moment I have written wrapper functions around ncdf. As in get.var.ncdf
you can subset which area to read, and which (continuous) time dimensions you
want to read. At the moment the functions (to correctly return time, lat and
long) are limited to output from the WRF-model (http://www.wrf-model.org/),
but it could easily be modified to other netcdf files as long you know the name
of lat, long and time.
If r-forge accepts my application (delivered yesterday) the methods should
appear there soon as rwrf.
Since ncdf works "out of the box" on Fedora I landed on that one. I also
tested it on windows XP, and no problems there either. It is a couple of years
since I tried RNetCDF, so I should not speak to strongly in favor for any of
them. Of course the bad thing with both is that they have not been updated for
a while. Still my impression is that they are both superior to GDALs NetCDF
support(?).
If the classes proves to be robust and the quality is sufficient it would be
natural to include them in the sp package in the future (instead of having to
dig around to look for sp classes). Currently spplot methods for
Spatial3dArray only support plotting of single times for levelplot,
contourplot, and wireframe (unless animation is requested). I have some
problems with the lattice graphics inside functions (also with print()) that
needs to be sorted out first.
Examples:
first convert Spatial3dArray to data.frame
print(levelplot(z~long*lat | time, data = tmp.df, aspect="iso"))
will cause the last time to over plot all frames inside a function.
What could work is:
for (i in 1:2) {
if (i == 1) {
a <- paste('c(slot(surf, "data")[ , , ', i, ',2])',
sep = "")
} else {
a <- paste(a, ' + c(slot(surf, "data")[ , , ', i, ',2])',
sep = "")
}
}
levelplot(parse(text=a)~
c(slot(surf, "coords")$long) * c(slot(surf, "coords")$lat),
strip = strip.custom(factor.levels=c("time1", "time2")))
where neither parse, print(a, quote = FALSE) or cat works as the z variable.
Probably this approach is not meant to work.
There are quite some issues that has to be solved before the class is
production-ready.
By the way, if someone has an idea on how to solve the last code snippets I
would be more than happy.
Best wishes
Torleif
On Thursday 19 November 2009 20:23:35 Michael Sumner wrote:
> Wow, this is great - I was thinking about this just yesterday.
>
> Torleif: do you have an opinion on which NetCDF path is the most
> useful for R with sp? RNetCDF or ncdf? GDAL is workable but takes
> extra effort to build and then reconstruct 3d/4d from 2d bands. (I use
> Windows mostly)
>
> I use the RNetCDF package a lot, mainly because that's the one I
> learnt to use first - there are binaries for Windows. It has some
> problems in terms of R-style but they could be easily fixed.
>
> Regards, Mike.
>
> On Fri, Nov 20, 2009 at 12:26 AM, Edzer Pebesma
>
> <edzer.pebesma at uni-muenster.de> wrote:
> > No problems; sp in csv now has this, the next release will have it.
> >
> > Torleif Markussen Lunde wrote:
> >> Hi
> >>
> >> To read netcdf data (or any other "gridded" spatial time data) I find it
> >> convenient to define new classes Spatial3dArray and Spatial4dArray.
> >>
> >> setClass("Spatial3dArray",
> >> representation("Spatial", data = "array", coords = "list",
> >> time = "character", btime = "character"),
> >> prototype= list(data = array(NA, c(1,1,1,1)),
> >> bbox=matrix(NA),
> >> proj4string = CRS(as.character(NA)),
> >> coords = list(1,1),
> >> time = "posix",
> >> btime = "posix"))
> >>
> >>
> >> ##########################################
> >> ###################EXAMPLE##################
> >> ##########################################
> >>
> >> x <- matrix(seq(-10, 10, length = 100), 100, 100,
> >> byrow = FALSE)
> >> y <- matrix(seq(-10, 10, length = 100), 100, 100,
> >> byrow = TRUE)
> >>
> >> tm <- 1:10
> >> tm.c <- as.character(seq(as.POSIXct("2002-01-01 06:00:00",
> >> "2002-01-01 15:00:00"),
> >> by="hours",
> >> length.out=10))
> >>
> >> z <- array(NA, c(dim(x)[1], dim(x)[2], length(tm.c), 1))
> >>
> >> for (i in 1:10) {
> >> z[,,i,] <- i * ( sin(sqrt(x^2+y^2)))
> >> }
> >>
> >> sin3dA <- new("Spatial3dArray",
> >> data = z,
> >> coords = list(x, y),
> >> bbox = matrix(c(min(x), min(y), max(x), max(y), 2, 2), 2, 2,
> >> dimnames = list(NULL, c("min","max"))),
> >> time = tm.c,
> >> btime = c(min(tm.c), max(tm.c)))
> >>
> >> dimnames(slot(sin3dA, "data")) = list(NULL,
> >> NULL,
> >> slot(sin3dA, "time"),
> >> c("a"))
> >> names(slot(sin3dA, "coords")) <- c("x", "y")
> >>
> >>
> >> ##########################################
> >>
> >> for the coordinates method I would like to have two options on how to
> >> return the coordinates; "list" or default "sp":
> >>
> >> coordinates.3dArray <- function (obj, type = "sp") {
> >> lat <- slot(obj, "coords")[[1]]
> >> long <- slot(obj, "coords")[[2]]
> >> if (type == "list") {
> >> return(list(long=long, lat=lat))
> >> } else if (type == "sp") {
> >> res <- as.matrix(cbind(c(long), c(lat)))
> >> dimnames(res) <- list(NULL, c("x1", "x2"))
> >> }
> >> }
> >> setMethod("coordinates", signature("Spatial3dArray"),
> >> coordinates.3dArray)
> >>
> >> This means that the default coordinates method in sp has to include the
> >> option "..." . Would it be possible to include this in a future release
> >> of sp?
> >>
> >> The reason I want to keep the list option is to use a matrix oriented
> >> approach in spplot, overlay, etc. methods. I also feel having a
> >> matrix/array approach with these kind of data makes sense. Allowing type
> >> = "sp" means overlay() will work more or less out of the box (however I
> >> would like to return a matrix), and still I could get the list/matrix
> >> when desired.
> >>
> >>
> >> Best wishes
> >> Torleif Markussen Lunde
> >> Centre for International Health
> >> Bjerknes Centre for Climate Research
> >> University of Bergen
> >> Norway
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at stat.math.ethz.ch
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > --
> > Edzer Pebesma
> > Institute for Geoinformatics (ifgi), University of Münster
> > Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> > 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de/
> > http://www.springer.com/978-0-387-78170-9 e.pebesma at wwu.de
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list