Learning Outcomes

After completing this lab, you should be able to:

Your tasks

Dataset Source

Data on the green python (Morelia viridis) in Cape York, was downloaded from https://www.movebank.org/cms/webapp?gwt_fragment=page=studies

INTRODUCTION

Animal movement analysis for tracking animals, seeking covariates with landscapes, and determining their range has received considerable attention over the last few decades. Not surprisingly, with the improvement in ways to track animals, such as smaller and more robust GPS collars, a greater number of studies on a greater variety of animals has followed. In my opinion there is little doubt that this will continue to exponentiate, and the spatial experts will be right there helping ecologists explore the movement processes of the Great White Shark, migratory birds, and Galapagos tortoises. What a remarkably interesting career that would be!

Movement data are just data that tracks the change of an object through time. To get you thinking, what about the movement patterns of burglars? Do they stick to the same patch or randomly burgle? Are they territorial? Some interesting work has been done in this area using the same principles as applied to animals, where the patterns of offences are analysed by different levels of offenders – from the beginner to the advanced. See Johnson (2014) for a description of the “offender as a forager”. More topical – how about tracking the Coronavirus? Contact tracing produces a similar dataset to that which we will use in this lab – and LOTS of it! The question to ask is, where did Person X spend most of his/her time whilst potentially shedding the virus. Those people may need to isolate.

The obvious difference to anything we have studied so far in this unit is that we are no longer considering individual points as individual events. Instead, we now have a point, which is connected by another point, which is connected by more points that are all related to one individual. These points are usually generated at regular intervals, such as every 30 minutes, once a day, once a week etc. Furthermore, often more than one individual is tagged, but rarely will you have data at the population level, so it is important to temper analyses to the sample you have access to. Analysis tools and techniques need to adapt to these differences – one individual animal or a sample of individuals? Whilst analysing multiple individuals at one time is possible, it is more common to analyse each with the same tool several times. Here, we will just look at one individual to give you the skills to repeat it and if desired find patterns you are able to generalise.

Animal tracking data is usually collected by one or both of two major methods – the Lagrangian and the Eulerian methods. The distinguishing factor between them is that under the Lagrangian method individuals are tagged using a radio tracker, GPS or special bands, which ping a location every so often. Under the Eulerian methods, many locations are continuously surveyed by camera traps, which are triggered if an animal walks past them; they are therefore limited to the observation arena they are set up for. Their analyses are different and I’ll leave it to you to research more about that. Here we will focus only on data collected with the Lagrangian method. If you ever have camera data at your disposal, have a look at camtrapR (Niedball et al., 2016).

EXERCISE 6.1 – GETTING STARTED

Start ArcGIS PRO and create a new map by clicking on the Map template icon. You will need to name your project and set a working folder location. Enter lab_06_pythonmovement for the name and then choose an output folder (e.g. Lab6_Data).

Remember to click on View>Catalog Pane to see your data sources and Insert> Add Folder if the one you want is not there. Browse to the folder where you placed the tutorial data off the Blackboard and connect to the EX61 folder. Add in the point and line data from this folder.

Right click on the point layer and symbolise Unique Values according to tag_ident:

As the basemap, choose Imagery Hybrid under Map and then just explore the habitat of these pythons. There’s clearly a couple of them that have ventured outside of their tree canopies and into the open spaces but in general, they have well-defined smallish ranges, relative to e.g. a migratory bird. Not being the author of the data, I currently assume that the association with dirt tracks to be a sampling artefact (“convenience sampling”) rather than a habitat preference.

Choose Symbology again and this time choose Heat Map. Leave all other values as default.

Zoom to layer and you note that most of the “heat” is around Iron Range. One of the individuals found in that cluster is tagged number “2”. I’m going to restrict our analyses to this number 2…

EXERCISE 6.2 – INSTALL R & R-STUDIO

If you are working in the Spatial labs, on Spatial computers, you can skip to the next section. Otherwise, navigate to https://cran.rstudio.com/ to download R. Under Download and Install R, choose the link for your OS (e.g. Download R for Windows). At this point we just want R, so choose “base” at the next screen. Then download the latest version of R (4.x.x for Windows). Run the executable and accept all defaults.

Now, navigate to https://rstudio.com/products/rstudio/ to download R-Studio. Choose the download RSTUDIO DESKTOP button from the Open Source Edition. Then click this

Now, install R-Studio by accepting all and clicking next and Finish. That’s it. Now Run R-Studio.

EXERCISE 6.3 – READING IN DATA

The data we received was in both shapefile and csv format. I have put the csv format data in the EX63 folder.

Open the CSV labelled Green_python2 in Excel just to inspect it. You will notice an event id, a timestamp, lat and long data, sensor type, tag id, an individual id, and the name of the study.

Go back to your R-Studio session now. The first thing we need to do is install a package. The easiest way to do this is click Packages, then install and then type sf (for simple features).

Click Install. After it runs through its installation procedure, find it in the list of packages and make sure it is ticked on. Note: Your list of packages may be different to mine.

You should generally set up a working directory for any project, especially if you wish to export some of the outputs.

Please adapt the path and change the wd using setwd as shown in the following example:

wd <- setwd("e:/AGIS_Todd/lab6/Lab6/EX63")

#If this does not show up in your values correctly, run the line again until it does. For some reason, it often takes two goes.  

#Now, read in the file. It should say pythons, 251 obs of 9 variables. 

python <- read.csv("Green_python2.csv", as.is=TRUE)

How does it look? Have a look at the first few rows with the head command.

head(python)  
##      event_id visible               timestamp     long       lat
## 1 13077481827    TRUE 2003-01-09 13:17:00.000 143.2838 -12.74287
## 2 13077481828    TRUE 2003-01-14 02:07:00.000 143.2838 -12.74289
## 3 13077481829    TRUE 2003-01-16 05:00:00.000 143.2837 -12.74297
## 4 13077481830    TRUE 2003-01-17 01:10:00.000 143.2835 -12.74294
## 5 13077481831    TRUE 2003-01-18 00:18:00.000 143.2836 -12.74283
## 6 13077481832    TRUE 2003-01-19 00:48:00.000 143.2835 -12.74259
##          sensor_typ  tag_local_ individ
## 1 radio-transmitter #1006146617       2
## 2 radio-transmitter #1006146617       2
## 3 radio-transmitter #1006146617       2
## 4 radio-transmitter #1006146617       2
## 5 radio-transmitter #1006146617       2
## 6 radio-transmitter #1006146617       2
##                                          study_name
## 1 Green python Morelia viridis, Cape York Australia
## 2 Green python Morelia viridis, Cape York Australia
## 3 Green python Morelia viridis, Cape York Australia
## 4 Green python Morelia viridis, Cape York Australia
## 5 Green python Morelia viridis, Cape York Australia
## 6 Green python Morelia viridis, Cape York Australia

Currently our timestamp is a character string that cannot yet be used to calculate movement. We need to convert it into a machine-readable time class known as POSIX – standing for Portable Operating System Interface format. This just requires unscrambling the timestamp into its major components: year, month, day and so on.

To learn more about this Date-Time class (or any command you want to try) use the ?command syntax:

?POSIXct

This will load the help description for POSIXct. Have a quick read of that, especially look at the examples of code at the bottom. Frequently, you can plug into those quite easily. Now make a new variable called pythontime in this format and a universal time zone.

pythontime <- as.POSIXct(python$timestamp, format="%Y-%m-%d %H:%M:%S", tz="UTC")

Now we will create a move object defined by the move package. Please use Packages and install “move”.

Once installed check that it is ticked in the package list. Copy ALL of the code below and paste into R-Studio.

library(move)
## Loading required package: geosphere
## Loading required package: sp
## Loading required package: raster
mypython <- move( 
  #set coordinates 
  x=python$long, y=python$lat, 
  #provide my time format in POSIX 
  time = pythontime, 
  #set the projection – this is PROJ4 format. See below for more details. 
  proj = "+proj=longlat +ellps=GRS80 +no_defs", 
  #set the id 
  animal=as.numeric(python$individ), 
  #add the rest of the data 
  data = python)

Now we will project our data into metres. If you do not know which UTM zone your projection is in then take your longitude, add 180 and divide by 6 and round to the next integer (ceiling). This will give you 53.8. Hence, our data is in UTM Zone 54. To make a projection, first you need to define the current coordinate system. We have done that above with the line +proj=longlat + ellps = GRS80 + no_defs. That means we can now project it to something else. Both defining a projection and doing the projection requires some basic knowledge of the PROJ4 syntax. This is my go to for that https://mgimond.github.io/Spatial/coordinate-systems-in-r.html.

Essentially, you need to specify that it is a UTM projection, the zone, south or north, the datum, and the units. No_defs means no defaults, i.e., nothing other than what we prescribe below.

mypython_prj <- spTransform(mypython, "+proj=utm +zone=54 +south +datum=WGS84 +units=m +no_defs")

Exercise 6.3.1 – Summary Statistics and Cleaning

In this next stage of the analysis, we are going to run our data through some simple statistics: distance and speed. The idea of this is twofold. First, you get to make inferences about your animal, such as whether it likes to hang around in the same spot (like a snake), or make large distance runs (like a salmon). Second, you can screen your data for unusual observations, i.e., outliers. Try out the following commands. I think you will find that our snake doesn’t really travel great distances, and never at great speed.

#distance, in m, between points.  
distance <- distance(mypython_prj) 
hist(distance) 

#speed is in m per second. 
speed <- speed(mypython_prj) 
hist(speed)

After looking at the data in this way, I do not think there are any obvious errors. We shall proceed as is. Another cleaning routine is to ensure your animal is not found in any locations it cannot be – e.g. a pig in deep ocean, a fish on land. We essentially did that when we viewed the data in ArcGIS.

EXERCISE 6.4 – THE HOME RANGE

The basic definition of an animal’s home range is the spatial extent of its movements that are required in its normal activities to maintain its needs, which include finding a meal, a partner, sleeping, residing and caring for young. It does not include migration and should seek to eliminate random “walks” – unusual and brief excursions that rarely, if ever, repeat.

We will begin with a minimum convex polygon (MCP) method using the adehabitatHR package. Install it and make sure it is checked. By ticking this on, it will tick on a few other dependencies. That’s fine.

library(adehabitatHR)
## Loading required package: ade4
## Loading required package: adehabitatMA
## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp
## 
## Attaching package: 'adehabitatMA'
## The following object is masked from 'package:raster':
## 
##     buffer
## Loading required package: adehabitatLT
## 
## Attaching package: 'adehabitatLT'
## The following object is masked from 'package:move':
## 
##     burst

The MCP method encloses the outer locations to make a home range polygon. Usually, around 90% (it’s up to the user) of the points are retained to remove random walks.

First, we convert to a SpatialPoints class and then we use the mcp command:

mypython_pts <- as(mypython_prj, "SpatialPoints") 
mypython_pts.mcp <- mcp(mypython_pts, percent = 90, unin = "m", unout = "ha") 
mypython_pts.mcp 
## Object of class "SpatialPolygonsDataFrame" (package sp):
## 
## Number of SpatialPolygons:  1
## 
## Variables measured:
##   id     area
## a  a 6.055367

So, it would appear, our snake has a home range of about 6 ha, according to the mcp method. It would certainly be interesting now, to see how that compares to the other snakes that have been monitored. However, for now, let’s view the home range on a Sentinel image. Install and check the package named “raster”.

library(raster)
#Change this to suit your path! 
wd <- setwd("e:/AGIS_Todd/lab6/Lab6/EX64")
TCC <- stack('TCC.tif')
plotRGB(TCC, 3, 2, 1, stretch = "lin")
points(mypython_prj, pch = 20, col = "yellow")
plot(mypython_pts.mcp, add = TRUE, border = "red", lwd = 3)

Our final MCP appears below. This is the approximate home range of python number 2. The points to the west are a random foray and therefore not part of its home range using this method.

Let’s save that as a shapefile so we can use it in ArcGIS. We can use the sf package. No need to reinstall.

library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
#We first convert our spatialpolygonsdataframe to an sf object. 
mcp <- st_as_sf(mypython_pts.mcp) 
st_write(mcp, dsn = "mcp_py2.shp", layer = "mcp_py2.shp", driver = "ESRI Shapefile", append = FALSE) 
## Deleting layer `mcp_py2' using driver `ESRI Shapefile'
## Writing layer `mcp_py2' to data source `mcp_py2.shp' using driver `ESRI Shapefile'
## Writing 1 features with 2 fields and geometry type Polygon.

One of the drawbacks of the MCP approach is that it is sensitive to sample size. In addition, the 90% threshold is just a rule of thumb to deal with outliers, but it may not exclude them all, or it may exclude too many. For example, one extreme point left in the solution that is still in the top 90% could inflate the convex polygon by many ha, especially if you are dealing with an animal that moves quickly. Therefore, another technique is to model the utilization distribution. In other words, the more an animal is found in the same spot, or close to it, the more it is likely to be suitable habitat and part of its home range. This concept is the same as the very quick heat maps we produced in ArcGIS earlier. It uses a moving kernel to estimate the utilisation of space. Here we will use the default “href” for the smoothing parameter, but you can change it if you like. The higher you make it, the more smoothing will occur and vice versa.

mypython_pts.kernel <- kernelUD(mypython_pts, h = "href") 
mypython_pts.kernelarea <- kernel.area(mypython_pts.kernel, unin = "m", unout = "ha") 
mypython_pts.kernelarea 
##        20        25        30        35        40        45        50        55 
##  1.270790  1.588488  2.117984  2.541580  2.965177  3.494673  4.130068  4.765463 
##        60        65        70        75        80        85        90        95 
##  5.400858  6.036253  6.883447  7.836539  8.895531 10.166321 12.072506 15.884877

The last function will output the different areas for different levels of utilisation. I’ve chosen to output a kernel of 65% and of 90%.

plotRGB(TCC, 3, 2, 1, stretch = "lin") 
points(mypython_prj, pch = 20, col = "yellow") 
plot(getverticeshr(mypython_pts.kernel, percent = 65), lwd = 3, border = "red", lwd = 3, add = TRUE)  
plot(getverticeshr(mypython_pts.kernel, percent = 90), lwd = 3, border = "grey", lwd = 3, add = TRUE)

The output of the two scenarios is found below.

For comparison sake, let’s output the second one as a shapefile where percent = 90.

kernel90 <- getverticeshr(mypython_pts.kernel, percent = 90) 
kernel90 <- st_as_sf(kernel90) 
st_write(kernel90, dsn = "kern90_py2.shp", layer = "kern90_py2.shp", driver = "ESRI Shapefile", append = FALSE) 
## Deleting layer `kern90_py2' using driver `ESRI Shapefile'
## Writing layer `kern90_py2' to data source 
##   `kern90_py2.shp' using driver `ESRI Shapefile'
## Writing 1 features with 2 fields and geometry type Polygon.

Now bring that into ArcGIS and compare to the MCP.

Neither of the two approaches examined so far are really taking advantage of the TIME factor, which makes this dataset so unique and powerful. Brownian bridge movement models do incorporate the temporal (time) structure of the data by integrating the order of locations and the amount of time between them. In effect, they model the utilisation distribution based on the movement path, rather than just how tightly clustered the recorded points are, as with the previous two models. For example, a foray might be identified as locations that were not stopped at for long, followed by more locations of the same, and none of them are frequently, if ever, revisited. A good paper on this is one by Kranstauber et al. (2012) who also go into the mathematics required to do it. They describe how the movement path between A and B are modelled by applying a conditional random walk, which considers how much variance can be expected when travelling in a straight line from A to B – the variance of the Brownian motion.

They extend the concept to make it more “dynamic” using a sliding window approach where the variance of the Brownian motion is computed iteratively to account for an animal moving in different ways such as when foraging, travelling between sites and periods of rest. For example, a migratory bird will move over a small range when breeding, but over a long range when migrating. Thus, one value for the Brownian motion parameter is insufficient at estimating the whole trajectory, underestimating some paths, and overestimating others. Underestimation leads to utilisation distributions that are too narrow and overestimation to wider utilisation distribution areas.

The parameters required to compute a dynamic Brownian bridge are:

Object - a move object, hence mypython_prj can be used.

Raster - is the resolution of the output probability surface. 10 m corresponds well to the Sentinel imagery, so I have chosen that, but it need not be.

Location error – Potential amount of error from the sender/receiver system. I have guessed at 19 m. A bit more than a pixel.

Margin - The margin at the start and end of the moving window to avoid edge effects.

Window size – The size of the moving window along the track. The larger the window, the more stable the estimates are of the Brownian motion variance, but poorer the ability to capture frequent behaviour changes. The window size must be odd.

myDBBM <- brownian.bridge.dyn(mypython_prj, raster = 10, location.error = 19, margin = 5, window.size = 21)  
## Computational size: 8.5e+08
plotRGB(TCC, 3, 2, 1, stretch = "lin") 
plot(myDBBM, add = TRUE)

You may like to write it out as a tif and view in ArcGIS.

writeRaster(myDBBM, filename="BBmodel.tif", options="INTERLEAVE=BAND", overwrite=TRUE)

Now, to export this as a polygon we have a little more work to do, since it’s currently a raster. Firstly, we need to convert the output to a probability, which we do using getVolumeUD. Then we need to convert that to contours, the contours to polygons, subset the 90% polygon and then export it.

UDmypy <- getVolumeUD(myDBBM) 
plot(UDmypy) 

UDmypy_conts <- rasterToContour(UDmypy, nlevels = 20) 
mycontours <- st_as_sf(UDmypy_conts) 
my_contours_poly <- st_polygonize(mycontours) 
my_contours_subset <- subset(my_contours_poly, level == 0.9) 
my_contours_dissolved <- st_union(my_contours_subset) 
st_write(my_contours_dissolved, dsn = "dbbm90_py2C.shp", layer = "dbbm90_py2C.shp ", driver = "ESRI Shapefile", append = FALSE, layer_options = "SHPT=POLYGON") 
## Deleting layer `dbbm90_py2C.shp ' failed
## Writing layer `dbbm90_py2C.shp ' to data source 
##   `dbbm90_py2C.shp' using driver `ESRI Shapefile'
## options:        SHPT=POLYGON 
## Writing 1 features with 0 fields and geometry type Geometry Collection.

Once again, bring that file into view in ArcGIS. I think you’d agree that this technique produces a more nuanced utilisation distribution.

So, there you have it. The points to the west are quite likely to be a result of foray and not where it spends a great deal of its time. The previous two methods also identified that but were less exclusive of the points in the centre and northwest, which were clustered but were not utilised for any great length of time. So, if we were looking to find Python2, say to refit a broken receiver, there appears to me that it has a few favourite trees in its relatively small range. Of course, there are obvious assumptions in the parameters chosen and it would be prudent to test other variations to determine their sensitivity in the solution. But, hopefully you get the idea!

References

Johnson, S. (2014). How do offenders choose where to offend? Perspectives from animal foraging. Legal and Criminological Psychology, 19(2), 193-210. doi: 10.1111/lcrp.12061

Kranstauber, B., Kays, R., LaPoint, S., Wikelski, M., & Safi, K. (2012). A dynamic Brownian bridge movement model to estimate utilization distributions for heterogeneous animal movement. Journal of Animal Ecology, 81(4), 738-746. doi: 10.1111/j.1365-2656.2012.01955.x

Niedballa, J., Sollmann, R., Courtiol, A., & Wilting, A. (2016). camtrapR: an R package for efficient camera trap data management. Methods In Ecology and Evolution, 7(12), 1457-1462. doi: 10.1111/2041-210x.12600

End of Session.