Kernel home ranges constrained by boundaries

James E Paterson

Kernel density estimators are commonly used for estimating the home range of animals (see my previous posts on estimating home ranges with kernels).

But, a common issue is that the kernel contour for an animal extends beyond a boundary the animal cannot cross.

Can we adjust the density estimator and limit kernel contours to not cross boundaries?

The code and data I use in this tutorial are available on my on my GitHub repo.

Some common examples of boundaries include:

  • steep cliffs,
  • fences,
  • shorelines (e.g. a fish can’t extend into terrestrial habitat)

In these cases, the estimated home range includes areas that are not used by the animal.

Why do boundaries create biased estimates? Because the contour (e.g. 95% contour line) smooths pixels in the grid that include both ‘used’ (inside boundary) and ‘unused’ (outside boundary) areas. This results in the contour lying beyond a boundary the animal can’t cross.

The clever solution used by Benhamou and Cornélis (2010) was to use ‘mirrors’ of points close to the boundary. Every point that is in a pixel that includes the boundary is reflected on to the other side of the boundary before estimating the kernel. The result is a contour lying directly on the boundary. This approach is incorporated into the kernelUD function of the adehabitatHR package, and is what we’ll use to work through an example.

Before we start, I’m going to warn you this won’t work in every situation because the boundary can’t be too tortuous (curvey). You might be able to get around this by simplifying the boundary and reducing the number of vertices.

For the this approach these two conditions must be met:

  • each segment of your boundary must be >= 3 x h (h = smoothing factor)
  • the angle between segments in boundary must be > \(\pi\)/2 or < -\(\pi\)/2

First, load the libraries we need.

library(sp) # Spatial objects
library(rgdal) # If you want to read boundary as a shapefile
library(adehabitatHR) # For kernels
library(dplyr) # For filtering data
library(magrittr) # For piping

Next, let’s try out boundary-limited kernels using a sample data set tracking_sample.csv with relocations of 5 turtles. For the example, we will only use two individuals (T003 and T004).

The two easiest ways to identify the boundary are:

  • read in a shapefile with the boundary as a line, or
  • construct the line in R using the coordinates of the vertices
turtles <- read.csv(file = "tracking_sample.csv",
                    stringsAsFactors = FALSE) %>%
  filter(id %in% c("T003", "T004"), # we will only use 2 individuals in example (T003 and T004)
         !is.na(x), # Remove NA coordinates
         !is.na(y)) 

# Make turtles a Spatial Point Dataframe
turtles.sp <- turtles
coordinates(turtles.sp) <- c("x", "y")

# Set CRS
turtles.sp@proj4string <- CRS("+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

# You can import boundary line as a shapefile
bound <- rgdal::readOGR(dsn = "boundaryfolder", # Files in "boundaryfolder"
                           layer = "boundary_shape", # name of shapefile
                           encoding = "ESRI Shapefile", # It is an ESRI shapefile
                        verbose = FALSE) # Change to TRUE and it prints source and number of features

# Or, you can create a boundary line with a vector of coordinates to form vertices (commented here)
# bound <- structure(list(x = c(349801.0, 350097.4, 350323.9, 350143.9),
#                            y = c(4944954, 4945324, 4945503, 4946280)), .Names = c("x", "y"))
# bound <- do.call("cbind",bound)
# Slo1 <- Line(bound)
# Sli1 <- Lines(list(Slo1), ID="myboundary")
# bound <- SpatialLines(list(Sli1))

# Make sure bound CRS same as point CRS
bound@proj4string <- turtles.sp@proj4string

As you can see, many points fall close to the boundary (T003 on the left site, and T004 on the right side). How does this affect our home range estimate?

# Make normal kernel
trad.kernel <- kernelUD(turtles.sp[,"id"], # Only pass the "id" column and the coordinates
                        h = "href", # Use the reference bandwidth
                        grid = 1000)

trad.kernel.poly <- trad.kernel %>%
  getverticeshr(., percent = 95) # Get the vertices for the 95% contour line

# Plot
plot(bound, 
     col = "blue",
     lwd = 5)
plot(trad.kernel.poly,
     add = TRUE)
plot(turtles.sp,
     col = as.factor(turtles.sp$id),
     pch = 16,
     add = TRUE)

It is clear our home range estimates include a bunch of area that the animals cannot access. Why does this matter? Well, if you’re studying:

  • habitat selection, you don’t want to include unused patches as used
  • social behaviour or territoriality, you don’t want to incorrectly assume home ranges are overlapping when they don’t
  • any number of other reasons related to home range size, shape, or overlap!

Now, let’s re-estimate the kernels by taking into account the boundary. Since I want to use the reference bandwidth (“href”) for each animal, I calculate the boundary limited kernels separately. The h argument I use is based on the “href” from the unbounded kernel above.

# Make a kernel with boundary for T003
bound.kernel.t003 <- turtles.sp[turtles.sp$id == "T003","id"] %>%
                              kernelUD(.,
                                       # h must be numeric 
                                       # here I give the numeric argument for "href" from our estimate above
                                       h = trad.kernel$T003@h$h, 
                                       # If you don't set grid, the default often won't follow the boundary
                                        # I set grid to have 1000 x 1000 pixels
                                       # You may have to adjust (will depend on your scale)
                                       grid = 1000, 
                                        # The boundary argument takes SpatialLines objects
                                       boundary = bound) 

# Get the vertices for T003 and extract as a polygon
bound.kernel.t003.poly <- bound.kernel.t003 %>%
  getverticeshr(., percent = 95)

# Make a kernel with boundary for T004
bound.kernel.t004 <- turtles.sp[turtles.sp$id == "T004","id"] %>%
  kernelUD(.,
           h = trad.kernel$T004@h$h,
           grid = 1000,
           boundary = bound) 

# Get vertices for T004
bound.kernel.t004.poly <- bound.kernel.t004 %>%
  getverticeshr(., percent = 95)

Now let’s have a look to see if the kernels respected the line boundary we set. The new kernels are filled in with gray, and the original kernels are filled in white.

You’ll notice that this correction results in a 95% contour that is NOT the same as just clipping the polygon against the boundary. That’s because it changes the underlying kernel utilization distribution by adding those mirrored points.

I hope this is helpful and, good luck applying it to your own data!