This code is rough and ugly, but it worked for my purpose. Since writing this in 2015, I've updated the actual script I use to be much more efficient and intuitive, but I haven't updated this blog post to reflect that. Contact me or comment with any questions.
install.packages("adehabitatHR") #install if needed
library(adehabitatHR) #load package for home range (HR) analyses
all.telemetry <- read.csv(file.choose()) #load all telemetry data
head(all.telemetry) #view data
all.telemetry$SnakeID <- as.factor(all.telemetry$SnakeID) #turn first column into categorical data if needed
s1<- subset(all.telemetry,SnakeID==1) # only look at Snake #1's data, OR: only look at Frank Underwood's data
# s1<- subset(all.telemetry,Name=="Frank Underwood")
head(s1)
s1xy<- SpatialPoints(s1[2:3]) #turn the xy data into coordinates
head(s1xy)
s1_homerange <- mcp(s1xy,percent=100)$area #find the 100% minnimum convex polygon value
s1_homerange
plot(mcp(s1xy,percent=100)) #graph the 100% mcp
plot(s1xy,cex=.75,pch=1,add=TRUE) #add the data points on top
ud <- kernelUD(s1xy,h="href") #calculate the utilization distribution
ver<-getverticeshr(ud,percent=95)$area #find the 95% utilization distribution
ver #print the 95% ud
kern <- kernel.area(ud, percent = seq(20, 95, by = 5)) #find the ud for 20 - 95%
kern #print the uds
plot(getverticeshr(ud,percent=95)) #view the 95% kernel to the plot
plot(getverticeshr(ud,percent=50),col="red",add=TRUE) #view the 50% kernel to the plot
plot(mcp(s1xy,percent=100),add=T) #add the 100% mcp
plot(s1xy,cex=1,pch=21,col="white",add=TRUE)
h <- unlist(ud@h[1],use.names=F) #find the smoothing factor from the utilization distribution
h #view the smoothing parameter
output <- matrix(ncol=2, nrow=h) #create an empty matrix to fill with home range and smoothing factor values
output <- data.frame(output) #turn it into a dataframe
colnames(output) <- c("h","hr") #add column names
for (i in 1:h){ #start the loop to find smoothing factor values
ud_loop <- tryCatch(kernelUD(s1xy,h=i),error=function(err) NA) #trycatch allows loop to progress despite NAs
output[i,1]=i #fill first column in dataframe
output[i,2]=tryCatch(getverticeshr(ud_loop)$area,error=function(err) NA) #fill second column in dataframe
} #end loop
output #view filled maxtrix
plot(output,type="l",lwd=2) #graph of home range against smoothing factor
abline(s1_homerange,0) #add a line at the 100% mcp homerange value
hvalue <- locator(n=1)$x #physically click on the intersection between the two lines on the graph, then finds the smoothing factor value associated with the 100% mcp
hvalue #print the smoothing factor
ud2 <- kernelUD(s1xy,h=hvalue) #calculate the new utilization distribution
ver2<-getverticeshr(ud2,percent=95)$area #find the new 95% utilization distirbution
s1_homerange
ver2 #print the new 95% ud
kern2 <- kernel.area(ud2, percent = seq(20, 95, by = 5)) #find the new uds for 20 - 95%
kern2 #print the uds
plot(getverticeshr(ud2,percent=95)) #view the new 95% kernel
plot(getverticeshr(ud2,percent=50),col="red",add=TRUE) #add the new 50% kernel to the plot
plot(mcp(s1xy,percent=100),add=T) #add the 100% mcp
plot(s1xy,cex=1,pch=21,col="white",add=TRUE)
getverticeshr(ud,percent=50)$area # old 50% ud
getverticeshr(ud2,percent=50)$area # new 50% ud
getverticeshr(ud,percent=50)$area - getverticeshr(ud2,percent=50)$area # difference
install.packages("adehabitatHR") #install if needed
library(adehabitatHR) #load package for home range (HR) analyses
all.telemetry <- read.csv(file.choose()) #load all telemetry data
head(all.telemetry) #view data
all.telemetry$SnakeID <- as.factor(all.telemetry$SnakeID) #turn first column into categorical data if needed
s1<- subset(all.telemetry,SnakeID==1) # only look at Snake #1's data, OR: only look at Frank Underwood's data
# s1<- subset(all.telemetry,Name=="Frank Underwood")
head(s1)
s1xy<- SpatialPoints(s1[2:3]) #turn the xy data into coordinates
head(s1xy)
s1_homerange <- mcp(s1xy,percent=100)$area #find the 100% minnimum convex polygon value
s1_homerange
plot(mcp(s1xy,percent=100)) #graph the 100% mcp
plot(s1xy,cex=.75,pch=1,add=TRUE) #add the data points on top
ud <- kernelUD(s1xy,h="href") #calculate the utilization distribution
ver<-getverticeshr(ud,percent=95)$area #find the 95% utilization distribution
ver #print the 95% ud
kern <- kernel.area(ud, percent = seq(20, 95, by = 5)) #find the ud for 20 - 95%
kern #print the uds
plot(getverticeshr(ud,percent=95)) #view the 95% kernel to the plot
plot(getverticeshr(ud,percent=50),col="red",add=TRUE) #view the 50% kernel to the plot
plot(mcp(s1xy,percent=100),add=T) #add the 100% mcp
plot(s1xy,cex=1,pch=21,col="white",add=TRUE)
h <- unlist(ud@h[1],use.names=F) #find the smoothing factor from the utilization distribution
h #view the smoothing parameter
output <- matrix(ncol=2, nrow=h) #create an empty matrix to fill with home range and smoothing factor values
output <- data.frame(output) #turn it into a dataframe
colnames(output) <- c("h","hr") #add column names
for (i in 1:h){ #start the loop to find smoothing factor values
ud_loop <- tryCatch(kernelUD(s1xy,h=i),error=function(err) NA) #trycatch allows loop to progress despite NAs
output[i,1]=i #fill first column in dataframe
output[i,2]=tryCatch(getverticeshr(ud_loop)$area,error=function(err) NA) #fill second column in dataframe
} #end loop
output #view filled maxtrix
plot(output,type="l",lwd=2) #graph of home range against smoothing factor
abline(s1_homerange,0) #add a line at the 100% mcp homerange value
hvalue <- locator(n=1)$x #physically click on the intersection between the two lines on the graph, then finds the smoothing factor value associated with the 100% mcp
hvalue #print the smoothing factor
ud2 <- kernelUD(s1xy,h=hvalue) #calculate the new utilization distribution
ver2<-getverticeshr(ud2,percent=95)$area #find the new 95% utilization distirbution
s1_homerange
ver2 #print the new 95% ud
kern2 <- kernel.area(ud2, percent = seq(20, 95, by = 5)) #find the new uds for 20 - 95%
kern2 #print the uds
plot(getverticeshr(ud2,percent=95)) #view the new 95% kernel
plot(getverticeshr(ud2,percent=50),col="red",add=TRUE) #add the new 50% kernel to the plot
plot(mcp(s1xy,percent=100),add=T) #add the 100% mcp
plot(s1xy,cex=1,pch=21,col="white",add=TRUE)
getverticeshr(ud,percent=50)$area # old 50% ud
getverticeshr(ud2,percent=50)$area # new 50% ud
getverticeshr(ud,percent=50)$area - getverticeshr(ud2,percent=50)$area # difference
https://johnpvanek.wordpress.com/2016/05/02/find-the-appropriate-smoothing-factor-h-for-kernel-densities-in-r