[R-sig-eco] simm.levy in adehabitat

[Clément Calenge]

Hi Tim,

>> p.sett
>>     
> *********** List of class ltraj ***********
> Type of the traject: Type II (time recorded)
> Regular traject. Time lag between two locs: 900 seconds
> Characteristics of the bursts:
>        id      burst nb.reloc  NAs          date.begin            date.end
> 1   p1801 1988038937        6    1 2006-08-09 09:00:00 2006-08-09 10:15:00
> 2   p1801 1988038938        3    0 2006-08-10 02:00:00 2006-08-10 02:30:00
> 3   p1801 1988038939        9    1 2006-08-11 03:00:00 2006-08-11 05:00:00
> 4   p1801 1988038940       17    5 2006-08-12 02:00:00 2006-08-12 06:00:00
> 5   p1801 1988038941       23   10 2006-08-13 02:00:00 2006-08-13 07:30:00
> ...
> 237 p1801 1988039304      126  110 2007-08-10 01:15:00 2007-08-11 08:30:00
> 238 p1801 1988039305        1    0 2007-08-11 02:45:00 2007-08-11 02:45:00
>
> Ideally what I'd like to do is define the simulation using the
> characteristics of a given 'burst' from my ltraj object. I'd appreciate some
> pointers on the syntax of calling the simulation using the the
> characteristics of each row in p.sett. Sample code follows:
>
>   
>> sim.levy<-simm.levy(date=1:126, id=p.sett$id=="1801",
>>     
> burst=p.sett$burst=="1988039304", x0=c(174.48255,-36.819058), typeII=T)
> Error in as.ltraj(data.frame(co, si), date, id, burst, typeII = typeII) :
>   id should be of the same length as xy, or of length 1
>
> How do I define the arguments of the function using the elements of a given
> burst?
>   


First note that the "burst" and "id" are stored as attributes of each 
element of the list of class ltraj, and not as elements of the list 
themselves. It means that p.sett$burst and p.sett$id do not exist. To 
get the burst or the id of the trajectories of the list, you can simply type

burst(p.sett)
id(p.sett)

See the help page of as.ltraj() for trajectory management with 
adehabitat. Now, it seems that you want to simulate a levy walk at the 
same times as your actual trajectory, beginning at the same point. In 
this case, you should use the following code:

res <- lapply(1:length(p.sett), function(i) {   
           simm.levy(date = p.sett[[i]]$date, id = id(p.sett)[i],
                            x0 = c(p.sett[[i]]$x[1], p.sett[[i]]$y[1]),
                            burst = burst(p.sett)[i], typeII = T, mu=2)
})
result <- do.call("c", res)
plot(result)

But note that you may have to change the parameter mu, depending on what 
you want to do. See the help page of simm.levy.
BTW, there is a bug in the function concerning the parameter x0. If you 
want to use parameter x0 (required only if you want to match the 
trajectory with maps), you should use the following fixed function, 
which will be included in the next version of adehabitat:

simm.levy <- function (date = 1:500, mu = 2, l0 = 1, x0 = c(0, 0), id = 
"A1",
    burst = id, typeII = TRUE)
{
    if (typeII)
        class(date) <- c("POSIX", "POSIXct")
    n <- length(date)
    dt <- c(diff(unclass(date)))
    if (all(dt - dt[1] > 1e-07))
        stop("the time lag between relocations should be constant")
    ang <- runif(n - 2, -pi, pi)
    v = dt * (l0 * (runif(n - 1)^(1/(1 - mu))))
    ang = cumsum(c(runif(1, 0, 2 * pi), ang))
    si = c(x0[2], x0[2]+cumsum(v * sin(ang)))
    co = c(x0[1], x0[1]+cumsum(v * cos(ang)))
    res <- as.ltraj(data.frame(co, si), date, id, burst, typeII = typeII)
    return(res)
}


Hope this helps,


Clément.

-- 
Clément CALENGE
Office national de la chasse et de la faune sauvage
Saint Benoist - 78610 Auffargis
tel. (33) 01.30.46.54.14