R/clusthr.r

In ClementCalenge/adehabitatHR: Home Range Estimation

### XX plot de la surface en fonction du level

spoldf2MCHu <- function(spdf, nam="a")
{
    if (!inherits(spdf,"SpatialPolygonsDataFrame"))
        stop("spdf should be of class SpatialPolygonsDataFrame")
    if ((ncol(spdf@data)!=2)||(!all(c("area","percent")%in%names(spdf@data))))
        stop("The @data component should contain two columns named percent and area")
    res <- list(spdf)
    names(res) <- "a"
    class(res) <- "MCHu"
    return(res)
}



MCHu2hrsize <- function(x, percent=seq(20,100, by=10), plotit=TRUE)
{
    if (!inherits(x,"MCHu")&!inherits(x, "SpatialPolygonsDataFrame"))
        stop("x should have been generated by LoCoH.* or clusthr")

    if (inherits(x, "SpatialPolygonsDataFrame")) {
        df <- as.data.frame(x)
        percent <- sort(percent)
        ind <- 1:nrow(df)
        lb <- unlist(lapply(percent, function(p) {
            min(ind[df$percent>=p])
        }))
        df <- df$area[lb]
        names(df) <- percent
        df <- data.frame(df)
        class(df)<-c("hrsize", "data.frame")
        return(df)
    } else {
        ar <- do.call("data.frame", lapply(x, MCHu2hrsize, percent=percent))
        names(ar) <- names(x)
        row.names(ar)<-percent
        class(ar)<-c("hrsize", "data.frame")
        if (plotit)
            plot(ar)
        return(ar)
    }
}

clusthr <-function(xy, unin = c("m", "km"),
                   unout = c("ha", "m2", "km2"),
                   duplicates=c("random","remove"), amount = NULL)
{
    ## Verifications
    if (!inherits(xy, "SpatialPoints"))
        stop("xy should inherit the class SpatialPoints")
    pfs <- proj4string(xy)
    if (inherits(xy, "SpatialPointsDataFrame")) {
        if (ncol(xy)!=1) {
            warning("xy contains more than one column. id ignored")
            m <- 1
            id<-rep("a", nrow(coordinates(xy)))
        } else {
            id <- xy[[1]]
            m <- 2
        }
    } else {
        m <- 1
        id<-rep("a", nrow(coordinates(xy)))
    }
    id<-factor(id)
    duplicates <- match.arg(duplicates)
    unin <- match.arg(unin)
    unout <- match.arg(unout)


    xy <- as.data.frame(coordinates(xy))


    ## splits the coordinates into one component per animal
    lixy<-split(xy, id)

    ## the output list
    res<-list()

    ## The function clubase is used to compute the length
    ## of the output vectors. It relies on the C function "longfacclust".
    ## It is needed to reserve memory for the main C function "clusterhrr"
    ## called below
    clubase <- function(xy)
    {
        nr <- as.integer(nrow(xy))
        xy <- as.double(t(as.matrix(xy)))
        len2 <- integer(1)
        toto <- .C("longfacclustr", as.double(xy),
                   as.integer(nr), integer(1),
                   PACKAGE = "adehabitatHR")[[3]]
        return(toto)
    }

    ## for each animal, an home range is desired
    for (i in 1:length(lixy)){

        x<-as.matrix(lixy[[i]])

        ## Management of duplicates
        levv <- factor(apply(x,1,paste, collapse=" "))
        if (duplicates=="remove") {
            x <- as.data.frame(x)
            x <- as.matrix(do.call("rbind", lapply(split(x, levv), function(z) z[1,])))
        } else {
            x <- as.data.frame(x)
            sam1 <- x[,1] - jitter(x[,1], 1, amount)
            sam2 <- x[,1] - jitter(x[,1], 1, amount)
            lixyt <- split(x, levv)
            lisam1 <- split(sam1, levv)
            lisam2 <- split(sam2, levv)
            x <- as.matrix(do.call("rbind",
                                   lapply(1:length(lixyt),
                                          function(a){
                                              z <- lixyt[[a]]
                                              if (nrow(z)>1) {
                                                   z[,1] <- z[,1]+lisam1[[a]]
                                                   z[,2] <- z[,2]+lisam2[[a]]
                                               }
                                              return(z)
                                          })))
        }




        ## computation of the output vectors
        len <- clubase(x)

        ## Computation of the tree: call to the C function "clusterhrr"
        toto <- .C("clusterhrr",
                   as.double(t(as.matrix(x))), as.integer(nrow(x)),
                   integer(len), integer(len), integer(len),
                   as.integer(len), PACKAGE = "adehabitatHR")

        step <- toto[[3]]              # contains indices of the step of the
                                        # algorithm (First step =1, second = 2)

        reloc <- toto[[4]]            # contains the indices of the
                                        # relocations clustered at the step
                                        # indicated by the factor step

        clust <- toto[[5]]              # contains the indices of the cluster
                                        # in which the relocations clustered
                                        # at the step i are clustered


        ## The relocations clustered at step one
        liclu <- list()
        liclu[[1]] <- reloc[step==1]

        ## The home range at step one
        nlocc <- 100*length(reloc[step==1])/nrow(x)
        poltot <- list()
        tmp <- as.matrix(x[reloc[step==1],])
        tmp <- rbind(tmp, tmp[1,])
        pc <- Polygons(list(Polygon(tmp)), 1)
        arre <- .arcpspdf(SpatialPolygons(list(pc)))

        ## poltot contains the home range
        poltot[1] <- list(pc)


        ## For each step:
        for (j in 2:max(step)) {

            ## The relocations clustered at step j
            relocj <- reloc[step==j]
            r1 <- relocj[1]

            ## are the relocations already clustered (in which case the step is
            ## a merging of two clusters)
            oussa <- unlist(lapply(1:length(liclu), function(o)
                                   r1%in%liclu[[o]]))

            ## If it is a merging of two clusters
            if (any(oussa)) {
                ## we merge the two sets of relocations
                liclu[[which(oussa)]] <- liclu[[which(oussa)]][-c(which(liclu[[which(oussa)]]%in%relocj))]
            }

            ## update liclu according to the new cluster
            liclu[clust[step==j][1]] <- list(c(unlist(liclu[clust[step==j][1]]), relocj))

            ## computes the convex polygons around all current clusters
            kkk <- lapply(1:length(liclu), function(n) {
                k <- liclu[[n]]
                if (length(k)>0) {
                    xy2 <- x[k,]
                    tmp <- as.matrix(xy2[chull(xy2[,1], xy2[,2]),])
                    tmp <- rbind(tmp, tmp[1,])
                    pol <- Polygon(tmp)
                    are <- .arcpspdf(SpatialPolygons(list(Polygons(list(pol), 1))))
                    return(list(pol=pol, are=are))
                }
            })
            arre[j] <- sum(unlist(lapply(kkk, function(x) x$are)))
            kkk <- lapply(kkk, function(x) x$pol)
            kkk <- kkk[!unlist(lapply(kkk, is.null))]
            kkk <- Polygons(kkk, j)


            ## Computes the number of relocations already clustered
            nlocc[j] <- 100*length(unique(unlist(liclu)))/nrow(x)

            ## the home range is stored in poltot
            poltot[j] <- list(kkk)
        }

        if (unin == "m") {
            if (unout == "ha")
                arre <- arre/10000
            if (unout == "km2")
                arre <- arre/1e+06
        }
        if (unin == "km") {
            if (unout == "ha")
                arre <- arre * 100
            if (unout == "m2")
                arre <- arre * 1e+06
        }


        res[[i]] <- SpatialPolygonsDataFrame(SpatialPolygons(poltot), data.frame(percent=nlocc, area=arre))
        if (!is.na(pfs))
            proj4string(res[[i]]) <- CRS(pfs)
    }

    if (m==1) {
        return(res[[1]])
    } else {
        names(res) <- names(lixy)
        class(res) <- "MCHu"
        return(res)
    }
}