How to calculate time spent on a world grid with circumpolar trips
Sorry to have missed this. There's really no way to do this in longitude / latitude without a lot of manual handling. You could split into two trips, one using the eastern hemisphere, the other the western and use a shared [-180,180] grid that you sum together. What to do at the split coordinate depends on your actual data, there are many options. If you are still looking for an answer I could have a closer look at your example code, but real data would be more helpful if you can share it. There are a number of problems in the code you posted that prevent it from being easily run. Cheers, Mike.
On Fri, Aug 6, 2010 at 8:10 PM, Pinaud David <pinaud at cebc.cnrs.fr> wrote:
Dear all,
I want to calculate a grid of time spent per cell with albatross data. The
goal is then to relate this to fisheries data set on a 1x1? grid.
I use the trip package with tripGrid(), but have some problems with birds
crossing the date line.
Usually I use a projection to do that but here I need to keep the 1x1? grid
in order to relate density tu some data (fisheries...).
I've tried with recenter() and nowrapSpatialLines(), without success.
Here a code:
library(trip)
library(maps) ? # to draw world maps
d <- data.frame(x=c(seq(50, 180, 10), seq(-170, 50, 10)), y=rnorm(n=37,
m=-50, sd=3), tms=Sys.time()-(37:1*60*60*24), id=rep(1, 37)) ?# an example
coordinates(d) <- ~x+y
tr <- trip(d, c("tms", "id"))
proj4string(at) <- CRS("+proj=longlat +ellps=WGS84")
map('world')
plot(tr, add=T)
points(coordinates(tr), t="l")
# grid
grid.base <- expand.grid(Long=seq(-180, 180, by=1), Lat=seq(-70, -15, by=1))
? # tab coords
coordinates(grid.base) <- ~ Long + Lat
proj4string(grid.base) <- CRS("+proj=longlat +ellps=WGS84")
gt <- makeGridTopology(grid.base, cellsize = c(1, 1)) ?# def grid topo
# time spent per cell
tppc <- tripGrid(tr, grid = gt)
tppc01 <- tppc ? ? ? ? ? ? ? # a copy for pres/abs
tppc01 at data$z <- ifelse(tppc at data$z > 0, 1, 0)
image(tppc, col=rev(heat.colors(500))) ?# seems ok
map('world', add=T)
points(coordinates(tr), t="l")
# but with presence/absence:
image(tppc01, col=c("white", "green"))
map('world', add=T)
points(coordinates(tr), t="l", col="blue")
Some ideas ?
Many thanks for your time.
David
--
***************************************************
David PINAUD
Ing?nieur de Recherche "Analyses spatiales"
Centre d'Etudes Biologiques de Chiz? - CNRS UPR1934
79360 Villiers-en-Bois, France poste 485
Tel: +33 (0)5.49.09.35.58
Fax: +33 (0)5.49.09.65.26
http://www.cebc.cnrs.fr/
***************************************************
__________ Information from ESET Mail Security, version of virus signature
database 5346 (20100806) __________
The message was checked by ESET Mail Security.
http://www.eset.com
_______________________________________________ R-sig-ecology mailing list R-sig-ecology at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
Michael Sumner Institute for Marine and Antarctic Studies, University of Tasmania Hobart, Australia e-mail: mdsumner at gmail.com