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).
- Exercise 4.1 - Download and extract zip folder into your preferred location
- Set working directory to the extracted folder in R under File - Change dir...
- 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) - Now open the script "HrefScript.R" and run code directly from the script
panther <- read.csv("pantherjitter.csv", header=T)
str(panther) - 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 - Calculation of the 95 percent home range
ver <- getverticeshr(ud, 95)
plot(ver, add=TRUE) - 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 - 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") - 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") - 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)
