Using R to Calculate KDE Home Ranges

One of the most common methods of calculating home ranges is to use kernel density estimates. While there are a multitude of different methods to do so, I’ve provided an R script below to help you out. The only thing you’ll need to begin is a CSV file of coordinates. While coordinates should be taken in UTM, its not a problem if its in lat long. There’s a line of code in the script that will transform it. Due to formatting issues, your lat long coordinates should be in decimal form.

Please note, you’ll need to know the epsg codes for your specific study area. EPSG codes can be referenced here. This also isn’t the only way to do it. Its one way to do it. When I wrote this script, I didn’t have that much experience with writing R scripts, and was in the midst of an undergraduate project. I suspect that if you’re here, you’re in the same boat.

Before I provide the script, I just want to let you know that there is another option available for you. You can always use Hawthorne Beyer’s geospatial modeling environment that can be downloaded from here. For his tool though, you’ll need to know a little bit about GIS and have ArcGIS installed. His tool hasn’t been updated for the most recent versions of ArcGIS, so hopefully you have an older version installed if you opt for that route. The last time I checked, GME was only available for ArcGIS versions up to 10.3. I believe the earliest version of ArcGIS that you can download now is v. 10.4. I contacted ESRI customer support and they said that they were unable to provide an earlier version. However, if you have a current license for their software, you can use that license for an earlier version should you somehow obtain it.

Required Software: R

Recommended Software: R Studio, either QGIS or ArcGIS

If you haven’t already, you’ll need to install R. I recommend R Studio, which is an IDE. For me, it just makes using R a lot easier. While technically you don’t need QGIS or ArcGIS, I do recommend them for making maps. QGIS is open source whereas ArcGIS requires a license. If you’re a student, there’s a chance that your school can provide you with a student evaluation license that’s good for a year. There is a learning curve with GIS software, so I only recommend it if you know how to use the software.

Here’s the script. You can copy and paste it into your R console. Please read it through, as I write notes in areas that you need to edit for your own specific work.

> rm(list=ls())

> library(ks)
> library(sp)
> library(rgdal)
> library(raster)
> library(spatialEco)
> library(rgeos)

##if these packages aren’t installed, you’ll need to run install.packages(“”), e.g. if you need the ks package, you’ll type in the function >install.packages(“ks”)

> coords <- read.csv(“”, header = TRUE) ##in between the quotation marks, you’ll need to enter the filename of your CSV which contains your coordinates

> coords.xy <- coords[c(“long”, “lat”)] ##change long and lat to whatever it is you named your columns containing your coordinates, please note that longitude corresponds to x and latitude corresponds to y
> coords.spatialpoints <- SpatialPoints(coords.xy) > proj4string(coords.spatialpoints) = CRS(“+init=epsg:4326) ##this is only necessary if your original coordinates are in latitude and longitude format

> coords.utm <- spTransform(coords.spatialpoints, CRS(“+init=epsg:12345) ##this transforms your coordinates into UTM format, which is what will allow you to calculate your home range; you’ll need to look up the epsg for your specific UTM zone of interest
> df <- data.frame(coords.utm)
> coords.kde <- kde(df) > coords.raster <- raster(coords.kde, crs = “+init=epsg:”)
> coords.95 <- raster.vol(coords.raster, p=.95) ##this is your 95% home range estimate, if you want another isopleth, just change your p value, e.g. 50% home range estimate will be p=.5 instead of p=.95
> coords.polygon <- rasterToPolygons(coords.95, dissolve = TRUE)
> area(coords.polygon) ##this will tell you the area, in the event that more than one area is returned, you can use the function sum(area(coords.polygon)) to get the total area

##if you want to convert from m^2 to km^2, the conversion factor is multiplying by 1e-6.

##end script

The featured image is of a kernel density estimate with some artistic edits done in ArcGIS, v. 10.4.

Creative Commons License
This work by Teng Keng Vang is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s