Home Range Analysis
Advanced Data Analytics

1 Introduction

This exercise will build on the skills acquired in the Mapping Basics and Species Distribution Modeling exercises to perform Home Range Analysis and map the results. We will start by examining the structure of data used in home range analyses, we will look at ways to perform QA/QC checks on the location data, and finally various ways to visualize the data. One of the new skills we will discuss in this exercise is to create functions and loops to replicate analyses.

1.1 Packages used in this exercise

There are several specialty packages that will be used in this exercise due to the specific nature of the analyses. Some of these packages you will need to install while several others we have used in previous exercises and should already be installed.

adehabitatHR | ctmm | pbapply | raster | sf | sp |
| terra | terrainr | tidyverse |

To begin, we will install the following:

pacman::p_load("adehabitatHR", "ctmm", "pbapply", "raster", "sf", "sp", "terra", "terrainr", "tidyverse")

2 Dataset

While we have seen that biological data can be obtained from sites such as Dryad and GBIF, we need to obtain data with specific variables that might not be available on those sites. So to start this exercise, we will download a *.csv file containing information for Ophiophagus hannah (King Cobra) from movebank.org that contains the appropriate variable to track movement from radio transmitter data. Because this data has to be downloaded we will need to go through the process of importing this information using the read.csv() command. Once we have imported the data we can view the structure.

If forked, this data will be available in the Data folder within the repository. Additionally, ESRI Shapefile and Google KMZ formats will be included if you are interested in working with those data types.

Import Dataset with read.csv

data <- read.csv("./Data/ophiophagus_hannah.csv")

head(data, 3)

Notice that this dataset contains more than species information and x,y coordinates. While that information is important, this dataset also contains time and date information (timestamp), sensor type, and UTM data. In order to be sure that the dataset contains no outliers, we can plot the data using ggplot and plotly to interactively view the data.

qaqc_plot <- ggplot() + geom_point(data=data, 
                                   aes(utm.easting,utm.northing,
                                       color=individual.local.identifier)) +
                        labs(x="Easting", y="Northing") +
                        guides(color=guide_legend("Identifier"))

plotly::ggplotly(qaqc_plot)

With plotly, similar to leaflet, we have the ability to examine the spread of the data points and additional information from various columns. From this plot we can see that there are two individuals, OPHA1 and OPHA2, that were tracked for this study and have very little overlap in their apparent range.

While we could continue with the current dataset to calculate home range for the entire population, to simplify this exercises we will instead focus on individuals.

MCP Analysis Example	KDE Analysis Example

Depending on the needs of your analysis, you may find a situation where a dataset being provided by a collaborator or through a repository contains a large number of individuals being recorded. In order to separate these larger datasets into individuals, we can create a function that uses the lapply command to apply a function over a list or vector dataset. Specifically, this function will take the original dataset, split it into separate files based on the individual identifier, and create new *.csv files using the identifier as the filename.

The anatomy of the function is as follows:

• lapply(), apply the function over a list
• split(), separates the data
• function(), compose a series of steps to be applied to the data
• write.csv(), write a csv file
• paste(), create a character string to be used for the file name

If you examine the root folder you will now find the addition of two new *.csv files: OPHA1 and OPHA2. We will use these files to run our home range analyses. Alternatively, if you are provided a number of individual *.csv files and need to quickly import the data into separate data frames you could write a for loop that would look something like this:

list <- gsub("\\.csv$","", list.files(pattern="\\.csv$"))

for(i in list){
  assign(i, read.csv(paste(i, ".csv", sep="")))
}

All of the *.csv files in your root directory will now be available in your global environment. However, the method we used above created individual *.csv files for each individuals, so instead we will make a list of those individual files as reference for the analyses when necessary.

files <- list.files(path = ".", pattern = "[OPHA]+[0-9]+", full.names = TRUE)

In the list.files command above, path = "." informs the location, in this case the root directory, pattern = describes the way the files are named, in this case OPHA followed by a number between 0-9, and full.names describes how the files will be listed.

3 Analysis

To date, this will be one of the more code intensive exercises we have attempted. While the code itself is not difficult, the implementation is tedious due to the varying data types and analyses performed. Additionally, while not run as a component of this exericse, I am providing a script that will walk you through the steps to automate an analysis for a number of individuals using a function and pbapply which produces a progress bar in the console detailing the level of completion.

3.1 Imagery

In previous exercises we used a number of different sources to obtain imagery or basemaps for the various visualizations. As we have discussed with packages like openstreetmap or ggmap there can be restrictions that complicate their usefulness. Additionally, some sources of raster imagery can be limited in their geographic scope. In the previous exercise we used Google Earth Engine to produce images that could be used as basemaps for our analyses. So for this exercise we will use GEE to obtain a Sentinel 2 image from Thailand and utilize it as the background for our homerange analyses.

Because this code is javascript I will simply use a code block to provide details on the script I used to locate, prepare, and download the image. For more information refer to the previous exercise or GEE documentation.

// Set variables for the location and imagery source
var thailand = /* color: #d63000 */ee.Geometry.Point([101.96660539060161, 14.546641122346026]),
sentinel2 = ee.ImageCollection("COPERNICUS/S2_HARMONIZED");
 
// Dates below used are closest available to data collection date
var image = ee.ImageCollection(sentinel2) // identify the sensor
  .filterDate("2016-04-01", "2016-10-30") // filter specific dates
  .filterBounds(thailand) // only include scenes that intersect with our point geometry
  .sort("CLOUD_COVERAGE_ASSESSMENT") // filter by cloudy pixel percentage
  .first(); // select the least cloudy scene 
 
// Set visualization parameters, these can be adjusted within the code editor
var naturalcolor = {
  bands: ["B4", "B3", "B2"],
  min: 0,
  max: 4000,
};
 
// Display images on the map
Map.addLayer(image, naturalcolor, "Natural Color Composite");

// Prepare the image and export in UTM 47; remember the data is from Thailand
var exportRGB = image.visualize({
  bands: ["B4", "B3", "B2"],
  min: 0,
  max: 4000,
});

// CRS was set to appropriate EPSG and maxPixels were set to allow for large file size 
Export.image.toDrive({
  image: exportRGB,
  scale: 10,
  description: "thailand32647",
  crs: "EPSG:32647",
  maxPixels: 1e13,
});

Either set the save location to the project folder, or mannually copy the imagery over to make it available in R. There is a rgee package you can experiment with to incorporate this step into workflow within R however, there is a little learning curve due to the majority of the script being python-based. In order to make the resulting *.geoTIFF file suitable for use in ggplot, we’ll need to create an image stack and convert it to a dataframe. While this step isn’t completely necessary, it does make the syntax of using the image in ggplot more straightforwrd. Additionally, because image swath from Sentinel 2 is 290km, we will also create two sizes of crops for the varying spatial constraints of the data. This can be adjusted throughout the analysis as needed to fit the specific visualizations.

ortho <- "./Data/thailand32647.tif"

image_stack <- raster::stack(ortho)

crop_extent <- raster::extent(min(data$utm.easting-2000), 
                              max(data$utm.easting+2000), 
                              min(data$utm.northing-2000), 
                              max(data$utm.northing+2000))

crop_full <- raster::extent(min(data$utm.easting-5000), 
                            max(data$utm.easting+5000), 
                            min(data$utm.northing-5000), 
                            max(data$utm.northing+5000))

image_crop <- crop(image_stack, crop_extent)
image_full <- crop(image_stack, crop_full)

thailand <- as.data.frame(image_crop, xy = TRUE) %>% setNames(c("x", "y", "red", "green", "blue"))
thailand.full <- as.data.frame(image_full, xy = TRUE) %>% setNames(c("x", "y", "red", "green", "blue"))

ggplot()+geom_spatial_rgb(data = thailand, aes(x=x, y=y, r=red, g=green, b=blue))+coord_sf()

3.2 Home Range Analysis

There are three basic types of home range analyses we will perform in this exercise: Minimum Convex Polygon (MCP), Kernel-Density Estimation (KDE), and Brownian Bridge Movement Model (BB). There are a number of different tuning parameters that can be applied these analyses, however in this exercise we will use the most basic versions of the analysis.

In the section above we use the lapply command to loop a function used to separate the original dataset into individual files. This is a useful tool, however, when the function loops through dozens or even hundreds of files, the process can take a long period of time to complete. Using the pblapply command adds a progress bar (i.e. pb) to the process which provides an estimated time for completion of the function. We will use a similar process to the one above, using the pblapply command to run MCP analysis on the individual *.csv files.

3.2.1 Minimum Convex Polygon

One of the most simple homerange analyses is a minimum convex polygon or mcp. A mcp draws a polygon which represents a specified minimum bounding geometry enclosing the data points. This essentially connects the points along the farthest geographic points enclosing them all in a single shape. For this analysis we will focus on the OPHA1 file we separated from the larger dataset using the adehabitatHR package.

A description of the steps within the following code will be discussed following the output. The process works by establishing the function and then running pblapply referring back to the files we created in the dataset portion of this exercise.

## Data for a single individual
opha1 <- read.csv('OPHA1.csv')
xy <- c(as.data.frame(opha1$utm.easting),as.data.frame(opha1$utm.northing))

## projecting the data to UTM
data.proj <- SpatialPointsDataFrame(xy,opha1, proj4string = CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs"))
xy.sp <- SpatialPoints(data.proj@coords)

## MCP analysis
mcp.out <- mcp(xy.sp, percent=100, unout="ha")

## converting mcp to sf object for ggplot and points to df
mcp.sf <- st_as_sf(mcp.out)

mcp.points <- cbind((data.frame(xy)),opha1$individual.local.identifier)
colnames(mcp.points) <- c("x","y", "identifier")

# MCP plot
mcp.plot <- ggplot() + theme_bw() +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5, vjust = 0.5)) + 
  theme(panel.grid = element_blank()) +
  geom_spatial_rgb(data = thailand, aes(x=x, y=y, r=red, g=green, b=blue)) +
  geom_sf(data=mcp.sf, alpha = 0.1) +
  geom_point(data=mcp.points, aes(x=x, y=y), color = "lightgreen") + 
  labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
  coord_sf(xlim = c(820000, 814250),  ylim = c(1604500, 1608000))
mcp.plot + annotate("text", x = 815000, y = 1607500, label = paste(round(mcp.out@data$area,2),"ha"), 
                    fontface = "bold", size = 5, color = "white")

The anatomy of the script above is as follows:

• read.csv(), read the data from a given *.csv file
• xy, create a coordinate file of the easting and northing data
• SpatialPointsDataFrame(), project the data to UTM Zone 47
• mcp(), computes home range using the Minimum Convex Polygon estimator
• st_as_sf, converts to a data frame for plotting with ggplot2
• colnames(), renames the columns
• annotate, adds text to the plot

The rest of the information can be found by examining the help menu for the various commands or looking at the display options for each portion of the script.

3.2.2 Kernel-Density Estimation

The same design above can be used to create KDE plots for each individual in the dataset. Because the process is essentially the same, only portions of the script that were not previously described will be discussed.

## Data for a single individual
opha2 <- read.csv('OPHA2.csv')
x <- as.data.frame(opha2$utm.easting)
y <- as.data.frame(opha2$utm.northing)
xy <- c(x,y)

## projecting the data to UTM
data.proj <- SpatialPointsDataFrame(xy,opha2, proj4string = CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs"))
xy.sp <- SpatialPoints(data.proj@coords)

## KDE analysis at 95, 85, 75
kde<-kernelUD(xy.sp, h="href", kern="bivnorm", grid=100)
ver95 <- getverticeshr(kde, 95)
ver85 <- getverticeshr(kde, 85)
ver75 <- getverticeshr(kde, 75)

## converting to sf object and df
kde.points <- cbind((data.frame(data.proj@coords)),opha2$individual.local.identifier)
colnames(kde.points) <- c("x","y","identifier")
kde95.sf <- st_as_sf(ver95)
kde85.sf <- st_as_sf(ver85)
kde75.sf <- st_as_sf(ver75)

##making the plot
kde.plot <- ggplot() + theme_bw() +
  theme(axis.text.y = element_text(angle = 90, hjust = 0.5, vjust = 0.5)) + 
  geom_spatial_rgb(data = thailand, aes(x=x, y=y, r=red, g=green, b=blue)) +
  geom_sf(data=kde95.sf, fill = "blue", alpha = 0.4) +
  geom_sf(data=kde85.sf, fill = "green", alpha = 0.3) +
  geom_sf(data=kde75.sf, fill = "yellow", alpha = 0.2) +
  geom_point(data=kde.points, aes(x=x, y=y), color = "lightgray") + 
  labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
  coord_sf(xlim = c(821900, 817000),  ylim = c(1605000, 1615000))
kde.plot + annotate("text", x = 818000, y = 1614500, label = paste("95% ",round(ver95$area,2),"ha"), 
                    fontface = "bold", size = 4, color = "white")

The anatomy of the function above is as follows:

• kernelUD(), estimation of kernel home-range
• getverticeshr(), extract home-range contour

The rest of the information can be found by examining the help menu for the various commands or looking at the display options for each portion of the script.

3.2.3 Automation Using Functions

Depending on your desired results, the steps above can be included in a function that will allow you to automate this process over a list of files and create mcp or kde analyses for all the individuals separately. Basically, we are taking all of the elements of the scripts above, making some of the values generic (e.g. filename) and creating the steps for the process that can be repeated over a list of files that you created earlier in the exercises. The object called files contains a list of all the *.csv files for each individual. So we will apply the function, over the list, to automate the analysis and visualizations for each individual. Because the kde process is essentially the same as above, only portions of the function that were not previously described will be discussed.

#run KDE creating plot for each individual with a progress bar
kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs"))
  xy.sp <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy.sp, h="href", kern="bivnorm", grid=100)
  ver95 <- getverticeshr(kde, 95)
  ver85 <- getverticeshr(kde, 85)
  ver75 <- getverticeshr(kde, 75)
  kde.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde95.sf <- st_as_sf(ver95)
  kde85.sf <- st_as_sf(ver85)
  kde75.sf <- st_as_sf(ver75)
  kde.plot <- ggplot() + theme_bw() + 
    theme(axis.text.y = element_text(angle = 90, hjust = 0.5, vjust = 0.5)) + 
    geom_spatial_rgb(data = thailand, aes(x=x, y=y, r=red, g=green, b=blue)) +
    geom_sf(data=kde95.sf, fill = "blue", alpha = 0.4) +
    geom_sf(data=kde85.sf, fill = "green", alpha = 0.3) +
    geom_sf(data=kde75.sf, fill = "yellow", alpha = 0.2) +
    geom_point(data=kde.points, aes(x=x, y=y), color = "lightgray") + 
    coord_sf(xlim = c(822000, 813950),  ylim = c(1616250, 1604500)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5))
  kde.plot + annotate("text", x = 815000, y = 1616000, label = paste("95% ",round(ver95$area,2)," ha"), 
                      fontface = "bold", size = 5, color = "white")
}

pblapply(files, kde_raster)

  |                                                  | 0 % ~calculating  
  |+++++++++++++++++++++++++                         | 50% ~00s          
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  
[[1]]

[[2]]

KDE Analysis OPHA1	KDE Analysis OPHA2

The anatomy of the script above is as follows: the kde_raster <- function(filename){... function runs for each file listed from the project directory in the files list complete the following commands

• read.csv(file = filename), there is no specific file name needed as it will use each *.csv in the list as a new file
• pbapply(), runs the function kde_raster over the list of *.csv in the files object

The rest of the information can be found by examining the help menu for the various commands or looking at the display options for each portion of the script.

Try it yourself! Using the example script above, and following the example functions, create a looping function with pblapply to run the mcp analysis for each individual in the files list.

3.3 Continuous-Time Movement Modeling

There is a more powerful package called ctmm that can be used for the creation of conventional kernel-density estimates as well as uniformly weighted kernel-density estimates. While I am less familiar with the variables, objects, and setting for each model, I was able to follow along with the vignette for Autocorrelated Kernel Density Estimation on the package details to create a similar kde analysis as above. All of the same methods automation methods can be employed, however, there are some difference to note.

For this portion of the exercises we will prepare the two OPHA *.csv files for analysis using akde in the `ctmm package.

## convert movebank to telemetry data
opha.telem <- as.telemetry(data, projection = "+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs")
opha1.telem <- as.telemetry('OPHA1.csv', projection = "+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs")
opha2.telem <- as.telemetry("OPHA2.csv", projection = "+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs")

You can examine the data in the environment to see how the data was restructured to be available for further analysis. Notice that a projection type is listed for each object. This was used to ensure the data would be geographically located rather than placed in general space by meters.

# calculate fit object
opha2.guess <- ctmm.guess(opha2.telem, interactive=FALSE)

# compute akde and get summary information
opha2.ouf <- ctmm.fit(opha2.telem, opha2.guess)
wAKDE.opha2 <- akde(opha2.telem, opha2.ouf, weights=TRUE)

summary(wAKDE.opha2)
# basic plot
plot(opha2.telem, UD=wAKDE.opha2)

You can see the details of each of the object created in your environment. The analysis create a class UD object that cannot be used in ggplot without conversion. Additionally, we will create an object specifically to use for the inclusion of area on the final plot.

# convert to SPDF to convert to sf object for plotting in ggplot
ud.sp <- SpatialPolygonsDataFrame.UD(wAKDE.opha2)
ud.sf <- st_as_sf(ud.sp)

opha2.sf <- st_as_sf(opha2, coords = c("utm.easting","utm.northing")) %>% st_set_crs(32647)

# adke summary data
akde.summary <- summary(wAKDE.opha2)
akde.summary[["CI"]][3]

[1] 124.8182

# ggplot for akde
ctmm.akde <- ggplot() + geom_spatial_rgb(data = thailand.full, aes(x=x, y=y, r=red, g=green, b=blue)) +
  geom_sf(data = ud.sf[3,], fill = "blue", alpha = 0.4) + 
  geom_sf(data = ud.sf[2,], fill = "green", alpha = 0.3) + 
  geom_sf(data = ud.sf[1,], fill = "yellow", alpha = 0.2) +
  geom_sf(data = opha2.sf, color = "lightgray") + 
  coord_sf(xlim = c(825500, 813500),  ylim = c(1618000, 1602000)) +
  labs(x="Longitude", y="Latitude", title = opha2.sf$individual.local.identifier) +
  theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) +
  theme_bw() + theme(axis.text.y = element_text(angle = 90, hjust = 0.5, vjust = 0.5))
ctmm.akde + annotate("text", x = 815000, y = 1617000, 
                     label = paste("Area ",round(akde.summary[["CI"]][3],2)," sqkm"),
                     fontface = "bold", size = 5, color = "white")

The details of this analysis are pretty new to me so I don’t yet have a solid grasp on the various components of the script, but if you want more information you can view:

Edelkind, J.L., Stalker, J.B., Beck, D.D., Denardo, D.F., Emplidge, P.G., Gallardo, L., Gentry, C.M., Goode, M.J., Jennings, R.D., Jones. J.L., Kwiatkowski, M.A., Nowak, E.M., Repp, R.A., Schuett, G.W., Sullivan, B.K., Tracy, C.R., and Gienger, C.M. (2025) Environmental Influences on Space Use of Gila Monsters (Heloderma suspectum). Herpetologica. https://doi.org/10.1655/Herpetologica-D-24-00053

This was simply an example of how to run the analysis and prepare the resulting objects for visualization.

4 YOUR TURN!

Now it’s your turn! Although some of this might not be applicable to your thesis research, try this information out on a dataset of your choice from any of the sources above. If you can apply this to your thesis, then add it to your website! Otherwise create a new repository for your presentation in the next class.