a - Smithsonian National Zoo and Conservation Biology Institute, Conservation Ecology Center, 1500 Remount Rd, Front Royal, VA 22630, USA.
b - Working Land and Seascapes, Conservation Commons, Smithsonian Institution, Washington, DC 20013, USA.
c - Boyd Orr Centre for Population and Ecosystem Health, Institute of Biodiversity, Animal Health & Comparative Medicine (IBAHCM), University of Glasgow, G12 8QQ, UK.
The following tutorial describes the code workflow presented in the manuscript “Enhancing animal movement analyses - Spatiotemporal matching of animal positions with remotely sensed data using Google Earth Engine and R.”
The code workflow allows you to find the closest image from an image collection to the time at which each GPS location was acquired and extract the pixel value.
To use this code it is necessary to have a Google Earth Engine account and to install the rgee package.
0.1 Loading and preparing tracking data
The first step is to read a csv file with the telemetry data.
For this example we randomly selected 100 GPS fixes from the entire wildebeest dataset used in the manuscript, obtained from Movebank. The data is provided in the repository within the folder Data.
library(sf)
library(dplyr)
trackingdata <- read.csv("./Data/Data.csv", header = T) #Load your dataset
head(trackingdata)## timestamp location.long location.lat
## 1 2010-08-10 15:00:00 35.21622 -1.164488
## 2 2011-06-26 12:00:00 36.83339 -1.462866
## 3 2011-10-19 03:00:00 36.90053 -1.404862
## 4 2011-11-20 12:00:00 36.91039 -1.430901
## 5 2010-08-29 14:00:00 35.30668 -1.345507
## 6 2010-12-20 05:00:00 36.91189 -1.433012In the next step, we need to convert the dataframe into an sf object. We also need to set the date as a string with a ‘YYYY-MM-DDTHH:MM:SS’. We will use this data to convert it into milliseconds since midnight on January 1, 1970, a format used in Google Earth Engine to manage dates.
trackingdata$Date <- as.POSIXct(trackingdata$timestamp, format = "%Y-%m-%d %H:%M:%S", tz="UTC") #Modify as necessary
trackingdata$Date <- as.factor(trackingdata$Date)
trackingdata$Date <- sub(" ", "T", trackingdata$Date) #Put in a format that can be read by javascript
trackingdata$ID <- seq(1:nrow(trackingdata)) # Add ID to each point (optional)
datasf <- st_as_sf(trackingdata, coords = c('location.long','location.lat'), crs = 4326) #Transform the dataframe into sf object. Make sure the name of the columns for the coordinates match. CRS needs to be in longlat WGS84.
head(datasf)## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 35.21622 ymin: -1.462866 xmax: 36.91189 ymax: -1.164488
## Geodetic CRS: WGS 84
## timestamp Date ID geometry
## 1 2010-08-10 15:00:00 2010-08-10T15:00:00 1 POINT (35.21622 -1.164488)
## 2 2011-06-26 12:00:00 2011-06-26T12:00:00 2 POINT (36.83339 -1.462866)
## 3 2011-10-19 03:00:00 2011-10-19T03:00:00 3 POINT (36.90053 -1.404862)
## 4 2011-11-20 12:00:00 2011-11-20T12:00:00 4 POINT (36.91039 -1.430901)
## 5 2010-08-29 14:00:00 2010-08-29T14:00:00 5 POINT (35.30668 -1.345507)
## 6 2010-12-20 05:00:00 2010-12-20T05:00:00 6 POINT (36.91189 -1.433012)0.2 Initialize rgee
Next, we need to load and initialize rgee.
library(rgee)
ee_Initialize()
ee_check()0.3 Define Google Earth Engine functions
The next set of functions will be used to match images to the data and extract pixel values.
Note that you can edit the maximum temporal window allowed to find a match.
#Function to add property with time in milliseconds
add_date<-function(feature) {
date <- ee$Date(ee$String(feature$get("Date")))$millis()
feature$set(list(date_millis=date))
}
#Join Image and Points based on a maxDifference Filter within a temporal window
#Set temporal window in days for filter. This will depend on the remote sensing data used.
tempwin <- 16
#Set the filter
maxDiffFilter<-ee$Filter$maxDifference(
difference=tempwin*24*60*60*1000, #days * hr * min * sec * milliseconds
leftField= "date_millis", #Timestamp of the telemetry data
rightField="system:time_start" #Image date
)
# Define the join. We implement the saveBest function for the join, which finds the image that best matches the filter (i.e., the image closest in time to the particular GPS fix location).
saveBestJoin<-ee$Join$saveBest(
matchKey="bestImage",
measureKey="timeDiff"
)
#Function to add property with raster pixel value from the matched image
add_value<-function(feature){
#Get the image selected by the join
img1<-ee$Image(feature$get("bestImage"))$select(band)
#Extract geometry from the feature
point<-feature$geometry()
#Get pixel value for each point at the desired spatial resolution (argument scale)
pixel_value<-img1$sample(region=point, scale=250, tileScale = 16, dropNulls = F)
#Return the data containing pixel value and image date.
feature$setMulti(list(PixelVal = pixel_value$first()$get(band), DateTimeImage = img1$get('system:index')))
}
# Function to remove image property from features
removeProperty<- function(feature) {
#Get the properties of the data
properties = feature$propertyNames()
#Select all items except images
selectProperties = properties$filter(ee$Filter$neq("item", "bestImage"))
#Return selected features
feature$select(selectProperties)
}0.4 Load image collection
In this example, we are using NDVI from MODIS Terra Vegetation Indexes 16-Day Global 250m data set. However, you can use any other remote sensing product of interest and filter to the desired dates.
One of the main advantages offered by Google Earth Engine is the enormous amount of data available to be used. The constantly growing database consists on more than a petabyte archive of publicly available remotely sensed imagery and other related data sets. In the study cases of the manuscript we used MOD13Q1 and ERA5_LAND/HOURLY, but data available includes other products from MODIS, data from other satellites such as Landsat, National Oceanographic and Atmospheric Administration Advanced Very High Resolution Radiometer (NOAA AVHRR), Sentinel 1, 2, 3 and 5-P, Advanced Land Observing Satellite (ALOS), and other products such as sea surface temperature data, CHIRPS climate data, topography data, and land cover data. The entire list of datasets is available at this link.
Note that all image collections in Earth Engine have a code that you can extract from the link provided above and use to import into the workflow by modifying this example.
We will set the start and end days to filter the image collections. Temporal availability depends on each dataset.
We will also create an object with the name of the band we are interested in working with. The name of the band is also specific to each image collection.
start<-"2010-01-01"
end<-"2013-01-01"
imagecoll<-ee$ImageCollection('MODIS/006/MOD13Q1')$filterDate(start,end)
band <- "NDVI" #Name of the band to use. You can change to EVI for instance when using MOD13Q1.0.5 Extract pixel value
A key function in this process is the ee_as_sf which converts the Google Earth Engine table in a sf object. This function provides three different options to convert the table (feature collection) into a sf object:
- getInfo: which is fast and direct but has a limit of 5000 features
- drive: which exports data through your Google Drive account
- gsc: which exports data through your Google Cloud Storage account
You can find more information about this function in the help: ?ee_as_sf
We use here the getInfo option given it is direct and simple. However, this option has a limit of 5000 features to convert. For that reason, we are going to run a loop, processing 1000 features (points) per time to avoid errors. If memory limit errors are display, then you can reduce the number of points to extract each time by changing the each argument on the rep function.
In this example we only have 100 points so the loop will only run once, but for larger datasets the loop may run multiple times.
datasf$uniq <- rep(1:1000, each=1000)[1:nrow(datasf)] #This is for up to 1 million points. To increase the max number of points, increase the value for max repetitions. To change the number of points to run per time, change the value in the argument each (up to 5000).
start_time <- Sys.time()
dataoutput <- data.frame()
for(x in unique(datasf$uniq)){
data1 <- datasf %>% filter(uniq == x)
# Send sf to GEE
data <- sf_as_ee(data1)
# Transform day into milliseconds
data<-data$map(add_date)
# Apply the join
Data_match<-saveBestJoin$apply(data, imagecoll, maxDiffFilter)
# Add pixel value to the data
DataFinal<-Data_match$map(add_value)
# Remove image property from the data
DataFinal<-DataFinal$map(removeProperty)
# Move GEE object into R
temp<- ee_as_sf(DataFinal, via = 'getInfo')
# Append
dataoutput <- rbind(dataoutput, temp)
}## Registered S3 method overwritten by 'geojsonsf':
## method from
## print.geojson geojsonend_time <- Sys.time()The time needed to run 100 points was:
end_time - start_time## Time difference of 5.600979 secsThe new sf data frame with the pixel values in now stored as the dataoutput object. You can use this for further analysis.
names(dataoutput)[4] <- band
dataoutput## Simple feature collection with 100 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 34.99598 ymin: -2.693471 xmax: 37.5994 ymax: -1.148282
## Geodetic CRS: WGS 84
## First 10 features:
## Date DateTimeImage ID NDVI date_millis timestamp
## 1 2010-08-10 15:00:00 2010_08_13 1 3011 1552145792 2010-08-10 15:00:00
## 2 2011-06-26 12:00:00 2011_06_26 2 3632 -875425280 2011-06-26 12:00:00
## 3 2011-10-19 03:00:00 2011_10_16 3 2957 438240128 2011-10-19 03:00:00
## 4 2011-11-20 12:00:00 2011_11_17 4 4952 -1059527168 2011-11-20 12:00:00
## 5 2010-08-29 14:00:00 2010_08_29 5 2699 -1104821504 2010-08-29 14:00:00
## 6 2010-12-20 05:00:00 2010_12_19 6 3518 36043904 2010-12-20 05:00:00
## 7 2010-07-16 18:01:00 2010_07_12 7 3227 -596994208 2010-07-16 18:01:00
## 8 2011-12-09 03:01:00 2011_12_03 8 4460 549732832 2011-12-09 03:01:00
## 9 2012-01-01 07:00:00 2012_01_01 9 2624 -1743694464 2012-01-01 07:00:00
## 10 2012-04-16 04:00:00 2012_04_22 10 5840 -1186029056 2012-04-16 04:00:00
## uniq geometry
## 1 1 POINT (35.21622 -1.164488)
## 2 1 POINT (36.83339 -1.462866)
## 3 1 POINT (36.90053 -1.404862)
## 4 1 POINT (36.91039 -1.430901)
## 5 1 POINT (35.30668 -1.345507)
## 6 1 POINT (36.91189 -1.433012)
## 7 1 POINT (35.24742 -1.324172)
## 8 1 POINT (37.565 -2.446763)
## 9 1 POINT (37.03885 -2.620116)
## 10 1 POINT (35.52517 -1.233949)0.6 Visualize locations
Pres <- dataoutputlibrary(tmap)
tmap_mode('view')
tm_shape(Pres) + tm_dots(col = 'blue', title = "Presence")0.7 Session
devtools::session_info()## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.4 (2021-02-15)
## os macOS Big Sur 10.16
## system x86_64, darwin17.0
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz Europe/London
## date 2021-10-25
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date lib source
## abind 1.4-5 2016-07-21 [2] CRAN (R 4.0.2)
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.2)
## base64enc 0.1-3 2015-07-28 [2] CRAN (R 4.0.2)
## bslib 0.2.5.1 2021-05-18 [2] CRAN (R 4.0.2)
## cachem 1.0.5 2021-05-15 [2] CRAN (R 4.0.2)
## callr 3.7.0 2021-04-20 [2] CRAN (R 4.0.2)
## class 7.3-19 2021-05-03 [2] CRAN (R 4.0.2)
## classInt 0.4-3 2020-04-07 [2] CRAN (R 4.0.2)
## cli 2.5.0 2021-04-26 [1] CRAN (R 4.0.4)
## codetools 0.2-18 2020-11-04 [2] CRAN (R 4.0.4)
## crayon 1.4.1 2021-02-08 [1] CRAN (R 4.0.2)
## crosstalk 1.1.1 2021-01-12 [2] CRAN (R 4.0.2)
## crul 1.1.0 2021-02-15 [2] CRAN (R 4.0.2)
## curl 4.3.1 2021-04-30 [2] CRAN (R 4.0.2)
## DBI 1.1.1 2021-01-15 [2] CRAN (R 4.0.2)
## desc 1.3.0 2021-03-05 [2] CRAN (R 4.0.2)
## devtools 2.4.2 2021-06-07 [2] CRAN (R 4.0.2)
## dichromat 2.0-0 2013-01-24 [2] CRAN (R 4.0.2)
## digest 0.6.27 2020-10-24 [2] CRAN (R 4.0.2)
## dplyr * 1.0.6 2021-05-05 [1] CRAN (R 4.0.2)
## e1071 1.7-7 2021-05-23 [1] CRAN (R 4.0.2)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.0.2)
## evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.1)
## fansi 0.5.0 2021-05-25 [2] CRAN (R 4.0.2)
## fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.0.2)
## foreign 0.8-81 2020-12-22 [2] CRAN (R 4.0.4)
## fs 1.5.0 2020-07-31 [2] CRAN (R 4.0.2)
## generics 0.1.0 2020-10-31 [2] CRAN (R 4.0.2)
## geojson 0.3.4 2020-06-23 [2] CRAN (R 4.0.2)
## geojsonio 0.9.4 2021-01-13 [2] CRAN (R 4.0.2)
## geojsonsf 2.0.1 2020-10-02 [2] CRAN (R 4.0.2)
## glue 1.4.2 2020-08-27 [2] CRAN (R 4.0.2)
## htmltools 0.5.1.1 2021-01-22 [2] CRAN (R 4.0.2)
## htmlwidgets 1.5.3 2020-12-10 [2] CRAN (R 4.0.2)
## httpcode 0.3.0 2020-04-10 [2] CRAN (R 4.0.2)
## jqr 1.2.1 2021-05-06 [2] CRAN (R 4.0.2)
## jquerylib 0.1.4 2021-04-26 [2] CRAN (R 4.0.4)
## jsonlite 1.7.2 2020-12-09 [2] CRAN (R 4.0.2)
## KernSmooth 2.23-20 2021-05-03 [2] CRAN (R 4.0.2)
## knitr 1.33 2021-04-24 [2] CRAN (R 4.0.2)
## lattice 0.20-44 2021-05-02 [2] CRAN (R 4.0.2)
## lazyeval 0.2.2 2019-03-15 [2] CRAN (R 4.0.2)
## leafem 0.1.6 2021-05-24 [2] CRAN (R 4.0.2)
## leaflet 2.0.4.1 2021-01-07 [2] CRAN (R 4.0.2)
## leaflet.providers 1.9.0 2019-11-09 [2] CRAN (R 4.0.2)
## leafsync 0.1.0 2019-03-05 [2] CRAN (R 4.0.2)
## lifecycle 1.0.0 2021-02-15 [1] CRAN (R 4.0.2)
## lwgeom 0.2-6 2021-04-02 [2] CRAN (R 4.0.2)
## magrittr 2.0.1 2020-11-17 [2] CRAN (R 4.0.2)
## maptools 1.1-1 2021-03-15 [2] CRAN (R 4.0.2)
## Matrix 1.3-4 2021-06-01 [2] CRAN (R 4.0.2)
## memoise 2.0.0 2021-01-26 [2] CRAN (R 4.0.2)
## pillar 1.6.1 2021-05-16 [1] CRAN (R 4.0.2)
## pkgbuild 1.2.0 2020-12-15 [2] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.2)
## pkgload 1.2.1 2021-04-06 [2] CRAN (R 4.0.2)
## png 0.1-7 2013-12-03 [2] CRAN (R 4.0.2)
## prettydoc 0.4.1 2021-01-10 [1] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.0.2)
## processx 3.5.2 2021-04-30 [1] CRAN (R 4.0.2)
## proxy 0.4-26 2021-06-07 [1] CRAN (R 4.0.2)
## ps 1.6.0 2021-02-28 [1] CRAN (R 4.0.2)
## purrr 0.3.4 2020-04-17 [2] CRAN (R 4.0.2)
## R6 2.5.0 2020-10-28 [2] CRAN (R 4.0.2)
## raster 3.4-10 2021-05-03 [2] CRAN (R 4.0.2)
## RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 4.0.2)
## Rcpp 1.0.6 2021-01-15 [2] CRAN (R 4.0.2)
## remotes 2.4.0 2021-06-02 [2] CRAN (R 4.0.2)
## reticulate 1.20 2021-05-03 [1] CRAN (R 4.0.2)
## rgee * 1.0.9 2021-04-24 [1] CRAN (R 4.0.2)
## rgeos 0.5-5 2020-09-07 [2] CRAN (R 4.0.2)
## rlang 0.4.11 2021-04-30 [2] CRAN (R 4.0.2)
## rmarkdown 2.8 2021-05-07 [2] CRAN (R 4.0.2)
## rprojroot 2.0.2 2020-11-15 [2] CRAN (R 4.0.2)
## rstudioapi 0.13 2020-11-12 [2] CRAN (R 4.0.2)
## s2 1.0.5 2021-06-01 [1] CRAN (R 4.0.2)
## sass 0.4.0 2021-05-12 [2] CRAN (R 4.0.2)
## sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.2)
## sf * 1.0-0 2021-06-09 [1] CRAN (R 4.0.2)
## sp 1.4-5 2021-01-10 [2] CRAN (R 4.0.2)
## stars 0.5-3 2021-06-08 [2] CRAN (R 4.0.2)
## stringi 1.6.2 2021-05-17 [2] CRAN (R 4.0.2)
## stringr 1.4.0 2019-02-10 [2] CRAN (R 4.0.2)
## testthat 3.0.2 2021-02-14 [2] CRAN (R 4.0.2)
## tibble 3.1.2 2021-05-16 [1] CRAN (R 4.0.2)
## tidyselect 1.1.1 2021-04-30 [2] CRAN (R 4.0.2)
## tmap * 3.3-1 2021-03-15 [1] CRAN (R 4.0.2)
## tmaptools 3.1-1 2021-01-19 [2] CRAN (R 4.0.2)
## units 0.7-2 2021-06-08 [1] CRAN (R 4.0.2)
## usethis 2.0.1 2021-02-10 [2] CRAN (R 4.0.2)
## utf8 1.2.1 2021-03-12 [1] CRAN (R 4.0.2)
## V8 3.4.2 2021-05-01 [2] CRAN (R 4.0.2)
## vctrs 0.3.8 2021-04-29 [1] CRAN (R 4.0.2)
## viridisLite 0.4.0 2021-04-13 [1] CRAN (R 4.0.2)
## withr 2.4.2 2021-04-18 [1] CRAN (R 4.0.2)
## wk 0.4.1 2021-03-16 [2] CRAN (R 4.0.2)
## xfun 0.23 2021-05-15 [1] CRAN (R 4.0.2)
## XML 3.99-0.6 2021-03-16 [2] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.2)
##
## [1] /Users/ramirocrego/Library/R/4.0/library
## [2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library