[AniMov] Re: AdeHabitat questions and suggestions!?

[Clément Calenge]

Skipped content of type multipart/alternative-------------- next part --------------
## fit the nearest-neighbor convex hull

NNCH<-function(xy, id=NULL, k=10, unin = c("m", "km"),
               unout = c("ha", "km2", "m2"))
  {
    xy <- xy[!is.na(xy[, 1]), ]
    xy <- xy[!is.na(xy[, 2]), ]
    unin <- match.arg(unin)
    unout <- match.arg(unout)
    class(xy[,1])<-"double"
    class(xy[,2])<-"double"
    if (!require(gpclib)) 
      stop("package gpclib required")
    if (is.null(id))
      id<-rep(1,nrow(xy))
    id<-factor(id)
    res<-list()
    opar<-par(mfrow=n2mfrow(nlevels(id)))
    on.exit(par(opar))
    
    for (kk in 1:nlevels(id)) {
      xyt<-xy[id==levels(id)[kk],]
      
      li<-list()
      li2<-list()
      lin<-list()
      lin2<-list()
      ar<-0
      
      dij<-as.matrix(dist(xyt))
      idt<-1:nrow(dij)
      
      for (i in 1:nrow(xyt)) {
        iid<-idt[order(dij[i,])][1:k]
        xytmp<-xyt[iid,]
        ch<-chull(xytmp[,1], xytmp[,2])
        li[[i]]<-as(xytmp[ch,], "gpc.poly")
        lin[[i]]<-iid
      }
    
      aa<-unlist(lapply(li, area.poly))
      li<-li[order(aa)]
      lin<-lin[order(aa)]
      idbis<-idt[order(aa)]
      li2[[1]]<-li[[1]]
      lin2[[1]]<-lin[[1]]
      
      for (i in 2:length(li)) {
        li2[[i]]<-union(li2[[i-1]], li[[i]])
        lin2[[i]]<-unique(c(lin2[[i-1]], lin[[i]]))
      }
      
      n<-unlist(lapply(lin2, length))/nrow(xyt)
      ar<-unlist(lapply(li2, area.poly))
      
      if (unin == "m") {
        if (unout == "ha") 
          ar <- ar/10000
        if (unout == "km2") 
          ar <- ar/1e+06
      }
      if (unin == "km") {
        if (unout == "ha") 
          ar <- ar * 100
        if (unout == "m2") 
          ar <- ar * 1e+06
      }

      plot(xyt, ty="n", asp=1, main=levels(id)[kk])
      gr<-grey(n)
      for (i in length(li2):1)
        plot(li2[[i]], poly.args=list(col=gr[i], border=0), add=T)
      
      names(li2)<-round(n*100)
      res[[levels(id)[kk]]]<-list(area=data.frame(levels=round(n*100,2),
                                    area=ar), polygons=li2, xy=xyt)

    }
    attr(res, "units") <- unout
    class(res)<-"NNCH"
    return(res)
  }

## print method

print.NNCH<-function(x)
  {
    cat("***********************************************\n")
    cat("***\n")
    cat("***      Nearest-neighbor convex hull\n\n")
    cat(paste("Home range available for", length(x), "animals:\n"))
    print(names(x))
    cat("\n\nEach animal is a component of the object. For each animal,")
    cat("\nthe following information is available:\n")
    cat("\n$area:       home-range size estimated at various levels")
    cat("\n$polygons:   objects of class \"gpc.poly\" storing the home-range limits")
    cat("\n$xy:         the relocations\n\n")
  }



## compute the home-range size
NNCH.area<-function(x, percent=c(95,90,80,70,60,50,40,30,20,10))
  {
    if (!inherits(x, "NNCH"))
      stop("x should be of class \"NNCH\"")

    res<-matrix(0,nrow=length(percent), ncol=length(x))

    for (kk in 1:length(x)) {
      for (i in 1:length(percent))
        res[i,kk]<-x[[kk]]$ar[max(which(x[[kk]]$area$levels<=percent[i])),2]
    }
    res<-as.data.frame(res)
    row.names(res)<-percent
    names(res)<-names(x)
    class(res) <- c("hrsize", "data.frame")
    attr(res, "units") <- attr(x, "units")
    return(res)
  }


## display the home ranges

plot.NNCH<-function(x, add.points=TRUE, pch=21, bg="white", col="black",
                    cex=0.7, add=FALSE, same4all=TRUE,...)
  {
    if (!inherits(x, "NNCH"))
      stop("x should be of class \"NNCH\"")
    opar<-par(mfrow=n2mfrow(length(x)))
    on.exit(par(opar))
    if (same4all) {
      xxx<-do.call("rbind", lapply(x, function(x) x$xy))
      rx<-range(xxx[,1])
      ry<-range(xxx[,2])
    }
    
    for (kk in names(x)) {
      if (!same4all) {
        rx<-range(x[[kk]]$xy[,1])
        ry<-range(x[[kk]]$xy[,2])
      }
      if(!add)
        plot(x[[kk]]$xy, ty="n", asp=1, main=kk,
             xlim=rx, ylim=ry,...)
      gr<-grey(x[[kk]]$area$levels/100)
      li2<-x[[kk]]$polygons
      for (i in length(li2):1)
        plot(li2[[i]], poly.args=list(col=gr[i], border=0), add=T)
      if (add.points)
        points(x[[kk]]$xy, pch=pch, bg=bg, col=col, cex=cex)
    }
  }


## example
library(adehabitat)
library(gpclib)
data(puechabon)
data(chamois)
xy2<-chamois$locs

aa<-NNCH(xy2)
aa
plot(aa)
plot(NNCH.area(aa))