Reference and Disclaimer: For this exercise, I have used explanation and codes from the book Spatial Ecology and Conservation Modeling by Robert Fletcher & Marie-Josée Fortin.

Load the necessary packages

rm(list = ls())

## First specify the packages of interest
packages = c("raster", "reshape2","rgdal","adehabitatLT", "adehabitatHR","adehabitatHS","survival")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)
par(mfrow = c(1, 1))

Import the required data and check for the projections

land <- raster("./HR_exercise/ccap_SE_FL_11class.tif")
#check projection
projection(land)
# Look inside gdb
fgdb <- "./HR_exercise/HR_exercise.gdb/"
ogrListLayers(fgdb)

#Add point data
point.data <- readOGR("./HR_exercise/HR_exercise.gdb","allpoints_XYTableToPoint")
projection(point.data)

#Re-project as land raster
crs.land <- projection(land)
point.data <- spTransform(point.data, crs.land)

Inspect the data

point.data@data <- na.omit(point.data@data)
point.data$IdNr <- as.factor(point.data$IdNr)
#str(point.data)

#convert to POSIXct object
point.data$Date <- as.POSIXct(point.data$GPSTime, tz="US/Central", format = "%Y/%m/%d %H:%M")
#summary(point.data)

head(point.data)
##   IdNr   GpsNumber                GPSTime                SMSTime Latitude
## 1    1 48797730423 2020/01/01 09:00:00+00 2020/01/01 13:02:59+00 26.48120
## 2    1 48797730423 2020/02/01 09:00:00+00 2020/02/01 10:01:00+00 26.48122
## 3    1 48797730423 2020/03/01 09:00:00+00 2020/03/01 10:31:00+00 26.48147
## 5    1 48797730423 2020/05/01 09:01:00+00 2020/05/01 11:32:00+00 27.54762
## 6    1 48797730423 2020/06/01 09:59:59+00 2020/06/01 12:32:00+00 27.54760
## 7    1 48797730423 2020/07/01 00:01:00+00 2020/07/01 02:00:59+00 27.54760
##   Longtitude BatteryVoltage GpsDescription Temperature GPSIntervals
## 1  -80.14165            381         WIUS41          23            b
## 2  -80.14165            388         WIUS41          25            b
## 3  -80.14167            389         WIUS41          20            b
## 5  -80.49375            390         WIUS41          24            b
## 6  -80.49370            384         WIUS41          29            b
## 7  -80.49403            388         WIUS41          27            b
##   VHFTelemetry Activity GSMSignal N_satellites Speed_knots Altitude Light Acc_X
## 1            0      183         9           10         0.0       13     0 -0.07
## 2            0      106        13           10         0.1       59     0  0.03
## 3            0      138        11           10         0.0        0     0 -0.01
## 5            0      658         9           11         0.2       77     0  0.02
## 6            0       29         7           10         0.0       19     0 -0.09
## 7            0       20         7            8         0.1        0   170  0.06
##   Acc_Y Acc_Z                Date
## 1  0.35 -0.91 2020-01-01 09:00:00
## 2  0.26 -0.95 2020-02-01 09:00:00
## 3  0.43 -0.90 2020-03-01 09:00:00
## 5  0.34 -0.92 2020-05-01 09:01:00
## 6  0.26 -0.91 2020-06-01 09:59:00
## 7  0.60 -0.77 2020-07-01 00:01:00

Plot

plot(land)
points(point.data, col=point.data$IdNr)

MCP Home range analysis

data_sub <- subset(point.data,GpsNumber == 48797732938)

mcp50 <- mcp(data_sub[,"GpsNumber"], percent = 50)
mcp90 <- mcp(data_sub[,"GpsNumber"], percent = 90)

head(mcp90@polygons)
## [[1]]
## An object of class "Polygons"
## Slot "Polygons":
## [[1]]
## An object of class "Polygon"
## Slot "labpt":
## [1] 1550620.9  570989.3
## 
## Slot "area":
## [1] 3949939840
## 
## Slot "hole":
## [1] FALSE
## 
## Slot "ringDir":
## [1] 1
## 
## Slot "coords":
##       coords.x1 coords.x2
## 2794    1567821  586139.0
## 1264    1583415  559501.6
## 1850    1590799  513386.9
## 2005    1590807  513269.4
## 4412    1589027  508817.6
## 4403    1588128  509033.8
## 2085    1503394  591215.8
## 2041    1503385  591246.2
## 2016    1503384  591273.9
## 201     1514960  608713.0
## 4335    1515131  608887.9
## 565     1515185  608928.3
## 723     1518096  610977.3
## 259     1527973  617694.7
## 489     1547386  608768.5
## 27941   1567821  586139.0
## 
## 
## 
## Slot "plotOrder":
## [1] 1
## 
## Slot "labpt":
## [1] 1550620.9  570989.3
## 
## Slot "ID":
## [1] "48797732938"
## 
## Slot "area":
## [1] 3949939840

Plot: MCP Homerange

plot(land, main="MCP Homerange")
plot(data_sub, add=TRUE, col=data_sub$IdNr)
plot(mcp90, add=TRUE)
plot(mcp50, add=TRUE, col="orange")

Local convex hull home range

data_sub$IdNr <- as.factor(data_sub$IdNr)

#subset
data_sub175 <- data_sub[data_sub$IdNr==175, ]
data_sub250 <- data_sub[data_sub$IdNr==250, ]

Plot:

plot(land, main="Local convex hull homerange")
plot(data_sub175, add=TRUE, col="blue")
plot(data_sub250, add=TRUE, col="red")

#initialize
k.int <- round(nrow(coordinates(data_sub175))^0.5,0)
a.int <- round(max(dist(coordinates(data_sub175))),0)

k.search <- seq(k.int, 10*k.int, by=50)
a.search <- seq(a.int, 2*a.int, by=3000)

#Parameter search for locoh-a
LoCoH.a.range <- LoCoH.a.area(SpatialPoints(coordinates(data_sub175)), unout="km2", arange=a.search)

#Parameter search for locoh-k
# LoCoH.k.range <- LoCoH.k.area(SpatialPoints(coordinates(data_sub175)), unout="km2", krange=k.search)

#plot
plot(LoCoH.a.range)

# plot(LoCoH.k.range)

#inspect
# a.search[21]
# k.search[11]

# #re-fit model
# LoCoH.k.61 <- LoCoH.k(SpatialPoints(coordinates(data_sub175)), k=k.search[11])

# plot(LoCoH.k.61)

#re-fit model
LoCoH.a.186259 <- LoCoH.a(SpatialPoints(coordinates(data_sub175)), a=a.search[21])
plot(LoCoH.a.186259)

Resource Selection

#Add IBIS data for first half of the year 2020
point.data <- readOGR("./HR_exercise/HR_exercise.gdb","WIBIS41_Jan1_June30_2020")
## OGR data source with driver: OpenFileGDB 
## Source: "D:\UGA\Academic Course Work\Fall-2022\FANR-8400-AdvGIS\Assignments\Ex4\HR_exercise\HR_exercise.gdb", layer: "WIBIS41_Jan1_June30_2020"
## with 2408 features
## It has 20 fields
# #check projection
# projection(point.data)
#label projection for later use
crs.land <- projection(land)
point.data <- spTransform(point.data, crs.land)
#convert to POSIXct object
point.data$Date <- as.POSIXct(point.data$GPSTime, tz="UTC", format = "%Y/%m/%d %H:%M")
#make trajectory object
point.ltraj <- as.ltraj(xy=coordinates(point.data), date=point.data$Date, id=point.data$GpsNumber, typeII=T)

Plot

par(mfrow = c(1, 3))
plot(point.ltraj)
#For kicks - what does the step length distribution look like?
plot(density(point.ltraj[[1]]$dist, na.rm=T), xlim=c(0,1000))
#turn angle distribution?
plot(density(abs(point.ltraj[[1]]$rel.angle), na.rm=T))

# FPT analysis

help(redisltraj)
## starting httpd help server ... done

Step 1: Create locations at evenly spaced intervals.

path.re <- redisltraj(point.ltraj,u=3600,nnew=10000)
plot(path.re)

Step 2: Analysis done for circles with radii values from 10m - 500m in 10m

coon.fpt <- fpt(path.re,radii=seq(from=3600,to=50000,by=1000))

Step 3: Variance plot to determine scale of Area-Restricted Search (ARS) behavior

varlogfpt(coon.fpt, graph=TRUE)

The scale value with the highest Variance is most appropriate scale of ARS behavior

Step 4: data frame of variance for each radii

coon.var <- varlogfpt(coon.fpt, graph = FALSE)

Step 5: plot FPT of a specific radii along the path

plot(coon.fpt, scale=8000)

Step 6: Append the FPT values for circles with r=40m to the dataframe of the locations

coon.r40 <- data.frame(path.re[[1]], fpt40 = coon.fpt[[1]]$r4)

Plot the path color-coded by ARS

plot(y~x, data=coon.r40, typ="l")
points(y~x, data=coon.r40, pch=19, col= ifelse(coon.r40$fpt >= 3700, "red", "green"),
       cex=ifelse(coon.r40$fpt >= 3700, 1, 0.5))
legend("topright", col="red", pch=19, legend="ARS")