Using the movegroup package

Install, help, background reading

To install on Linux systems it may be required to type, in terminal:

sudo apt install libgeos-dev
sudo apt install libproj-dev
sudo apt install libgdal-dev

then manually install rgeos and rgdal in R/RStudio. However be advised these packages are being depreciated in late-2023.

Also see each script’s Details section in the manual pages, as these frequently contain tips or common bugfixes.

Academic background and rationale

An animal’s home range can be defined as the area traversed by an animal within which it engages in normal activities e.g., foraging, mating (Burt, 1943). The most basic home range estimation method is the minimum convex polygon (MCP), which draws a polygon enclosing all relocations (Mohr, 1947). However, due to subjectivity of the home range definition, a new term was coined: the utilization distribution (UD; Winkle, 1975).

The UD is an improvement upon the MCP-derived home range estimation because this definition not only quantifies the size of a tracked animal’s home range, but also the intensity with which sub-areas within the home range are used (Worton, 1989). The kernel utilization distribution (KUD) is the first method that incorporated UDs for home range estimation. This method applies a bivariate kernel function over each location fix, after which resulting values are averaged. However, the KUD assumes no spatial and temporal autocorrelation in the data, which is of course unrealistic in the context of animal movement data that are inherently spatially and temporarily autocorrelated.

Relocations and UD
Relocations and UD

Horne et al. (2007) describe the use of movement models that incorporate Brownian motion (Brownian bridge moment model, BBMM), which offers a more sophisticated way to estimate space use. The traditional BBMM integrates the temporal component of tracking data by explicitly modeling movements between consecutive relocations. This is achieved by accounting for both the order of subsequent relocations as well as the travel time between them. The BBMM reconstructs the movement path by computing biased random walk iterations, creating a probability density distribution or “bridge” between two consecutive relocations where the animal could have been when it was not detected.


A minimal time difference among relocations is then suggestive of a straight-line movement, and therefore the variance in Brownian motion (the associated error of the movement path) would be small; conversely, a larger time difference is suggestive of a tortuous movement path, which would result in a larger Brownian motion variance of the bridge.


Kranstauber et al. (2012) later introduced the dynamic BBMM (dBBMM), which improves upon the traditional BBMM in its calculation of utilization distribution (UD) statistics. While the Brownian motion variance in the BBMM is static i.e., constant throughout the movement track, the dBBMM allows this variance to vary to allow for changes in in behaviour along the movement trajectory (e.g., foraging, travelling, resting, etc.). The result is improved performance in predicting animal locations. The dBBMM is also better equipped in dealing with irregular sampling of tracks, and would therefore be applied more appropriately to telemetry data for which you can anticipate gaps in detection (e.g., due to tags requiring to breach the ocean surface to transmit their location to a satellite (‘Smart Position Only Tag’ or SPOT), which depends on the animal’s behaviour, or tag-equipped aquatic animals leaving a fixed acoustic receiver array, or a tag’s line of sight with a satellite being obstructed by physical structures, etc.).

The dBBMM is calculated using the ‘move’ package. While this package is great for calculating UDs, the package can only calculate a model and output a UD for a single individual. Grouping individuals together to create a group-level/aggregated UD is not realistic, because the brownian.bridge.dyn() function requires a chronological movement path as input; grouping multiple individuals together would imply that individuals can teleport. This package that builds on the move package by being able to handle multiple individuals simultaneously, and aggregates individual UDs in a single group-level UD, offering significant advancements in the investigation of group-/population level space use estimation of telemetered animals. Of additional benefit is the ability to incorporate heterogeneous survey design e.g. unbalanced numbers of receivers across multiple arrays. All functions are designed to maximally automate the typical methodological pipeline, offloading the workload and technical skill required to (e.g.) scale and re-project multiple movement tracks to an optimal shared projection and extent, and plot output maps containing various disparate elements.

We strongly recommend that you download papers:

Kranstauber, B., Kays, R., LaPoint, S. D., Wikelski, M. and Safi, K. (2012) A dynamic Brownian bridge movement model to estimate utilization distributions for heterogeneous animal movement. Journal of Animal Ecology.

Kranstauber, B., M. Smolla & A. K. Scharf. (2019) Move: visualizing and analyzing animal track data. R package version 4.2.4 (at 2023-08-15).

van Zinnicq Bergmann, M. P. M., Guttridge, T. L., Smukall, M. J., Adams, V. M., Bond, M. E., Burke, P. J., Fuentes, M. M. P. B., Heinrich, D. D. U., Huveneers, C., Gruber, S. H., and Papastamatiou, Y. P. (2022) Using movement models and systematic conservation planning to inform marine protected area design for a multi-species predator community. Biological Conservation.

Strickland, B. A., Gastrich, K., Beauchamp, J. S., Mazzotti, F. J. and Heithaus, M. R. (2021) Effects of hydrology on the movements of a large-bodied predator in a managed freshwater marsh. Hydrobiologia

Also it’s imperative you read the R help files for each function before you use them. In RStudio: Packages tab, scroll to movegroup, click its name, the click the function to see its man (manual) page. Read the whole thing. Function man pages can also be accessed from the console by typing ?function.


Visualizing and Quantifying Space Use Data for Groups of Animals

Automates dynamic Brownian bridge movement model calculation for utilization distribution (UD) estimation for multiple individuals simultaneously, using functions in the ‘move’ package. The authors are indebted to the move package authors Bart Kraunstauber, Marco Smolla, and Anne K Scharf, and to Sarah Becker for seed code which inspired the development of the movegroup::movegroup function.

The movegroup function takes a dataframe of positions, datetimes, and IDs, e.g. from marine species, such as sharks or turtles being tracked by satellite or acoustic tags. These data are filtered to remove individuals with too few data to process, a group-level basemap raster is created upon which to calculate and later plot dBBMM utilisation distribution hotspots, then each individual has a movement track calculated using the move package, data gaps are assessed and the track split into multiple segments if gaps are too long, those segments have a dynamic Brownian Bridge Movement Model and variance calculated on them again using the move package, then the 50 and 95% utilisation distribution volume areas are calculated. Finally the outputs are saved on a per-unique-ID basis.

For this and all functions, see the function’s help for specifics on the function parameters, as well as errors and their origins.

Running the function is as simple as:

mysavedir <- "/your/directory/here/"
 data = TracksCleaned,
 ID = "Shark",
 Datetime = "Datetime",
 Lat = "Lat",
 Lon = "Lon",
 savedir = mysavedir)

However please appraise yourself of the meaning of the various parameters as they pertain to your data, most notably movement error distance (see moveLocErrorCalc below), and variously buffpct, rasterExtent, rasterResolution, centre for sizing your data and rasters, which has implications for later plotting. It might be that you need to adjust these elements after later seeing the resulting plots from plotraster, or if an overly large raster has caused a crash.


Scales Individual Utilization Distribution Rasters and Volume Area Estimates

Scales individual-level utilization distribution (UD) rasters from 0 to 1 to facilitate interpretation as relative intensity of utilization (as opposed to absolute), making comparisons across individuals and interpretations at the group level more straightforward. Subsequently, scaled individual-level rasters are aggregated to create a single group-level UD raster. See for issues, feedback, and development suggestions. There is an option to account for bias in acoustic receiver array spatial representation (see Details).

The process employed by the function is:

  1. Scale rasters: Individual-level UD rasters are scaled from 0 to 1 by dividing each raster by the maximum probability density value occurring within the raster set.

  2. Aggregate into a group-level raster: Scaled individual-level rasters are summed to create a single group-level UD raster.

  3. Re-scale to 0 to 1: The group-level raster is divided by its own maximum value.

  4. Weight raster (optional): The scaled group-level UD raster is divided by the specified weighting factor(s). Note that this is only useful if you want to account for an unbalanced listening station (e.g., acoustic receivers) array and have split up the study site and receivers in regions, and have run movegroup for each regional dataset separately. See van Zinnicq Bergmann et al. 2022 for example. If not applicable, leave as 1, the default.

  5. Standardize raster: Standardize the potentially weighted and scaled group-level UD raster so that its values sum to 1.

  6. Export as geographical (latlon) CRS file: Change crs to latlon for plotting and calculation purposes, save file, continue.

  7. Estimate 50 and 95pct contour volume areas: For each scaled individual-level UD raster, estimate 50 and 95pct contour volume areas, as well as their mean and standard deviation. Additionally, the 50 and 95pct volume area is estimated for the group-level UD raster.

  8. Export the projected-CRS group-level raster.

Having run the movegroup function above, scaleraster is run with:

scaleraster(path = mysavedir)

To weigh by number of positions per ID, fewer locations = lower weighting value = higher final values after dividing by weighting. This scales all IDs up to match the group max.

Weighting <- TracksCleaned |>
   dplyr::group_by(Shark) |>
   dplyr::summarise(N = n()) |> 
   dplyr::filter(N > 23) |> 
   dplyr::mutate(N = N / max(N, na.rm = TRUE)) |> 
scaleraster(path = mysavedir, weighting = Weighting)


Combines Region-Specific Group-Level UD Rasters into a Single Raster

Extends the spatial extent of each area-specific group-level raster to the spatial extent shared by all rasters. This will only be required if you have multiple individuals (e.g. different sharks) divided amongst a few discrete areas (e.g. around different islands) and the effort (e.g. receiver coverage) is different among islands. Not required for multiple individuals all within the same region or sampling regime.

To loop movegroup and scaleraster through tide subsets:

tide <- c("H", "M", "L")
for (i in tide) {
 dir.create(paste0(mysavedir, i))
    data = TracksCleaned[TracksCleaned$T.Ph == i, ],
    ID = "Shark",
    Datetime = "Datetime",
    Lat = "Lat",
    Lon = "Lon",
    savedir = paste0(mysavedir, i, "/"))
  scaleraster(path = paste0(mysavedir, i),
              crsloc = paste0(mysavedir, i))

alignraster(folderroots = paste0(mysavedir, tide),
            foldernames = tide,
            savefolder = paste0(mysavedir, "Aligned"))


Plots a Group-Level Utilization Distribution

Plots 50 and 95pct contours of a group-level utilization distribution raster on a spatial map background. Contains functionality to also visualize geographic locations of individual listening stations (e.g., acoustic receivers) as well as the entire surface UD.

This function plots the outputs of scaleraster, and individual movegroup rasters if desired.

For plottitle, you can use the term ‘home range’ when an animal can be detected wherever it goes i.e. using GPS, satellite or acoustic telemetry whereby it is known that acoustic receivers cover the entire home range of the study species. This term is problematic when applied to a passive acoustic telemetry setting where an array of non-overlapping receivers are used to assess local space use patterns i.e. the home range is bigger than the coverage by the acoustic array.

See the function’s help for specifics on the function parameters, as well as errors and their origins.

To get Google map basemaps:

(from here):

  1. Sign up with dev console
    1. You must enter credit card details, but won’t be charged if your daily API requests stay under the limit.
    2. Follow the link.
    3. Sign up for Google cloud account (it may auto populate your current gmail), click agree and continue.
    4. Click the navigation email in the top left corner and click on Billing.
    5. Create a billing account – they will NOT auto charge after trial ends.
    6. Enter information, click on ‘start my free trial’. They may offer a free credit for trying out their service. More pricing details here.
    7. Click “Select a Project” then “New project” in the top right corner.
    8. Enter Project Name, leave Location as is, click “Create”.
    9. You should now see your project name at the top, where the drop-down menu is.
  2. Enable Maps and Places API
    1. Click ‘Library’ on the left.
    2. In the search field type “Maps”.
    3. Scroll down, click “Maps Java Script API”.
    4. Click Enable.
    5. Click ‘Library’ again, search “Places”, click on “Places API”.
    6. Click Enable.
  3. Create Credentials for API Key
    1. Return to ‘APIs & Services’ page.
    2. Click on Credentials.
    3. At the top click ‘Create Credentials > API Key’.
    4. API key should pop up with option to copy it.
    5. You can restrict the key if you want by following steps 4 & 5 here.

Having run the movegroup and scaleraster function examples:

 x = paste0(mysavedir, "Scaled/All_Rasters_Scaled_Weighted_UDScaled.asc"),
 mapsource = "stamen",
 maptype = "terrain",
 savedir = paste0(mysavedir, "Plot"),
 xlatlon = paste0(mysavedir, "Scaled/All_Rasters_Scaled_Weighted_LatLon.asc"),
 locationpoints = TracksCleaned |> dplyr::rename(lat = "Lat", lon = "Lon"),
 pointsincontourssave = paste0(mysavedir, "Scaled/pointsincontours.csv"))

If you’ve setup your Google maps API, you can expect graphics such as: