Set up

Load packages

#clear environment
rm(list = ls()) 

# install.packages("tlocoh", dependencies=T, repos=c("http://R-Forge.R-project.org" ,
# "http://cran.cnr.berkeley.edu"))
library(tlocoh)
library(sp)
library(rgdal)
library(rgeos)
library(sf)
library(ggmap)
library(dplyr)

#get some data
data(toni) #attach data

The data is from the GPS collar of a single buffalo (“Toni”) in South Africa

sf.toni <- st_as_sf(toni, coords = c("long","lat"))
head(sf.toni)
## Simple feature collection with 6 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 31.73874 ymin: -24.1695 xmax: 31.75345 ymax: -24.15329
## CRS:           NA
##         id           timestamp.utc                   geometry
## 17930 toni 2005-08-23 06:35:00.000  POINT (31.75345 -24.1695)
## 17931 toni 2005-08-23 07:34:00.000 POINT (31.73884 -24.15402)
## 17932 toni 2005-08-23 08:34:00.000 POINT (31.73969 -24.15359)
## 17933 toni 2005-08-23 09:35:00.000 POINT (31.73874 -24.15329)
## 17934 toni 2005-08-23 10:34:00.000 POINT (31.73946 -24.15336)
## 17935 toni 2005-08-23 11:35:00.000 POINT (31.73898 -24.15363)

The essential parts of GPS tracking data is the GPS coordinates (long, lat etc), a time stamp (date, time), and ID of the animal (especially important if there are multiple individuals being tracked)

Visualize:

sf.toni$date <-as.Date(as.POSIXct(sf.toni$timestamp.utc, tz="UTC")) #convert times to POSIX

bbox = c(left=min(toni$long - 0.05),right=max(toni$long) + 0.05, #longitude
             top=max(toni$lat) + 0.05,bottom=min(toni$lat) -0.05) #latitude


#Open street maps
basemap <- ggmap(get_stamenmap(maptype='terrain', 
                            bbox=bbox,
                            zoom=12)) 
basemap+
  theme_void()+
  geom_sf(data=sf.toni, aes(color=date), inherit.aes = FALSE)

Kernel density estimation

Estimates density of spatial points within a given area; a smooth curve is fit over each point and then the height of these curves are summed over the space.

refs:

library(ks)

coords <- st_coordinates(sf.toni)
H.scv <- Hscv(x=coords) #Smoothed cross-validation (SCV) bandwidth selector
NumGridLines =10

min.x <- bbox[["left"]]
min.y <- bbox[["bottom"]]
max.x <- bbox[["right"]]
max.y <- bbox[["top"]]

fhat <- ks::kde(x=coords, H=H.scv, gridsize=NumGridLines, binned=TRUE,
                    xmin=c(min.x, min.y), 
                    xmax=c(max.x, max.y)) #kernel density estimater

#Function to extract contour polygons
Fx_ExtractContour <-  function(fhat, ContourLine){
  #Get a few other contours 
  contour.x <- with(fhat,contourLines(x=eval.points[[1]],y=eval.points[[2]],
                                      z=estimate,levels=cont[ContourLine])[[1]]) 
  poly.x <- with(contour.x, data.frame(x,y))
  poly.x <- rbind(poly.x,poly.x[1,])    # close polygon
  poly.x <- as.data.frame(poly.x) %>%
    mutate(contour = ContourLine) %>%
    sfheaders::sf_polygon(x="x", y="y", polygon_id="contour")
  
  return(poly.x)
}

contours <- rbind(
  Fx_ExtractContour(fhat, "99%"),
  Fx_ExtractContour(fhat, "95%"),
  Fx_ExtractContour(fhat, "75%"),
  Fx_ExtractContour(fhat, "50%"),
  Fx_ExtractContour(fhat, '10%')
  )

Pal <- c("99%"="red4","95%"="red2", "75%"="orange3", "50%"="goldenrod3", "10%"="lightgoldenrod3")
basemap + 
  geom_sf(data=contours, aes(fill=contour), inherit.aes = FALSE, alpha = 0.6, color=NA)+
  scale_fill_manual(values = Pal, name= "Contour density") +
    geom_sf(data=sf.toni, inherit.aes = FALSE, size =0.01)

t-LoCoH

This is the code provided in the t-locoh R package starting tutorial. For more detail, visit the t-locoh pacakge website.

Refs:

Note: the package we are using is based on the ‘sp’ r spatial pacakge

Format data for t-locoh

#convert to sp object, set projection and timezones
toni.sp.latlong <- SpatialPoints(toni[ , c("long","lat")], proj4string=CRS("+proj=longlat +ellps=WGS84")) #sp object
toni.sp.utm <- spTransform(toni.sp.latlong, CRS("+proj=utm +south +zone=36
+ellps=WGS84")) #convert to utm
toni.mat.utm <- coordinates(toni.sp.utm) #extract coords
colnames(toni.mat.utm) <- c("x","y")
toni.gmt <- as.POSIXct(toni$timestamp.utc, tz="UTC") #convert times to POSIX
local.tz <- "Africa/Johannesburg" #specify timezon
toni.localtime <- as.POSIXct(format(toni.gmt, tz=local.tz), tz=local.tz)

#convert to locoh object
toni.lxy <- xyt.lxy(xy=toni.mat.utm, dt=toni.localtime, id="toni",
proj4string=CRS("+proj=utm +south +zone=36 +ellps=WGS84"))

Create hullsets:

toni.lxy <- lxy.nn.add(toni.lxy, s=0.003, k=25) #find nearest neighbors
toni.lhs <- lxy.lhs(toni.lxy, k=15, s=0.003) #create hulls

Create isopleths (aggregated hulls):

toni.lhs <- lhs.iso.add(toni.lhs) #add iso layers
toni.lhs.k15 <- lhs.select(toni.lhs, k=15) #select the k=15 for analysis
toni.lhs.k15 <- lhs.ellipses.add(toni.lhs.k15) #add ellipses
toni.lhs.k15 <- lhs.iso.add(toni.lhs.k15, sort.metric="ecc") #calculate eccentricity
toni.lhs.k15 <- lhs.visit.add(toni.lhs.k15, ivg=3600*12) #cycle period

Visualize how the algorthm works

plot(toni.lhs.k15, hulls=T, ellipses=T, allpts=T, nn=T, ptid="auto")

Density isopleths (aggregated hulls by density of points)

plot(toni.lhs, iso=T, k=15, allpts=T, cex.allpts=0.1, col.allpts="gray30",
ufipt=F); plot(toni.lhs, iso=T, k=15, ufipt=F)

Eccentricity isopleths (aggregated hulls by eccentricity of ellipses)

plot(toni.lhs.k15, iso=T, iso.sort.metric="ecc")

Revistiation rate (frequency of visit to sites)

#hist(toni.lhs.k15, metric="nsv")
plot(toni.lhs.k15, hpp=T, hpp.classify="nsv", ivg=3600*12, col.ramp="rainbow")

Duration rate (length of time at each site)

 plot(toni.lhs.k15, hpp=T, hpp.classify="mnlv", col.ramp="rainbow")

Revistiation vs duration

hsp <- lhs.plot.scatter(toni.lhs.k15, x="nsv", y="mnlv", col="spiral",
bg="black"); plot(toni.lhs.k15, hpp=T, hsp=hsp, hpp.classify="hsp")