Re: [R-sig-Geo] Spatial3dArray - coordinates method
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@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.pebe...@wwu.de
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial3dArray - coordinates method
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.pebe...@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@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.pebe...@wwu.de
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
___
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial3dArray - coordinates method
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.pebe...@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:
