Dear Leila,
This is indeed a bug. I have corrected it and I will upload the
correction to CRAN asap. It will be available in a few days. Meanwhile,
sourcing the attached file (i.e. typing source("kernelbb.R")) will allow
to use the function in this context.
Best,
Clément
On 09/19/2012 08:38 AM, Leila Brook wrote:
I'm sorry for cross postings, as I'm not sure where exactly to ask
this question:
I am trying to estimate a Brownian Bridge kernel home range using
adehabitatHR. I created
I created a traj file as below:
bobtraj<-as.ltraj(xy=bob_locs_ss[c("x","y")], id=bob_locs_ss$id,
date=bob_locs_ss$newdate, burst=bob_locs_ss$burst2)
Then tried to create a kernelBB as follows, with a reference grid
around the animal's home range:
bob_bbss<- kernelbb(bobtraj, sig1=4.0811, sig2=50, grid = refgrid,
byburst = TRUE, extent = 0.2)
I then got this error:
Error in as.vector(gr) : no method for coercing this S4 class to a vector
I didn't get the error and the function worked when I ran kernelBB()
without byburst=TRUE. However, I would like to keep this specification
in the model, as my data were collected at high-frequency with two
days on and two days off, and without the burst=TRUE, the kernelBB95
looks quite oversmoothed.
Can anyone provide any advice on how to fix this error?
Thanks for your help,
Leila Brook
--
Clément CALENGE
Cellule d'appui à l'analyse de données
Direction des Etudes et de la Recherche
Office national de la chasse et de la faune sauvage
Saint Benoist - 78610 Auffargis
tel. (33) 01.30.46.54.14
"kernelbb" <- function(ltr, sig1, sig2, grid = 40,
same4all=FALSE, byburst=FALSE, extent = 0.5,
nalpha=25)
{
## verifications
if (!inherits(ltr, "ltraj"))
stop("tr should be of class \"ltraj\"")
x <- do.call("rbind",ltr)[,c("x","y")]
bu <- burst(ltr)
x$burst <- unlist(lapply(1:length(bu), function(i) {
rep(bu[i], nrow(ltr[[i]]))
}))
idd <- id(ltr)
x$id <- unlist(lapply(1:length(idd), function(i) {
rep(idd[i], nrow(ltr[[i]]))
}))
x$date <- unlist(lapply(ltr, function(y) y$date))
## Bases
sorties <- list()
gr <- grid
x <- x[!is.na(x$x),]
x <- x[!is.na(x$y),]
xy<-x[,c("x","y")]
sig12<-sig1^2
sig22<-sig2^2
h<-c(sig1, sig2)
names(h)<-c("sig1","sig2")
fac<-x$burst
if (!byburst)
fac<-x$id
fac<-factor(fac)
lixy<-split(x,fac)
if (is.list(grid)) {
if (is.null(names(grid)))
stop("when grid is a list, it should have named elements")
nn <- names(grid)
lev <- levels(fac)
if (length(lev) != length(nn))
stop("the length of the grid list should be equal to the number of
levels of id/burst")
if (!all(lev%in%nn))
stop("some levels of id/bursts do not have corresponding grids")
}
if (same4all) {
if (inherits(grid, "SpatialPoints"))
stop("when same4all is TRUE, grid should be a number")
grid <- adehabitatHR:::.makegridUD(xy, grid, extent)
}
for (i in 1:length(lixy)) {
dft<-lixy[[i]]
df<-dft[,c("x","y")]
if (!is.list(gr)) {
if (!inherits(gr, "SpatialPoints")) {
if (length(as.vector(gr)) == 1) {
if (!is.numeric(gr))
stop("non convenient grid")
if (!same4all) {
grid <- adehabitatHR:::.makegridUD(xy, gr, extent)
}
}
}
} else {
grid <- gr[[names(lixy)[i]]]
}
gridded(grid) <- TRUE
fullgrid(grid) <- TRUE
pfs <- proj4string(grid)
grrw <- gridparameters(grid)
if (nrow(grrw) > 2)
stop("grid should be defined in two dimensions")
if (is.loaded("adehabitatMA")) {
opteps <- adehabitatMA::adeoptions()$epsilon
} else {
opteps <- 1e-08
}
if ((grrw[1, 2] - grrw[2, 2])> opteps)
stop("the cellsize should be the same in x and y directions")
xyg <- coordinates(grid)
xg <- unique(xyg[,1])
yg <- unique(xyg[,2])
date<-as.double(dft$date)-min(as.double(dft$date))
toto<-.C("kernelbb", double(length(xg)*length(yg)), as.double(xg),
as.double(yg), as.integer(length(yg)),as.integer(length(xg)),
as.integer(nrow(df)), as.double(sig12), as.double (sig22),
as.double(df$x), as.double(df$y), as.double(date),
as.integer(nalpha),
PACKAGE="adehabitatHR")
too <- c(t(matrix(toto[[1]], nrow=length(xg), byrow=TRUE)))
UD <- data.frame(ud=too)
coordinates(UD) <- expand.grid(yg,xg)[,2:1]
gridded(UD) <- TRUE
UD <- new("estUD", UD)
slot(UD, "h") <- list(values=c(sig12=sig12, sig22=sig22),
meth="BB-specified")
slot(UD, "vol") <- FALSE
if (!is.na(pfs))
proj4string(UD) <- CRS(pfs)
sorties[[names(lixy)[i]]] <- UD
}
class(sorties) <- "estUDm"
if (length(sorties)==1)
sorties <- sorties[[1]]
return(sorties)
}
_______________________________________________
AniMov mailing list
AniMov@faunalia.it
http://lists.faunalia.it/cgi-bin/mailman/listinfo/animov
