Introduction

Example: White-bearded wildebeest (2010-2012)

Today, we will use part of a dataset collected during Jared Stabach’s field research in Kenya from 2010-2012. This dataset is freely available on Movebank and consist of GPS tracking data from 12 adult white-bearded wildebeest (Connochaetes taurinus). The goal of this exercise is to demonstrate how to work with animal movement data and combine it with various types of remote sensing data.

In this exercise, you will:

Import and clean an animal movement trajectory
Annotate animal movement data using existing environmental data from remote sensing, including land cover, elevation, and NDVI (Normalized Difference Vegetation Index)
Annotate animal movement with customized environmental information. In this practice, we will calculate the distance to the nearest urban area and various terrain indices as an example.

References

You can refer to these data in your research using the following citations:

Stabach JA, Hughey LF, Crego RD, Fleming CH, Hopcraft JGC, Leimgruber P, Morrison TA, Ogutu JO, Reid RS, Worden JS, Boone RB. 2022. Increasing anthropogenic disturbance restricts wildebeest movement across East African grazing systems. Frontiers in Ecology and Evolution. doi.org/10.3389/fevo.2022.846171

Stabach JA, Hughey LF, Reid RS, Worden JS, Leimgruber P, Boone RB. 2020. Data from: Comparison of movement strategies of three populations of white-bearded wildebeest. Movebank Data Repository. doi:10.5441/001/1.h0t27719

Setup

We have already covered data cleaning in previous sessions of our workshop. Therefore, we will start by loading a pre-processed white-beared wildebeest dataset. This will allow us to focus on the integration and analysis of remote sensing data.

Load Libraries

First, we need to load the required libearies and clear any existing data from R’s memory. This ensures that we start with a clean workspace and have all the necessary tools for our analysis.

# Clear all objects from the workspace to ensure a clean environment
rm(list=ls())

# Load required libraries
library(sf)         # Simple features for R
library(rgdal)      # Geospatial data abstraction library
library(sp)         # Spatial data classes
library(rgeos)      # Interface to geometry engine
library(dplyr)      # Data manipulation
library(raster)     # Raster data manipulation
library(viridis)    # Color scales
library(ggplot2)    # Data visualization
library(mapview)    # Interactive maps
library(terra)      # Raster and vector data manipulation
library(lubridate)  # Date and time manipulation

##########################################
#Version 2: only use sf and terra package

# Load required libraries
#library(sf)         # Simple features for R
#library(dplyr)      # Data manipulation
#library(terra)      # Raster and vector data manipulation
#library(ggplot2)    # Data visualization
#library(mapview)    # Interactive maps
#library(lubridate)  # Date and time manipulation

Animal Data

Next, we will load the animal movement data, filter out any rows with missing coordinates and set the coordinate reference system (CRS) to UTM zone 37S. This ensures that all spatial data is in the same projection, making it easier to integrate and analyze.

# Load the pre-processed animal movement data from an RDS file
WB <- readRDS("wildebeest_3hr_adehabitat.rds")

# Check the basic structure and class of the data to understand its format
class(WB)

## [1] "data.frame"

str(WB)

## 'data.frame':    52847 obs. of  17 variables:
##  $ x        : num  281442 282340 281930 282040 282187 ...
##  $ y        : num  9808391 9809132 9808660 9808734 9808773 ...
##  $ date     : POSIXct, format: "2010-10-19 03:00:00" "2010-10-19 06:00:00" ...
##  $ dx       : num  898.3 -410.7 109.9 147.7 -23.7 ...
##  $ dy       : num  740.4 -471.6 74.1 38.3 97.9 ...
##  $ dist     : num  1164 625 133 153 101 ...
##  $ dt       : num  10800 10800 10800 10800 10800 10800 10800 10800 10800 10800 ...
##  $ R2n      : num  0 1355089 310028 474624 700726 ...
##  $ abs.angle: num  0.689 -2.287 0.593 0.254 1.808 ...
##  $ rel.angle: num  NA -2.977 2.88 -0.339 1.555 ...
##  $ id       : Factor w/ 12 levels "Karbolo","Kikaya",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ burst    : chr  "Karbolo" "Karbolo" "Karbolo" "Karbolo" ...
##  $ sex      : Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...
##  $ DOP      : num  3.4 2.8 1.8 1.6 1.8 2.2 4.4 3.6 1.6 3.2 ...
##  $ fixType  : num  3 3 3 3 3 3 3 3 3 3 ...
##  $ temp     : num  17 17 25 30 28 25 22 18 15 17 ...
##  $ speed    : num  0.10779 0.0579 0.01227 0.01413 0.00933 ...

# Filter out rows with missing coordinates to ensure all data points are valid
WB <- WB %>%
  filter(!is.na(x) & !is.na(y))

#Convert the data frame to an sf object and set the CRS to UTM zone 37 S
WB_utm <- WB %>%
  st_as_sf(coords = c("x", "y"), crs = 32737)

Remote Sensing Data

To analyze the environment where the animal move, we will use remote sensing data. We can download this data from NASA AppEEARS platform, which requires creating an area sample to upload. Note that you need to register for an Earthdata account and log in to the website.

Create An Area Sample

Firstly, we need to create a minimum convex polygon (MCP) that encompasses the area where the animals move. This MCP will serve as the area sample to upload. The projection must be WGS 84.

# Convert the projection to WGS 84 (EPSG:4326)
WB_84 <- st_transform(WB_utm, crs = 4326)

# Create a minimum convex polygon (MCP) to cover the moving areas of all individuals
WB_84_mcp <- st_convex_hull(st_union(WB_84))

# Save the MCP as a shapefile
st_write(WB_84_mcp, "Boundary.shp")

## Writing layer `Boundary' to data source `Boundary.shp' using driver `ESRI Shapefile'
## Writing 1 features with 0 fields and geometry type Polygon.

# Load the MCP shapefile to confirm is was saved correctly
MCP=st_read("Boundary.shp")

## Reading layer `Boundary' from data source 
##   `/Users/Qianr/Documents/GitHub/Movement_workshop/Movement_workshop/animove-2023-remote-sensing-main/Boundary.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 36.78754 ymin: -1.865096 xmax: 37.11519 ymax: -1.338212
## Geodetic CRS:  WGS 84

Download Data

Next, we need to combine all the “Boundary” files into a single zip file and upload it to NASA’s AppEEARS platform.

Step 1

To start, navigate to the AppEEARS website, click “Extract” in the top left corner, and select “Area”. On the next page, choose “Start a new request”.

Fig 1.

Step 2

From the top to the bottom, set a name for your project, upload the zip file that contains Boundary.shp, Boundary.dbf, Boundary.prj, and Boundary.shx files (the shapefile documents we just created), and select the start date as 10-15-2010, and the end date as 12-04-2012.

Fig 2.

Step 3

Next, we will select the layers we need for our analysis:

Landcover Type
- Type “MCD12Q1.061” or “Combined MODIS Land Cover Type” - “LC_Type1”
Elevation Data
- Type “SRTMGL1_NC.003” or “SRTM Elevation (DEM)” - “SRTMGL1_DEM”
NDVI (Normalized Difference Vegetation Index)
- Type “MOD13Q1.061” or “Terra MODIS Vegetation Indices (NDVI & EVI)” - “_250M_16_days_NDVI”

Fig 3.

Fig 4.

Step 4

After selecting the environmental layers, we need to choose the format and projection for the output. For the format, choose “GeoTiff”, and for the projection, choose “Geographic”, which is WGS 84. Finnaly, click “Submit”.

Fig 6.

Step 5

Once submitted, you will receive a confirmation email indicating that the data is ready for download. After downloading, create three folders (Landcover, Elevation and NDVI) to store the data downloaded from NASA AppEEARS, as shown in the following picture.

Fig 7.

Land Cover Data

We used the Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) Version 6.1 data product. It utilized data from both Terra and Aqua satellites and provide global land cover types at yearly intervals (2001-2022). The land cover types are derived from the International Geosphere-Biosphere Programme (IGBP). And the spatial resolution is 500 meters.

In this section, we will visualize the land cover types and calculate the distance from each GPS location to the nearest urban areas.

Visualize Land Cover Types

We will load, stack, rename each layers of the land cover data according to the Annual International Geosphere-Biosphere Programme (IGBP) classification, and then visualize the land cover types for 2010.

Fig 8. Landcover classification

# Load land cover data for different years
landcover_2010 <- raster("./Landcover/MCD12Q1.061_LC_Type1_doy2010001_aid0001.tif")
landcover_2011 <- raster("./Landcover/MCD12Q1.061_LC_Type1_doy2011001_aid0001.tif")
landcover_2012 <- raster("./Landcover/MCD12Q1.061_LC_Type1_doy2012001_aid0001.tif")

# Stack the land cover data
landcover <- stack(landcover_2010, landcover_2011, landcover_2012)

##########################################
#Version 2: only use sf and terra package

# Load land cover data for different years
#landcover_2010 <- rast("./Landcover/MCD12Q1.061_LC_Type1_doy2010001_aid0001.tif")
#landcover_2011 <- rast("./Landcover/MCD12Q1.061_LC_Type1_doy2011001_aid0001.tif")
#landcover_2012 <- rast("./Landcover/MCD12Q1.061_LC_Type1_doy2012001_aid0001.tif")

# Stack the land cover data
#landcover <- c(landcover_2010, landcover_2011, landcover_2012)

# Define land cover types
landcover_levels <- c(
  "Evergreen needleleaf forests", "Evergreen broadleaf forests", "Deciduous needleleaf forests",
  "Deciduous broadleaf forests", "Mixed forests", "Closed shrublands", "Open shrublands",
  "Woody savannas", "Savannas", "Grasslands", "Permanent wetlands", "Croplands",
  "Urban and built-up lands", "Cropland/natural vegetation mosaics", "Snow and ice",
  "Barren", "Water bodies"
)

# Convert stacked land cover data to data frame
landcover_df <- as.data.frame(landcover, xy = TRUE, na.rm=TRUE)

# Rename land cover types for each year
landcover_df$landcover_2010 <- factor(landcover_df$MCD12Q1.061_LC_Type1_doy2010001_aid0001, levels = 1:17, labels = landcover_levels)
landcover_df$landcover_2011 <- factor(landcover_df$MCD12Q1.061_LC_Type1_doy2011001_aid0001, levels = 1:17, labels = landcover_levels)
landcover_df$landcover_2012 <- factor(landcover_df$MCD12Q1.061_LC_Type1_doy2012001_aid0001, levels = 1:17, labels = landcover_levels)

# Plot land cover types for 2010
ggplot(landcover_df) +
  geom_raster(aes(x = x, y = y, fill = landcover_2010)) +
  geom_sf(data = MCP, inherit.aes = FALSE, fill = NA) +
  scale_fill_manual(values = c(
    "Evergreen needleleaf forests" = "#05450a", "Evergreen broadleaf forests" = "#086a10",
    "Deciduous needleleaf forests" = "#54a708", "Deciduous broadleaf forests" = "#78d203",
    "Mixed forests" = "#009900", "Closed shrublands" = "#c6b044", "Open shrublands" = "#dcd159",
    "Woody savannas" = "#dade48", "Savannas" = "#fbff13", "Grasslands" = "#b6ff05",
    "Permanent wetlands" = "#27ff87", "Croplands" = "#c24f44", "Urban and built-up lands" = "#a5a5a5",
    "Cropland/natural vegetation mosaics" = "#ff6d4c", "Snow and ice" = "#69fff8",
    "Barren" = "#f9ffa4", "Water bodies" = "#1c0dff"
  )) +
  labs(title = "Land Cover 2010", fill = "Land Cover Type", x = "Longitude", y = "Latitude") +
  theme_minimal()+
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1))

From the above plot, we can see that there are only two land cover types in this area: Grassland and Urban area.

Calculating Distances to Urban Areas

We downloaded the land cover for three years (2010-2012). Since there is almost no change in the land cover during these three years, we will use the land cover data from 2010 to calculate the shortest distance from each GPS location to the nearest urban area.

# Since the animal data projection is UTM 37S (EPSG:32737), we use WB_utm.
WB_sp <- as(WB_utm, "Spatial")

# Identify urban areas (land cover type 13)
urban_mask <- calc(landcover_2010, fun = function(x) { x == 13 })
urban_cells <- Which(urban_mask, cells = TRUE)
keyXY <- xyFromCell(landcover_2010, urban_cells, spatial = TRUE)

# Convert to UTM 37S (EPSG:32737)
utm_crs <- CRS("+proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs")
keyXY_utm <- spTransform(keyXY, utm_crs)

# Calculate the shortest distance from each GPS location to the nearest urban pixel
distances <- gDistance(WB_sp, keyXY_utm, byid = TRUE)
min_distances <- apply(distances, 2, min)

# Add distances to the original data
WB_utm$min_distance_to_urban <- min_distances

# View the results
head(WB_utm)

## Simple feature collection with 6 features and 16 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 281442.1 ymin: 9808391 xmax: 282340.4 ymax: 9809132
## Projected CRS: WGS 84 / UTM zone 37S
##                  date         dx         dy      dist    dt       R2n
## 1 2010-10-19 03:00:00  898.30439  740.36329 1164.0827 10800       0.0
## 2 2010-10-19 06:00:00 -410.65289 -471.61806  625.3474 10800 1355088.6
## 3 2010-10-19 09:00:00  109.92265   74.08026  132.5552 10800  310028.0
## 4 2010-10-19 12:00:00  147.73146   38.28599  152.6119 10800  474624.2
## 5 2010-10-19 15:00:00  -23.69448   97.91944  100.7454 10800  700726.4
## 6 2010-10-19 18:00:00 -233.04236   43.83217  237.1287 10800  750193.2
##    abs.angle  rel.angle      id   burst sex DOP fixType temp       speed
## 1  0.6893111         NA Karbolo Karbolo   2 3.4       3   17 0.107785437
## 2 -2.2872039 -2.9765150 Karbolo Karbolo   2 2.8       3   17 0.057902538
## 3  0.5930148  2.8802187 Karbolo Karbolo   2 1.8       3   25 0.012273627
## 4  0.2535805 -0.3394343 Karbolo Karbolo   2 1.6       3   30 0.014130734
## 5  1.8082120  1.5546315 Karbolo Karbolo   2 1.8       3   28 0.009328282
## 6  2.9556780  1.1474660 Karbolo Karbolo   2 2.2       3   25 0.021956358
##                   geometry min_distance_to_urban
## 1 POINT (281442.1 9808391)              29456.10
## 2 POINT (282340.4 9809132)              29021.30
## 3 POINT (281929.7 9808660)              29342.04
## 4 POINT (282039.6 9808734)              29304.88
## 5 POINT (282187.4 9808773)              29314.03
## 6 POINT (282163.7 9808871)              29213.60

##########################################
#Version 2: only use sf and terra package

# Identify urban areas (land cover type 13)
#urban_mask <- classify(landcover_2010, cbind(13, 1), others = NA)

# Convert urban_mask to SpatRaster and then to SpatVector
#urban_points <- as.points(urban_mask, values=FALSE)

# Ensure the projection is UTM 37S (EPSG:32737)
#utm_crs <- "+proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs"

# Ensure urban_points CRS is the same as WB_sf
#urban_points <- project(urban_points, utm_crs)

# Convert urban_points to sf object
#urban_points_sf <- st_as_sf(urban_points)

# Ensure both sf objects have the same CRS
#urban_points_sf <- st_transform(urban_points_sf, crs = st_crs(WB_utm))

# Calculate the shortest distance from each GPS location to the nearest urban pixel
#distances <- st_distance(WB_utm, urban_points_sf, by_element = FALSE)

# Find the minimum distance for each GPS point
#min_distances <- apply(distances, 1, min)

# Add distances to the original data
#WB_utm$min_distance_to_urban <- min_distances

# View the results
#head(WB_utm)

Elevation Data

The elevation data we used here is SRTMGL1 v003, also known as the Shuttle Radar Topography Mission Global 1 arc-second dataset (version 003). It is a high-resolution digital elevation model (DEM) that provides global topographic data. The spatial resolution is 1 arc-second, which is approximately 30 meters. The unit is meter.

Load and Reproject Elevation Data

First, we will load the elevation data and reproject it to UTM 37S (EPSG: 32737).

# Load elevation data
elev <- raster("./Elevation/SRTMGL1_NC.003_SRTMGL1_DEM_doy2000042_aid0001.tif")

# Reproject to UTM 37S (EPSG:32737)
elev_utm <- projectRaster(elev, crs = utm_crs)

# Check the spatial resolution of the elevation data
res(elev_utm)

## [1] 30.9 30.7

##########################################
#Version 2: only use sf and terra package

# Load elevation data
#elev <- rast("./Elevation/SRTMGL1_NC.003_SRTMGL1_DEM_doy2000042_aid0001.tif")

# Reproject to UTM 37S (EPSG:32737)
#elev_utm <- project(elev, "EPSG:32737")

# Check the spatial resolution of the elevation data
#res(elev_utm)

The raster::cellStats function is useful for finding broad statistics about a layer such as sum, mean, min, max, and stanard deviation.

For Version 2 (only use sf and terra package), we use terra::global function to achieve the same goal.

# Calculate the mean and min value for the whole raster layer
cellStats(elev_utm, stat='mean',na.rm=TRUE)

## [1] 1633.514

cellStats(elev_utm, stat='min',na.rm=TRUE)

## [1] 1483.308

##########################################
#Version 2: only use sf and terra package

# Calculate the mean and min value for the whole raster layer
#global(elev_utm, c("mean", "min"), na.rm=TRUE)

The raster::aggregate function allows us to reduce the spatial resolution of the raster data.

For Version 2 (only use sf and terra package), we use terra::aggregate function.

# Load elevation raster
plot(elev_utm, main = "Original Elevation Data (res:30m)")

# Aggregate raster data to 300m resolution
elev_utm_300m <- aggregate(elev_utm, fact=10, fun=mean, expand=TRUE, na.rm=TRUE)
plot(elev_utm_300m, main = "Aggregated Elevation Data (res:300m)")

Calculate Terrain Indices

Topographical indices can be obtained by using raster::terrain function. They are according to Wilson et al. (2007). Belowed are explanations of the indices we calculate:

Slope: it measures the steepness or the degree of incline of a surface, representing the rate of change in elevation. The unit here is degree. Higher slope values indicate steeper terrain.

Aspect (orientation): the compass direction that a slope faces. The unit is degree, ranging from 0 to 360. 0 ° = North; 90° = East; 180° = South; 270° = West.

TPI (Topographic Position Index): it measures the relative position of a cell on the terrain by comparing its elevation to the mean elevation of its 8 surround cells. Positive TPI values indicate that the point is higher than its surroundings, such as ridges. Negative TPI values indicate that the point is lower than its surroundings, such as valleys.

TRI (Terrain Ruggedness Index): it measures the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells. Higher TRI values indicate more rugged terrain.

Roughness: the difference between the maximum and the minimum value of a cell and its 8 surrounding cells. Higher roughness values indicate more uneven terrain.

And flowdir returns the flow direction of water across the terrain.

Terrain Indices in 30-m resolution

# Calculate terrain indices
terrain_elev <- terrain(elev_utm, opt=c("slope", "aspect", "TPI", "TRI", 
                                         "roughness", "flowdir"), unit='degrees')
# Set up plotting layout for 2 rows and 3 columns
par(mfrow = c(2, 3))

# Plot each terrain index with individual titles using a loop
for (i in 1:nlayers(terrain_elev)) {
  plot(terrain_elev[[i]], main = names(terrain_elev)[i])
}

# Reset plotting layout to default
par(mfrow = c(1, 1))

##########################################
#Version 2: only use sf and terra package

# Calculate terrain indices
#terrain_elev <- terrain(elev_utm, v=c("slope", "aspect", "TPI", "TRI", 
#                                         "roughness", "flowdir"), unit='degrees')
# Set up plotting layout for 2 rows and 3 columns
#par(mfrow = c(2, 3))

# Plot each terrain index with individual titles using a loop
#for (i in 1:nlyr(terrain_elev)) {
#  plot(terrain_elev[[i]], main = names(terrain_elev)[i])
#}

# Reset plotting layout to default
#par(mfrow = c(1, 1))

Terrain Indices in 300-m resolution

# To observe the pattern more clearly, we can plot the elevation data in 30 m res as well.
terrain_elev_utm_300m <- terrain(elev_utm_300m, opt=c("slope", "aspect", "TPI", "TRI", 
                                         "roughness", "flowdir"), unit='degrees')
# Set up plotting layout for 2 rows and 3 columns
par(mfrow = c(2, 3))

# Plot each terrain index with individual titles using a loop
for (i in 1:nlayers(terrain_elev)) {
  plot(terrain_elev_utm_300m[[i]], main = names(terrain_elev_utm_300m)[i])
}

# Reset plotting layout to default
par(mfrow = c(1, 1))

##########################################
#Version 2: only use sf and terra package

# To observe the pattern more clearly, we can plot the elevation data in 30 m res as well.
#terrain_elev_utm_300m <- terrain(elev_utm_300m, v=c("slope", "aspect", "TPI", "TRI", 
#                                         "roughness", "flowdir"), unit='degrees')
# Set up plotting layout for 2 rows and 3 columns
#par(mfrow = c(2, 3))

# Plot each terrain index with individual titles using a loop
#for (i in 1:nlyr(terrain_elev)) {
#  plot(terrain_elev_utm_300m[[i]], main = names(terrain_elev_utm_300m)[i])
#}

# Reset plotting layout to default
#par(mfrow = c(1, 1))

Save Elevation Data and Specific Terrain Indices

# Create output directory
dir.create("output")

#Subset specific terrain index layers
slope <- subset(terrain_elev, "slope")
aspect <- subset(terrain_elev, "aspect")
TPI <- subset(terrain_elev, "tpi")
roughness <- subset(terrain_elev, "roughness")

# Save processed elevation data
# A raster layer
writeRaster(slope, filename="output/slope.tif", format="GTiff", overwrite=TRUE)
# A layer from a raster stack/brick
writeRaster(terrain_elev[[1]], filename="output/tri", format="GTiff", overwrite=TRUE)
# A whole stack/brick
writeRaster(terrain_elev, filename="output/terrain_elev", format="GTiff", overwrite=TRUE)
# A whole raster stack/brick by layer
writeRaster(terrain_elev, filename=paste("output/", names(terrain_elev), sep = ""), format="GTiff", bylayer=T, overwrite=TRUE)

Extract Elevation Data and Specific Terrain Indices

We will use terra::extract to extract the terrain indices at each individual’s GPS location and lapply to create a loop

For Version 2 (only use sf and terra package), we use terra::extract function.

# Put all the environmental layer into a list
raster_list <- list(
  elevation = elev_utm,
  TPI = TPI,
  slope = slope,
  aspect = aspect,
  roughness = roughness
)
# Extract based on individual's position
extracted_values <- lapply(raster_list, function(r) extract(r, WB_sp))

# Add the value back to the movement data
for (name in names(extracted_values)) {
  WB_utm[[name]] <- extracted_values[[name]]
}

# Check the modified animal movement data 
head(WB_utm)

## Simple feature collection with 6 features and 21 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 281442.1 ymin: 9808391 xmax: 282340.4 ymax: 9809132
## Projected CRS: WGS 84 / UTM zone 37S
##                  date         dx         dy      dist    dt       R2n
## 1 2010-10-19 03:00:00  898.30439  740.36329 1164.0827 10800       0.0
## 2 2010-10-19 06:00:00 -410.65289 -471.61806  625.3474 10800 1355088.6
## 3 2010-10-19 09:00:00  109.92265   74.08026  132.5552 10800  310028.0
## 4 2010-10-19 12:00:00  147.73146   38.28599  152.6119 10800  474624.2
## 5 2010-10-19 15:00:00  -23.69448   97.91944  100.7454 10800  700726.4
## 6 2010-10-19 18:00:00 -233.04236   43.83217  237.1287 10800  750193.2
##    abs.angle  rel.angle      id   burst sex DOP fixType temp       speed
## 1  0.6893111         NA Karbolo Karbolo   2 3.4       3   17 0.107785437
## 2 -2.2872039 -2.9765150 Karbolo Karbolo   2 2.8       3   17 0.057902538
## 3  0.5930148  2.8802187 Karbolo Karbolo   2 1.8       3   25 0.012273627
## 4  0.2535805 -0.3394343 Karbolo Karbolo   2 1.6       3   30 0.014130734
## 5  1.8082120  1.5546315 Karbolo Karbolo   2 1.8       3   28 0.009328282
## 6  2.9556780  1.1474660 Karbolo Karbolo   2 2.2       3   25 0.021956358
##                   geometry min_distance_to_urban elevation        TPI     slope
## 1 POINT (281442.1 9808391)              29456.10  1652.120  0.2408092 1.0618554
## 2 POINT (282340.4 9809132)              29021.30  1647.122  1.2992623 5.1687238
## 3 POINT (281929.7 9808660)              29342.04  1651.082 -0.7596832 0.8049379
## 4 POINT (282039.6 9808734)              29304.88  1654.829  2.5042965 2.5757412
## 5 POINT (282187.4 9808773)              29314.03  1645.105 -1.0588152 0.7801751
## 6 POINT (282163.7 9808871)              29213.60  1654.851  2.7578916 1.6830921
##      aspect roughness
## 1 237.04533  3.976277
## 2 346.10567  7.870842
## 3 277.38227  2.169356
## 4 171.71715  5.657715
## 5  44.39388  4.109192
## 6 233.28078  4.468331

##########################################
#Version 2: only use sf and terra package

# Put all the environmental layers into a list
#raster_list <- list(
#  elevation = elev_utm,
#  TPI = terrain_elev$TPI,
#  slope = terrain_elev$slope,
#  aspect = terrain_elev$aspect,
#  roughness = terrain_elev$roughness
#)

# Extract based on individual's position
#extracted_values <- lapply(raster_list, function(r) terra::extract(r, WB_utm, ID=FALSE))

# Combine extracted values into a data frame
#extracted_df <- do.call(cbind, extracted_values)

# Add column names to the extracted data frame
#colnames(extracted_df) <- names(raster_list)

# Ensure the extracted values are in the correct format
#head(extracted_df)

# Add the values back to the movement data
#WB_utm <- cbind(WB_utm, extracted_df)

# Check the modified animal movement data
#head(WB_utm)

NDVI Data

The Normalized Difference Vegetation Index (NDVI) is a measure of the amount and vigor of vegetation on the land surface and NDVI spatial composite images are developed to more easily distinguish green vegetation from bare soils. In general, NDVI values range from -1.0 to 1.0, with negative values indicating clouds and water, positive values near zero indicating bare soil, and higher positive values of NDVI ranging from sparse vegetation (0.1 - 0.5) to dense green vegetation (0.6 and above).

NDVI is also directly related to the:

“leaf area index” (LAI), which is often used in crop growth models
herbaceous or total green biomass (tons/ha) for given vegetation types
photosynthetic activity of the vegetation
percent ground cover

Indirectly, NDVI has been used to estimate the cumulative effective of rainfall on vegetation over a certain time period, rangeland carrying capacity, crop yields for different crop types, and the quality of the environment as habitat for various animals, pests and diseases.

The remote sensing product we used here is the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1) Version 6.1 data, which is generated every 16 days at 250 meter (m) spatial resolution.

Load NDVI Data

In this section, we will load NDVI data from the MODIS dataset, reproject it to UTM 37S, and prepare it for analysis.

# Load NDVI data
ndvi_files <- list.files(path = "./NDVI", pattern = "*_16_days_NDVI_.*\\.tif$", full.names = TRUE)

# Function to extract dates from NDVI filenames
extract_date_from_filename <- function(filename) {
  date_str <- sub(".*doy([0-9]{7}).*", "\\1", filename)
  year <- as.numeric(substr(date_str, 1, 4))
  day_of_year <- as.numeric(substr(date_str, 5, 7))
  date <- ymd(paste0(year, "-01-01")) + days(day_of_year - 1)
  as.POSIXct(date, tz="EAT")
}

# Extract dates from NDVI filenames
# Note: the time zone of MODIS data is UTC
ndvi_dates <- sapply(ndvi_files, extract_date_from_filename)
readable_ndvi_dates <- as.character(as.POSIXct(ndvi_dates, origin="1970-01-01", tz="UTC"))
ndvi_info <- data.frame(filename = basename(ndvi_files), date = readable_ndvi_dates)

# Load NDVI raster stack
ndvi_stack <- stack(ndvi_files)

# Reproject to UTM 37S (EPSG:32737)
ndvi_stack_utm <- projectRaster(ndvi_stack, crs = utm_crs)

##########################################
#Version 2: only use sf and terra package

# Load NDVI data
#ndvi_files <- list.files(path = "./NDVI", pattern = "*_16_days_NDVI_.*\\.tif$", full.names = TRUE)

# Function to extract dates from NDVI filenames
#extract_date_from_filename <- function(filename) {
#  date_str <- sub(".*doy([0-9]{7}).*", "\\1", filename)
#  year <- as.numeric(substr(date_str, 1, 4))
#  day_of_year <- as.numeric(substr(date_str, 5, 7))
#  date <- ymd(paste0(year, "-01-01")) + days(day_of_year - 1)
#  as.POSIXct(date, tz="EAT")
#}

# Extract dates from NDVI filenames
# Note: the time zone of MODIS data is UTC
#ndvi_dates <- sapply(ndvi_files, extract_date_from_filename)
#readable_ndvi_dates <- as.character(as.POSIXct(ndvi_dates, origin="1970-01-01", tz="UTC"))
#ndvi_info <- data.frame(filename = basename(ndvi_files), date = readable_ndvi_dates)

# Load NDVI raster stack
#ndvi_stack <- rast(ndvi_files)

# Reproject to UTM 37S (EPSG:32737)
#ndvi_stack_utm <- project(ndvi_stack, utm_crs)

In this part, we first define a function extract_date_from_filename to extract the date information from the filenames of the NDVI files. We then load the NDVI data and stack them together for easier handling. Finally, we reproject the stacked NDVI data to the UTM 37S coordinate reference system.

Match NDVI Data with Animal Movement Data

Next, we will match the NDVI data with the timestamps from the animal movement data to get the NDVI value at each animal location for the nearest avalible data.

# Ensure timestamp is POSIXct format with EAT timezone
WB_utm$date <- as.POSIXct(WB_utm$date, format = "%Y-%m-%d %H:%M:%S", tz="EAT")

# Function to find the nearest NDVI date
find_nearest_date <- function(date, readable_ndvi_dates) {
  diffs <- abs(difftime(readable_ndvi_dates, date, units = "days"))
  readable_ndvi_dates[which.min(diffs)]
}

# Find the nearest NDVI date for each timestamp
nearest_dates <- sapply(WB_utm$date, find_nearest_date, readable_ndvi_dates)

# Extract NDVI values for each timestamp
ndvi_values_list <- lapply(1:nlayers(ndvi_stack_utm), function(i) {
  extract(ndvi_stack_utm[[i]], WB_sp)
})

# Add NDVI values to the animal movement data
WB_utm$NDVI <- NA

for (i in 1:nrow(WB_utm)) {
  ndvi_layer_index <- which(readable_ndvi_dates == nearest_dates[i])
  if (length(ndvi_layer_index) > 0) {
    WB_utm$NDVI[i] <- ndvi_values_list[[ndvi_layer_index]][i]
  }
}

# Multiple the scale factor (0.0001) of the NDVI product MOD13Q1 V061
WB_utm$NDVI=WB_utm$NDVI*0.0001
# Check the data we made
head(WB_utm)

## Simple feature collection with 6 features and 22 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 281442.1 ymin: 9808391 xmax: 282340.4 ymax: 9809132
## Projected CRS: WGS 84 / UTM zone 37S
##                  date         dx         dy      dist    dt       R2n
## 1 2010-10-19 03:00:00  898.30439  740.36329 1164.0827 10800       0.0
## 2 2010-10-19 06:00:00 -410.65289 -471.61806  625.3474 10800 1355088.6
## 3 2010-10-19 09:00:00  109.92265   74.08026  132.5552 10800  310028.0
## 4 2010-10-19 12:00:00  147.73146   38.28599  152.6119 10800  474624.2
## 5 2010-10-19 15:00:00  -23.69448   97.91944  100.7454 10800  700726.4
## 6 2010-10-19 18:00:00 -233.04236   43.83217  237.1287 10800  750193.2
##    abs.angle  rel.angle      id   burst sex DOP fixType temp       speed
## 1  0.6893111         NA Karbolo Karbolo   2 3.4       3   17 0.107785437
## 2 -2.2872039 -2.9765150 Karbolo Karbolo   2 2.8       3   17 0.057902538
## 3  0.5930148  2.8802187 Karbolo Karbolo   2 1.8       3   25 0.012273627
## 4  0.2535805 -0.3394343 Karbolo Karbolo   2 1.6       3   30 0.014130734
## 5  1.8082120  1.5546315 Karbolo Karbolo   2 1.8       3   28 0.009328282
## 6  2.9556780  1.1474660 Karbolo Karbolo   2 2.2       3   25 0.021956358
##                   geometry min_distance_to_urban elevation        TPI     slope
## 1 POINT (281442.1 9808391)              29456.10  1652.120  0.2408092 1.0618554
## 2 POINT (282340.4 9809132)              29021.30  1647.122  1.2992623 5.1687238
## 3 POINT (281929.7 9808660)              29342.04  1651.082 -0.7596832 0.8049379
## 4 POINT (282039.6 9808734)              29304.88  1654.829  2.5042965 2.5757412
## 5 POINT (282187.4 9808773)              29314.03  1645.105 -1.0588152 0.7801751
## 6 POINT (282163.7 9808871)              29213.60  1654.851  2.7578916 1.6830921
##      aspect roughness      NDVI
## 1 237.04533  3.976277 0.2156442
## 2 346.10567  7.870842 0.2054445
## 3 277.38227  2.169356 0.2102625
## 4 171.71715  5.657715 0.2102625
## 5  44.39388  4.109192 0.2073611
## 6 233.28078  4.468331 0.2069555

In this step, we ensure the timestamps of the animal movement data are in the correct format and timezone. We then define a function find_nearest_date to find the nearest NDVI date for each timestamp in the animal movement data. We use this function to get the NDVI values for each location and add these values to our dataset WB_utm.

Save Data

Finally, we save the combined dataset with various environmental variables as RData file for future use. For example, you can examine animal resource use as what you will learn in the afternoon.

# Save the combined data
save(WB_utm, file = "WB_environmental_information.RData")

By following steps, you will have integrated NDVI data with the animal movement data, allowing you to analyze the relationship between vegetation and animal movement patterns.

Task after class

The time interval of NDVI data we used previously is 16-day. Is there any way to obtain finer temporal resolution data? There’re multiple ways to achieve this. You can combine Terra and Aqua MODIS data, use spatial-temporal interpolation, or even calculate daily NDVI data using daily updated MODIS bands.

Tips:

Option 1: Spatial-Termporal Interpolation

Option 2: Calculate daily NDVI

Calculate daily NDVI using the MCD43A4 v061 product. The MCD (MODIS Combined Data) reﬂectance is generated by combining data from Terrra and Aqua spacecrafts. The MCD43A4 product provides a Bidirectional Reﬂectance Distribution Function (BRDF) froma nadir view in the 1–7 MODIS bands, such as red, near infrared (NIR), blue, green, short wave infrared-1(SWIR-1), short wave infrared-2 (SWIR-2), and middle wave infrared (MWIR)). It is produced daily using 16 days of Terra and Aqua MODIS data at 500-meter resolution.

As we known, NDVI can be calculated from satellite imagery whereby the satellite’s spectrometer or radiometric sensor measures and stores reflectance values for both red and NIR bands on two separate channels or images. Kriegler, et al. (1969) were the first to propose NDVI and it is calculated by subtracting the red channel from the near-infrared (NIR) channel and dividing their difference by the sum of the two channels:

NDVI= (NIR - RED ) / (NIR + RED )

where, RED = the red portion of the electromagnetic spectrum (0.6-0.7 μm) NIR = the near infrared portion of the electromagnetic spectrum (0.75-1.5 μm).

Reference

Kriegler, F.J., Malila, W.A., Nalepka, R.F. and Richardson, W., 1969, Preprocessing transformations and their effect on multispectral recognition, in: Proceedings of the sixth International Symposium on Remote Sensing of Environment, University of Michigan, Ann Arbor, MI, pp. 97-131

Examples of Remote Sensing Data and Animal Movement Analyses

Qianru Liao (qianru@umd.edu), University of Maryland, College Park

2024-06-26