Share

4.1 Kernel Density Estimation (KDE) with reference bandwidth selection (href)

In KDE, a kernel distribution (i.e. a three-dimensional hill or kernel) is placed on each telemetry location. The height of the hill is determined by the bandwidth of the distribution, and many distributions and methods are available (e.g. fixed versus adaptive, univariate versus bivariate bandwidth). We will focus here on "fixed kernel" but will alter the bandwidth selection. Datasets for avian and mammalian species can include as many as 10,000 locations and only the reference or default bandwidth (href) was able to produce KDE in both Home Range Tools and adehabitat or adehabitatHR (Calenge 2007, 2011). Estimation with (href) typically is not reliable for use on multimodal datasets because it results in over-smoothing of home ranges and multimodal distribution of locations is typical for most species (Worton 1995, Seaman et al. 1999).

  1. Exercise 4.1 - Download and extract zip folder into your preferred location
  2. Set working directory to the extracted folder in R under File - Change dir...
  3. First we need to load the packages needed for the exercise
    library(adehabitat) 
    NOTE: adehabitatHR package requires Spatial Points instead of data frame so code below will not work if adehabitatHR is also loaded
    library(sp)
    library(rgdal)
    library(raster)

  4. Now open the script "HrefScript.R" and run code directly from the script
    panther <- read.csv("pantherjitter.csv", header=T)
    str(panther)
  5. Now we can run fixed kernel home range with href bandwidth selection

    #Note below the code uses the original adehabitat package to run home range

    loc <- panther[, c("x", "y")]
    ## Estimation of UD for each animal separately
    id <- panther[, "id"]
    udbis <- kernelUD(loc, id, h = "href")
    ud <- kernelUD(loc, id = id, h = "href", grid = 40, same4all = FALSE,
    hlim = c(0.1, 1.5), kern = c("bivnorm"), extent = 0.5)
    image(ud) ## Note that the contours are under the locations

  6. Calculation of the 95 percent home range
    ver <- getverticeshr(ud, 95)
    plot(ver, add=TRUE)

  7. We estimated size of home range using href now we need to output the size of these home ranges that were estimated 

    #Look at the estimates of home range by contour
    cuicui1 <- kernel.area(loc, id)
    plot(cuicui1)
    cuicui1

    #Results for 2 animals by probability contour (i.e., 20-95 percent)
    FP121 FP128
    20 224.7252 346.8085
    25 304.0400 468.6602
    30 396.5739 609.2582
    35 508.9366 759.2295
    40 634.5183 946.6936
    45 786.5383 1162.2773
    50 971.6062 1424.7270
    55 1189.7218 1743.4159
    60 1440.8854 2165.2101
    65 1731.7062 2652.6167
    70 2062.1845 3205.6357
    75 2445.5394 3824.2671
    80 2928.0377 4546.0038
    85 3529.5082 5464.5778
    90 4316.0465 6748.7067
    95 5545.4257 8820.1848

  8. We can export these estimates with the previous code
    write.table(cuicui1,"output.csv", row.names=TRUE, sep=" ", col.names=TRUE, quote=TRUE, na = "NA")

  9. We can also output the respective shapefiles with code below:

    #######
    # OVERRIDE the default kver2spol function so that we can include the projection info
    #######

    kver2spol <- function(kv,projstr){
    x <- kv
    if (!inherits(x, "kver"))
    stop("x should be of class \"kver\"")
    if (!require(sp))
    stop("sp package needed")
    lipols <- lapply(1:length(x), function(i) {
    y <- x[[i]]
    class(y) <- c("data.frame", "list")
    res <- split(y[, 2:3], y[, 1])
    lipol <- lapply(res, function(z) {
    if (sum(abs(z[1, ] - z[nrow(z), ])) > 1e-16)
    z <- rbind(z, z[1, ])
    Polygon(as.matrix(z))
    })
    pols <- Polygons(lipol, ID = names(x)[i])
    return(pols)
    })
    return(SpatialPolygons(lipols, proj4string=CRS(as.character(projstr))))}

    #############
    # Function to export specific levels of isopleths of a "kv" object
    #############

    #Code creates contours only for animal 1 at each level so need to repeat for each animal if needed
    kv<-list()
    class(kv) <- "kver"
    kvtmp <- getverticeshr(udbis, lev = 99)
    kv$KHR99<- kvtmp[[1]]
    kvtmp <- getverticeshr(udbis, lev = 95)
    kv$KHR95<- kvtmp[[1]]
    kvtmp <- getverticeshr(udbis, lev = 90)
    kv$KHR90<- kvtmp[[1]]
    kvtmp <- getverticeshr(udbis, lev = 75)
    kv$KHR75<- kvtmp[[1]]
    kvtmp <- getverticeshr(udbis, lev = 50)
    kv$KHR50<- kvtmp[[1]]
    kvtmp <- getverticeshr(udbis, lev = 25)
    kv$KHR25<- kvtmp[[1]]

    spolTmp <- kver2spol(kv,"+proj=utm +zone=17N +ellps=WGS84")
    dfTmp <- data.frame(Isopleth=c("99","95","90","75","50","25"), row.names=c("KHR99","KHR95","KHR90","KHR75","KHR50","KHR25"))
    spdfTmp <- SpatialPolygonsDataFrame(spolTmp, dfTmp, match.ID = TRUE)
    writeOGR(spdfTmp,"HREF","FP048HREF", "ESRI Shapefile")
  10. We can also use package adehabitatHR to determine home range but requires different code


    panther <- subset(panther, panther$CatID == "143")
    panther$CatID <- factor(panther$CatID)
    loc <- data.frame("x"=panther$x,"y"=panther$y)
    proj4string <- CRS("+proj=utm +zone=17N +ellps=WGS84")
    cats <- SpatialPointsDataFrame(loc,panther, proj4string = proj4string)
    udbis <- kernelUD(cats[,1], h = "href")
    image(udbis)
    ver <- getverticeshr(udbis, standardize = FALSE)
    ver50 <- getverticeshr(udbis, percent=50)
    ver80 <- getverticeshr(udbis, percent=80)
    ver90 <- getverticeshr(udbis, percent=90)
    ver95 <- getverticeshr(udbis, percent=95)
    ver99 <- getverticeshr(udbis, percent=99)
    ver
    plot(ver99, col="grey",axes=T);plot(ver95, add=T);plot(ver90, add=T);
    plot(ver80, add=T);plot(ver50, add=T)
    points(cats)