Update: The code for using the adehabitatHR package is given below. Beneath the code will follow my original blog post for informative purposes. Please note, if you didn’t do a transformation, I’ve included a line of code to do a transformation for you. If your coordinates are already in UTM or another projected coordinate, you can ignore line 7 (the line with the spTransform function). ver95 is the 95% contour and ver50 is the 50% contour.
library(adehabitatHR)
library(rgdal)
coords <- read.csv("ID5593.csv", header = TRUE)
coords.xy <- coords[c("x", "y")]
coords.spatialpoints <- SpatialPoints(coords.xy)
proj4string(coords.spatialpoints) = CRS("+init=epsg:4326")
coords.utm <- spTransform(coords.spatialpoints, CRS("+init=epsg:32614"))
ud <- kernelUD(coords.utm)
ver95 <- getverticeshr(ud, 95, unout = "km2")
ver50 <- getverticeshr(ud, 50, unout = "km2")
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)
##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 coords.xy coords.spatialpoints proj4string(coords.spatialpoints) = CRS(“+init=epsg:4326) ##this is only necessary if your original coordinates are in latitude and longitude format
> coords.utm df coords.kde coords.raster coords.95 coords.polygon 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.
This work by Teng Keng Vang is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
Teng Keng Vang 
When I run the script for 95%KDE, I get the exact same values as for 50%KDE. Can you please help me? I altered your R script a bit, but the overarching steps are the same. But see below:
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”)
##make sure there are no NA values or blank rows, otherwise it will not run
coords <- read.csv("E:/PortMansfieldRedSnapperMovements/PM_HomeRange_R/ID5593.csv", header = TRUE) ##chg fish ID here ##in between the quotation marks, you'll need to enter the filename of your CSV which contains your coordinates
coords.79 <-subset(coords, TimeNum==0.79)
coords.81 <-subset(coords, TimeNum==0.81)
coords.82 <-subset(coords, TimeNum==0.82)
coords.19 <- rbind(coords.79,coords.81,coords.82)
coords.xy <- coords.19[c("x", "y")] ##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:32614")) ##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)
coords.95 <- raster.vol(coords.raster, p=0.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=0.50
coords.polygon <- rasterToPolygons(coords.95, dissolve = TRUE)
coords.area <- 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
coords.sumarea <- sum(coords.area)
coords.sumarea
areatime.19 <- data.frame(hour=19,kde=coords.sumarea,typ=95)
areatime.19
##if you want to convert from m^2 to km^2, the conversion factor is multiplying by 1e-6.
Hi Catheline. There’s a line of code that says >coords.95 coords.50 <- raster.vol(coords.raster, p=.5)
So with the change in nomenclature, you'll need to update the coords.90 variable to coords.50 where ever it appears in the rest of the script, which is only one more time I believe. Its basically the next line of code.
Hope this helps.
Sorry about that Catheline, it looks like my comment got jumbled a little bit.
Look for the line of code that says >coords.95 <- raster.vol(coords.raster, p=0.95)
You need to change the p=.95 to p=.5
Thank you for your prompt reply.
If I understand correctly, you are suggesting I update that section with this:
coords.50 <- raster.vol(coords.raster, p=0.50) ##this is your 50% 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=0.50
coords.polygon <- rasterToPolygons(coords.50, dissolve = TRUE)
However I still end up with the same area.
Any other suggestions by chance?
How many data points are there in the data?
I have 275 rows and 2 columns (x and y)
Sure! What email should I send it to?
Hi Catheline, I sent you an email to the email account associated with your wordpress acount. Please check your .edu account.
Hello,
This is an extremely useful and well written guide. Thank you for providing it. I just have two questions;
1. I have over 50 animals I need to calculate KDE for, is it possible to tweak this code to recognise and calculate for multiple IDs within the same CSV?
2. How can I transform the end result into a shapefile for ArcGIS?
Very best,
Georgia
Hi Georgia
I’d probably just run a for loop to run it for multiple IDs. You can export your results using the writeOGR() function. The function is found in the gdal package. If you look at the bottom of this post (https://tengkengvang.com/2017/08/18/local-convex-hull-or-locoh-in-r/), I explain how to do it.
I haven’t had time to update the post, but I’d recommend using the adeHabitatHR package to run KDE home ranges instead. It requires only one line of code to do so (assuming your coords are UTM already).
Hello!
Thanks very much for this guide! I have gotten my KDE area values using the getverticshr function and am now trying to overlay the kernel density onto a map using the getNOAA.bathy() function and geom_contour() function. However, I notice that for the code to work, I need to have the option of selecting lat/lon columns from the saved kernel density value of choice (ie. kernel_poly if the example is kernel_poly <- getverticeshr(kernel,
percent = 90,
unin = "m",
unout = "km2")
However, after I use the getverticeshr() function the only option I have is the column of the animals identification "id" and the calculated "area" of KDE. Do you have any idea to why this is? Or any codes that you use to graph kernel density?
Thanks so much!!
Best,
Katherine
Hi Katherine,
If I had to guess, the problem is probably because you are trying to mix projected and unprojected data. The KDE values should’ve been projected into an equal area projection (I recommend UTM as its the standard in spatial ecology). However, when you call bathymetric data from NOAA, it looks like that data is unprojected, as it is in lat long. It is possible to do this all programmatically (I’m thinking something with raster$extent()), but you would have to experiment with that. You could adhoc it by plotting your points in something like QGIS and estimating the extents you would like to use.
As far as geom_contour is concerned, the input data has to be in vector format. getverticeshr returns a raster. You would need to polygonize the raster as we say in GIS. You might find this helpful https://r-spatial.github.io/stars/articles/stars5.html
I just saw your email address, and my first thought was, she must be in the Shaffer Lab, and lo and behold, there you are. Scott helped me out quite a bit when I was an undergrad at UCSC and I always appreciated it. I normally don’t do this, but out of respect to Dr. Shaffer, if you need additional help, please use the contact me form and we can discuss more in private.
Hi there,
Many thanks for this guide. I seem to be having some issues with the code throwing an error message when I use it. Where am I going wrong please?
> coords.utm df coords.kde coords.raster coords.95 coords.polygon area(coords.polygon)
Gives me:
Error: unexpected symbol in “coords.utm df”
Many thanks!
Elli
For whatever reason, the script at the bottom of my post had the syntax mangled by wordpress. You should be using the code at the very top of the post. The one at the bottom is just left for reference purposes.
The script at the top assumes a csv with point coordinates, one column x and one column y. X is your longitude and Y is your latitude. You read it in as a csv.