Introduction
The move package provides a series of functions to import, visualize and analyse animal movement data. With the move package movement data associated with individual animals or a group of animals can be imported as a single object. The minimum requirements are timestamps, coordinates, and a unique ID per animal. Further attributes associated with the individuals (e.g. sex, age, species, etc), the locations (e.g. temperature, behavior, battery voltage, etc.), coordinates (projection) and the study (e.g. citation, license terms) can be also included.
Movebank (www.movebank.org) is a free online database of animal tracking data, where data owners can manage their data and have the option to share it with colleagues or the public. Movebank users retain ownership of their data and choose whether and with whom to share it. If the public or a registered user has permission to see a study, the study can be downloaded as a .csv file and imported directly into R. Those with Movebank accounts can also log in, browse and download data they have access to from Movebank directly within the Move package (see the ”browsing Movebank” vignette for more information). Besides importing data from Movebank, it is also possible to use personal data sets, as long as they contain the minimum required information: timestamps, coordinates, and animal IDs.
A series of functions allow you to plot, summarize, and analyse the movement data. Individual Move objects can be stacked to a MoveStackobject, which includes a series of animals and the additional data that are associated with them. This allows you to work with multiple animals at the same time. Tracks can be also bursted by a specific characteristic, e.g.specific behaviors (migratory, non-migratory, resting behavior), which facilitates analyzing segments of the trajectory belonging to a specific behavior or variable. This package also provides functions to calculate the dynamic Brownian Bridge Movement Model (dBBMM) and the utilization distribution(UD) of a trajectory.
This vignette gives an overview of the main functions of the move library. For more functions and more details about functions we recommend viewing the corresponding help files.
Data import
The import of tracking data into R can happen in three ways:
- Importing data directly from Movebank using the web interface of the move package (with
getMovebankData(), getMovebankLocationData(), getMovebankNonLocationData()).
- Importing a .csv file downloaded from Movebank (or .zip files from EnvDATA) using the
move() function.
- Importing non-Movebank data using the
move() function.
In all cases multiple individuals can be imported simultaneously.
Import data directly from Movebank
See the “browsing Movebank” vignette for detailed explanation.
Import data downloaded from Movebank
Files from Movebank can be imported with move(x="path"), where “path” is a character string that locates the file on the hard drive. The files can also be compressed csv or .zip files downloaded with the EnvDATA tool. It may look like:
Import non-Movebank data
Any data set containing movement data (position and time information) can be imported as a move object. It is mandatory to indicate the columns that contain the coordinates and the timestamp column in the correct format. The projection of the coordinates should be indicated, the data frame that includes all of the imported data can be added (if it contains additional information associated to the locations) and if the data correspond to a single individual, its name/ID can be indicated as a character, but if the data correspond to several individuals the column with the animal names/IDs has to be indicated.
Any data set containing movement data (position and time information) can be imported as a move object. You must define the columns that contain the coordinates, the projection of the coordinates, and the format of the timestamp column. If the data correspond to a single individual, its name/ID can be indicated as a character, but if the data correspond to several individuals, the column with the animal names/IDs must be indicated. If the data frame contains additional columns, these additional attributes can be included.
In the following example a move object is created from a custom file read in as a data frame:
data <- read.csv(system.file("extdata","leroy.csv.gz",package="move"))
leroy <- move(x=data$location.long, y=data$location.lat,
time=as.POSIXct(data$timestamp, format="%Y-%m-%d %H:%M:%OS", tz="UTC"),
proj=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"),
data=data, animal=data$individual.local.identifier, sensor=data$sensor)
By entering the object an overview is provided. This overview indicates the object class (Move), number of coordinates, the extent of the coordinates, the coordinates reference system (projection), the number of columns of the imported data.frame, and the names of the columns of the data.frame, as well as the first and last timestamp and the duration of the observation.
## class : Move
## features : 919
## extent : -73.93067, -73.84366, 42.70898, 42.7687 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 17
## names : timestamp, eobs.battery.voltage, eobs.horizontal.accuracy.estimate, eobs.key.bin.checksum, eobs.speed.accuracy.estimate, eobs.start.timestamp, eobs.status, eobs.temperature, eobs.type.of.fix, eobs.used.time.to.get.fix, ground.speed, heading, height.above.ellipsoid, utm.easting, utm.northing, ...
## min values : 1234354605, 3596, 3.07, 3258904, 0.27, 2009-02-11 12:14:59.000, A, 13, 3, 4, 0.01, 0, -169.6, 587507.837877134, 4729143.16566605, ...
## max values : 1236158219.998, 3666, 97.02, 4291715164, 33.04, 2009-03-04 09:15:01.000, A, 35, 3, 119, 31.71, 359.79, 349, 594679.382402579, 4735720.47868847, ...
## timestamps : 2009-02-11 12:16:45 ... 2009-03-04 09:16:59 Time difference of 21 days (start ... end, duration)
## sensors : gps
## indiv. data : eobs.fix.battery.voltage, manually.marked.outlier, visible, sensor.type, individual.taxon.canonical.name, tag.local.identifier, individual.local.identifier, study.name, study.timezone
## indiv. value: 0 NA true gps Martes pennanti 74 Leroy Urban fisher GPS tracking Eastern Standard Time
## unused rec. : 1071
## study name : Urban fisher GPS tracking
## date created: 2020-06-24 15:04:52
Import movement objects of other packages
The following objects containing movement data can also be read in by the move() function:
- object of class
ltraj from the library adehabitatLT
- object of class
telemetry or a list of telemetry objects from the library ctmm
- object of class
track from the library bcpa
- object of class
track_xyt from the library amt
- object of class
binClstPath or binClstStck from the library EMbC
- object of class
data.frame containing data downloaded from Movebank webpage or with getMovebankLocationData
Handling multiple animals
If the data are from Movebank, the function will automatically recognize that there are multiple individuals contained in the file. If the custom dataset includes more than one animal, you must indicate the column with the animal IDs in the argument animal. When multiple individuals are imported, the move() function automatically creates a stack of Move objects called MoveStack. A MoveStack can also be created from a list of multiple Move or MoveStack objects with the function moveStack(). In this case all objects have to be in the same projection, and their timestamps have to be in the same time zone. Most functions of the Move package are capable of working with MoveStacks.
In the following example a MoveStack is created from two Move objects:
In the overview the information of all individuals is combined.
## class : MoveStack
## features : 9877
## extent : -73.94032, -73.84366, 42.70898, 42.851 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## variables : 20
## names : timestamp, eobs.battery.voltage, eobs.horizontal.accuracy.estimate, eobs.key.bin.checksum, eobs.speed.accuracy.estimate, eobs.start.timestamp, eobs.status, eobs.temperature, eobs.type.of.fix, eobs.used.time.to.get.fix, ground.speed, heading, height.above.ellipsoid, utm.easting, utm.northing, ...
## min values : 1234354605, 3449, 2.05, 524362, 0.22, 2009-02-11 12:14:59.000, A, -12, 3, 3, 0, 0, -315.6, 586587.171471437, 4729143.16566605, ...
## max values : 1270056686, 3740, 97.02, 4294919150, 48.67, 2010-03-31 17:30:02.000, A, 35, 3, 120, 43.12, 359.79, 436.6, 594679.382402579, 4744838.73355732, ...
## timestamps : 2009-02-11 12:16:45 ... 2010-03-31 17:31:26 Time difference of 413 days (start ... end, duration)
## sensors : gps
## indiv. data : eobs.fix.battery.voltage, manually.marked.outlier, visible, sensor.type, individual.taxon.canonical.name, tag.local.identifier, individual.local.identifier, study.name, study.timezone, behavioural.classification
## min ID Data : NA, NA, true, gps, Martes pennanti, 74, Leroy, Urban fisher GPS tracking, Eastern Standard Time, NA
## max ID Data : NA, NA, true, gps, Martes pennanti, 1016, Ricky.T, Urban fisher GPS tracking, Eastern Standard Time, NA
## individuals : Leroy, Ricky.T
## unused rec. : 2959
## date created: 2022-05-31 14:03:34
It is also possible to break down a MoveStack into a list of Move objects using the split() function. This is often useful to do calculations with functions that require lists as an input as e.g. lapply(). This list of Move objects can be transformed easily back into a MoveStack with the function moveStack().
Handling duplicated timestamps
Given that an animal cannot be at two places at the same time, Move* objects account for this and cannot be created if there are duplicated timestamps present in the data. There are several options to remove these duplicated timestamps:
If the study is from Movebank and it contains duplicated timestamps you can set the argument removeDuplicatedTimestamps=TRUE when reading in the .csv file with the move() function. This will retain the first of multiple records with the same animal ID and timestamp, and remove any subsequent duplicates. In some cases, one duplicate record might contain more complete information or a better location estimate than the other. In case you want to control which of the duplicate timestamps are kept and which are deleted, we recommend to download the data as a .csv file from Movebank or to use the function getMovebankLocationData(), find the duplicates using e.g. getDuplicatedTimestamps(), decide which of the duplicated timestamp to retain, and than create a move/moveStack object with the function {move(). Another option is to edit the records in movebank and mark the appropriate records as outliers.
One option is to remove the duplicated timestamps with the function duplicated(), nevertheless in this case, depending on the settings, the first or the last location will be retained. In some cases, one duplicate record might contain more complete information or a better location estimate than the other. In case you want to control which of the duplicate timestamps are kept and which are deleted, you can find the duplicates using e.g. getDuplicatedTimestamps(), decide which of the duplicated timestamp to retain, and than create a move/moveStack object with the function move().
Bursting a track
It can be interesting to compare different parts of the track with each other. For example, how do data points differ between winter and summer, or between behaviors like migrating and non-migrating, or between day and night. To indicate which point of the data set belongs to which category, a track is ’bursted’ with the function burst(). A track is bursted by supplying a vector with the specific behavior/variable associated with each location. The length of this vector has to be one shorter than the number of coordinates. The returned object belongs to the class MoveBurst.
In the following example, the trajectory of ‘leroy’ is bursted according to day and night.
The MoveBurst object has a very similar structure to a Move object. The MoveBurst object has an additional slot @burstId which contains the information of the assignment of each location to each behavior/category (see str(leroy.burst)). In the overview of the MoveBurst there is an additional information feature under “bursts” with the information how many locations fall within each category. MoveBurst objects cannot be stacked into a MoveStack.
The resulting MoveBurst can be plotted highlighting the different categories in different colors. Note that R by default only uses 8 different colors, if there are more than 8 categories colors will be recycled, in this case we recommended specifying the colors.
The plotBursts() function plots a circle at the midpoint of each burst segment (consecutive coordinates that belong to a single burst). The size of the cycles can have different meanings, depending on what is calculated in the sizeFUN argument (for details see ?plotBursts). The default size refers to the relative time of the burst segment compared to the whole track.
### in the default the size of the cicles is relative to the total time spent within each burst segment
plotBursts(leroy.burst,breaks=5, col=c("red", "black"), pch=19, add=F,main="Size of points: total time spent in burst segment", asp=1)
legend("bottomleft",legend=c("day","night"), col=c("red", "black"), pch=19)
## here, e.g., the size of the circles is relative to the total distance moved within each burst segment
plotBursts(object=leroy.burst, breaks=5, sizeFUN=function(x){sum(distance(x))},col=c("red", "black"), pch=19, add=F,main="Size of points: total distance moved in burst segment", asp=1)
legend("bottomleft",legend=c("day","night"), col=c("red", "black"), pch=19)
A MoveBurst object can be split into a list of Move objects, one per each burst segment:
leroyB.split <- split(leroy.burst)
par(mfrow=c(1,4))
plot(leroyB.split[[1]], type="b", pch=20, main=paste0(names(leroyB.split[1])," 1"))
plot(leroyB.split[[2]], type="b", pch=20, main=paste0(names(leroyB.split[2])," 1"))
plot(leroyB.split[[3]], type="b", pch=20, main=paste0(names(leroyB.split[3])," 2"))
plot(leroyB.split[[4]], type="b", pch=20, main=paste0(names(leroyB.split[4])," 2"))
Corridor behavior
The corridor() function identifies sections of the track that suggest corridor use behavior (i.e. parallel and fast movements). For details see LaPoint et al. (2013) Animal Behavior, Cost-based Corridor Models, and Real Corridors. Landscape Ecology. https://doi.org/10.1007/s10980-013-9910-0.
The speed threshold and how parallel the segments should be can be adjusted in the function. The proportion of speeds that are high enough to be a valid corridor point is be defined in the argument “speedProp”. The default value selects speeds that are greater than 75% of all speeds. The proportion of the circular variances low enough to be a valid corridor point is defined in the argument “circProp”. Low values of the circular variance indicate that the segments are (near) parallel. The default value selects variances that are lower than 25 % of all variances.
LeroyCorr <- corridor(leroy)
plot(LeroyCorr, type="l", xlab="Longitude", ylab="Latitude", col=c("black", "grey"), lwd=c(1,2))
legend("bottomleft", c("Corridor", "Non-corridor"), col=c("black", "grey"), lty=c(1,1), bty="n")
The resulting object is of class MoveBurst which in the @data slot contains new attributes associated to the locations, where “segMid_x” and “segMid_y” are the coordinates of the midpoint of each segment, and “speed”, “azimuth”, “pseudoAzimuth” and “circVar” are the variables used to identify the segments that suggest corridor behavior.
Dynamic Brownian Bridge Movement Model (dBBMM)
The dBBMM is a method to calculate the occurrence distribution of a given track, i.e. it estimates where the animal was located during the observed (tracking) period. It calculates the probability landscape for the transition between any two known consecutive locations given the amount of time it had available (assuming a conditional random (Brownian) motion between locations). It is “dynamic” because it allows the variance to change along the trajectory, using the behavioral change point analysis in a sliding window. For details see Kranstauber et. all (2012) A dynamic Brownian bridge movement model to estimate utilization distributions for heterogeneous animal movement. Journal of Animal Ecology. https://doi.org/10.1111/j.1365-2656.2012.01955.x
The Move* object used to calculate the dBBMM needs to be projected and, the location error, margin, window size and raster options need to be specified. Note that in most cases different values of window size and the margin do not have a great effect on the results. For details on the different ways to specify the raster options see the ‘Details’ section of ?brownian.bridge.dyn. In the following examples, the calculation is done on a raster with a pixel size of 100m:
The resulting object will be of class DBBMM, DBBMMStack, or DBBMMBurstStack depending on the input class. All of them contain a raster or a rasterStack were all pixels of each raster sum up to 1 (besides the DBBMMBurstStack, for details see ‘Calculating a dBBMM from a MoveBurst’ section)
Plotting the result provides limited visual information, as the values correspond to the probability the animal was present in a given pixel during the observed period. These results can be used to estimate probability of the animal being present in a certain area during the tracking period (see ‘Utilization distribution’ section), or to calculate the probability of the animal being present in a certain environment during the tracking period by overlaying the resulting raster over a categorical raster (e.g. of land cover) and summing up the dBBMM pixel values that fall within each category.
Utilization distribution (UD)
From the Dynamic Brownian Bridge object, the UD can be calculated, which represents the minimum area in which an animal has some specified probability of being located.
UDleroy <- getVolumeUD(dBB.leroy)
par(mfrow=c(1,2))
plot(UDleroy, main="UD")
## also a contour can be added
plot(UDleroy, main="UD and contour lines")
contour(UDleroy, levels=c(0.5, 0.95), add=TRUE, lwd=c(0.5, 0.5), lty=c(2,1))
The obtained UD can be subset to a specific probability:
Calculating a dBBMM from a MoveBurst
The dBBMM can also be calculated for a MoveBurst, for all burst types, or only for a subset of them. The result will provide one raster per burst segment in an object of class DBBMMBurstStack. By default it is calculated for all burst types:
The UD can also be calculated for each raster layer of the resulting DBBMMBurstStack. In this case to calculate the UD, each raster layer needs to be standardized so the sum of all pixels is 1. This can be done with the function UDStack()
To calculate the UD for each burst type, the layers of each burst type have to be summed up and converted to the .UD class.
To calculate the dBBMM for a specific or a subset of burst types, these have to be stated in the argument “burstType”. In this case the dBBMM is calculated only taking into account the burst segments “Night”:
From the resulting object, the UD per layer, or the total UD can be calculated as above.
Issues when calculating the dBBMM
Below some common issues that may arise while calculating the dBBMM and possible solutions to them.
Calculation takes very long
There can be several reasons why the calculation is taking very long:
- long trajectories of thousands or tens of thousands locations will take some time
- the pixel size for which the calculation is done will also be reflected on the calculation time. The smaller the pixel size, the longer the calculation time. Always make sure that the pixel size of the raster is larger than the location error.
- the time lag between locations is very small,
- because the data were collected at a very high resolution or
- because sometimes, there is a location that is just seconds or milliseconds apart from the next, even though the GPS schedule was set to minutes or hours.
By default the dBBMM will take the shortest time lag (in minutes) of the trajectory divided by 15 as the time interval taken for every integration step. If the shortest time lag is e.g. 15min, than the trajectory will be divided into steps of 1min. But if for some reason there is one time lag of 1sec (0.0167mins), the trajectory will be divided into steps of 0.001 mins, increasing the calculation time immensely.
Using the argument time.step the integration step time can be set manually. As a general rule of thumb one can take the intended GPS schedule in minutes divided by 15.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.200 1.875 2.117 9.422 6.883 60.700
For trajectories collected at a very high resolution, e.g. 1Hz, or trajectories containing high resolution bursts, the time.step argument can be set to 0.017 min (1sec), as there is no need to estimate where the animal was at a shorter time interval, because it is known where it was every 1sec.
In any case giving higher values to the time.step argument will speed up the calculation time.
Error: extent is not large enough or resulting UD is one “big blob”
This problem is mostly due to data that contains large chunks of missing data. During these large time gaps, the uncertainty of where the animal may have been is very large, or in other words the animal could have been anywhere.
Here is an example:
## Error in .local(object, raster, location.error = location.error, ext = ext, : Lower x grid not large enough, consider extending the raster in that direction or enlarging the ext argument
## making the extent bigger seems to solve the problem
dbb <- brownian.bridge.dyn(leroyWithGap_p, raster=100, location.error=20, ext=4)
## but than the UD is not very informative
ud <- getVolumeUD(dbb)
par(mfrow=c(1,2))
plot(ud, main="UD")
contour(ud, levels=c(0.5, 0.95), add=TRUE, lwd=c(0.5, 0.5), lty=c(2,1))
plot(ud, main="UD with locations")
points(leroyWithGap_p, col="red", cex=.5, pch=20)
contour(ud, levels=c(0.5, 0.95), add=TRUE, lwd=c(0.5, 0.5), lty=c(2,1))
And the solution of such a case:
The solution is to remove the variance of the segments corresponding to the large time gaps. For this first the variance is calculated with the brownian.motion.variance.dyn() function, then the segments corresponding to the large time gaps are set to FALSE, and finally the dBBMM is calculated.
## calculate the dynamic brownian motion variance of the gappy track
dbbv <- brownian.motion.variance.dyn(leroyWithGap_p, location.error=20, window.size=31, margin=11)
## the intended GPS fix rate of leroy was 15min, so we will ignore for example all segments that have a larger time lag than 5hours. The 'dBMvariance' object resulting from the function above, contains the slot '@interest' in which those segments marked as FALSE won't be included in the calculation of the dBBMM. Therefore we set all segments with time lag larger than 300mins to false
dbbv@interest[timeLag(leroyWithGap_p,"mins")>300] <- FALSE
## then we use the 'dBMvariance' object to calculate the dBBMM
dbb.corrected <- brownian.bridge.dyn(dbbv, raster=100, ext=.45,location.error=20)
## now the UD makes more sense
ud.corrected <- getVolumeUD(dbb.corrected)
par(mfrow=c(1,2))
plot(ud.corrected, main="UD")
contour(ud.corrected, levels=c(0.5, 0.95), add=TRUE, lwd=c(0.5, 0.5), lty=c(2,1))
plot(ud.corrected, main="UD with locations")
points(leroyWithGap_p, col="red", cex=.5, pch=20)
contour(ud.corrected, levels=c(0.5, 0.95), add=TRUE, lwd=c(0.5, 0.5), lty=c(2,1))
