[AniMov] Re: kNNCH code questions?

Clément Calenge calenge at biomserv.univ-lyon1.fr
Wed Jan 19 18:14:33 CET 2005


Dear Sander,

I am now preparing the next version of adehabitat, which will include 
several new functions,
among them the kNNCH method. So I thank you for your message which points 
out some
gaps in the package.

>x <- rnorm(100)
>y <- rnorm(100)
>
>XYnnch <- data.frame(X=x, Y=y)
>
>aa1<-NNCH(XYnnch, k=5)
>aa2<-NNCH(XYnnch, k=10)
>aa3<-NNCH(XYnnch, k=15)
>print(aa1)
>NNCH.area(aa1)
>par(mfrow=c(2,2)) # a 2 by 2 matrix of plots
>plot(x,y)
>plot(aa1, add.points=FALSE)
>plot(aa2, add.points=FALSE)
>plot(aa3, add.points=FALSE)
>
>Problems, suggestions and questions:
>1) I would like to get insight into the data by plotting several plots in 
>a single graph, but it won't work (had the same problem with plots from 
>hexbin package as well). Maybe I am missing something basic about R plotting!?

Thanks for this suggestion. Actually, the function NNCH estimates one home 
range per animal,
and stores the results in a list of class "NNCH". The function plot.NNCH 
plots several graphs on
the same window, i.e. one plot per animal (as for the results of kernelUD). 
To do that, the function
sets the graphical parameter mfrow (number of rows and columns of graphs on 
the window).
When you have only one animal, the function plot.NNCH therefore sets the 
window so that there
is only one graph, and this cannot be changed.

However, the function can be changed. To proceed, first edit the function:

plot.NNCH <- edit(plot.NNCH)

Then, delete the lines of the function plot.NNCH:

opar<-par(mfrow=n2mfrow(length(x)))
on.exit(par(opar))

(third and fourth line of the function)
And finally, using your example:

x <- rnorm(100)
y <- rnorm(100)

XYnnch <- data.frame(X=x, Y=y)

aa1<-NNCH(XYnnch, k=5)
aa2<-NNCH(XYnnch, k=10)
aa3<-NNCH(XYnnch, k=15)
print(aa1)
NNCH.area(aa1)
par(mfrow=c(2,2)) # a 2 by 2 matrix of plots
plot(x,y)
plot(aa1, add.points=FALSE)
plot(aa2, add.points=FALSE)
plot(aa3, add.points=FALSE)

it works fine.
But be carefull that each object aa contains only one home range !
Otherwise, it won't work. I will think to a better way to deal with such
objects, and add it in the package.

>2) The method is indeed rather slow (for big data sets), but I can just 
>make myself some coffee.
>3) I would like to save the kNNCH polygons to a GIS file. Are the polygons 
>available somewhere? Preferably I would be able to save individual layers 
>for each k, or the polygons would have the k level as one of their 
>attributes. This would allow for powerful GIS analysis later on.

The function is still very slow for big data sets. It would be quicker if
it was programmed in C, interfaced with an R function (as with the function
kernelUD, by example). However, it would be very, very long to program
this method in C, because the home range may contain holes. The R functions
of the package gpclib may deal with holes, but it is much harder to do in C.

I have written a new function to get the vertices of the home range at a 
given level, and
to store it in an object (and a plot function to display it). Because the 
home range may
contain holes, the function stores the objects in a list of class 
"gpc.poly" (the vector
format "area" is not intended to store vector shapes with holes). Export 
such polygons
toward a GIS is much harder than with "area". Maybe by using the functions 
from other
packages... The package gpclib also contains a function "write.polyfile", 
to export the
object, but as stated in the help page of the function: "write.polyfile 
does not return
anything useful". Maybe contacting Roger D. Peng - the maintainer of the
package - may help...

Another way would be to write a function that converts objects of class 
"gpc.poly"
(library gpclib) to lists of class "shapefile" (library shapefiles). It 
would be a long work...
The exportation of objects of class area - a quite simple class - toward 
GRASS has
posed many problems to users. Can you imagine the problems posed with the 
format
"gpc.poly"?

Note that the package will include a new function, NNCH.rast, which 
rasterize the home range
at a given percentage level. The result is a raster map of class "asc", 
which can be
exported toward a GIS (see export.asc). After that, in the GIS, a 
vectorization may
be done... Or the vectorization can be done in R, using the function 
getcontour, but
this function does not take into account the holes inside the shapes...
Still unsolved problem, sorry...

>4) I assume that the NNCH.area function puts out various percentages of 
>the total area covered by all polygons in the kNNCH. Would it be possible 
>to provide the area covered by each individual level of k? I think there 
>might be some ecological relevance in the the way area increases with 
>increasing k! Following Figure 3 in Getz et al 2004 Ecography 27:489.
>Alternatively I could write a loop for each k and store the total area, 
>but than I recalculate all the previous k's each time as well. With the 
>method being slow, that is quite a penalty.

Thanks again for this suggestion. This function is indeed necessary.
Maybe I should take inspiration from the Arcview extension Locoh
options for these functions (many thanks to Wayne for this clue !)...
Note that, at first sight, it seems that a loop will be unavoidable...

>5) I also think there will be ecological relevance in the changes in home 
>range size through time. Thus I have been trying to run the NNCH analysis 
>in a moving window function. The aim is to end up with a dataframe with 
>the window ID, the associated NNCH area, and the mean 'date' of the moving 
>window, so I can plot NNCH area against time. However, hacking other 
>people's code is risky if you have no clue what you are hacking. Can you 
>help me out with this (maybe you know of a better moving window function. 
>I came across this one on the R list):
>
>movingWindowNNHC <- function(data, width, k, units) {
>    mCall = match.call()
>    mCall$width = NULL
>    mCall[[1]] = as.name("lm")
>    mCall$x = mCall$y = TRUE
>    bigfit = eval(mCall, parent.frame())
>    ncoef = length(coef(bigfit))
>    nr = nrow(data)
>    width = as.integer(width)[1]
>    stopifnot(width >= ncoef, width <= nr)
>    y = bigfit$y
>    x = bigfit$x
>    terms = bigfit$terms
>    inds = embed(seq(nr), width)[, rev(seq(width))]
>    multiNNCH <- lapply(seq(nrow(inds)),
>                     function(st) {
>                         aa1<-NNCH(data, k=k, units=units)
>                         areaNNCH <- NNCH.area(aa1, percent=100)
>                         areaNNCH[,1]
>
>                     })
>}
># parameters: data, width of window, k, units
>test <- movingWindowNNHC(XYnnch, 5, 5, "m")
>multiNNCH
>
>I'm sure I will come up with more later on, but this is it for now!

Indeed, a moving window may be useful to analyse the changes in home-range 
size
(it would be also useful for other home-range estimation methods). However, 
you use a "call"
to the method "lm" (linear model) which is not necessary here (line 1 to 13 
of the program).
Moreover, the core of your function

>multiNNCH <- lapply(seq(nrow(inds)),
>                     function(st) {
>                         aa1<-NNCH(data, k=k, units=units)
>                         areaNNCH <- NNCH.area(aa1, percent=100)
>                         areaNNCH[,1]
>
>                     })

applies the function to the same dataset (data) at each iteration.
Try the following, with the same arguments (plus one argument "percent",
indicating the level of the estimation) :

movingWindowNNHC <- function(data, width, k, units, percent) {
         multiNNCH <- lapply(seq(nrow(data)-width),
                        function(st) {
                              aa1<-NNCH(data[st:(st+width),], k=k, unin=units)
                              areaNNCH <- NNCH.area(aa1, percent=percent)
                              return(areaNNCH[,1])
                    })
         return(unlist(multiNNCH))
}

test <- movingWindowNNHC(XYnnch, 20, 5, "m", 100)
plot(test)

However, be patient, because this function is very long (63 seconds on my
computer, with 100 relocations, and 20 neighbours).
Note that this function requires that:
- the relocations are sorted in the data frame data
- the time-intervals between successive relocations are of equal length.
In case of non-equal time interval in the data, this function could not be
used because:
- a vector indicating the date (class POSIXt) needs to be included
   in the function, which will become more complex.
- as the estimate may be affected by the number of relocations used,
   some form of weighting is needed.
Some work is required, when unequal intervals are under study (deep
studies of the properties of the methods, of the validity of the "moving
window" approach, etc)

Finally, note that the next version of the package adehabitat will be 
available
on CRAN in a couple of weeks. I can send a copy of the version in development
of the package (not "clean", but compiling without problems) on request, 
but i do
not attach it to the present message, as the file is rather large (2 MB).

Hope this helps,
Best wishes


Clément.

======================================
UMR CNRS 5558 - Equipe "Ecologie Statistique"
Laboratoire de Biométrie et Biologie Evolutive
Université Claude Bernard Lyon 1
43, Boulevard du 11 novembre 1918
69622 Villeurbanne Cedex
FRANCE
tel. (+33) 04.72.43.27.57
fax. (+33) 04.72.43.13.88
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.faunalia.com/pipermail/animov/attachments/20050119/1f6f02b7/attachment-0001.html


More information about the AniMov mailing list