Dynamic Brownian Bridge Movement Models
Example script
Libraries
library(raster)## Loading required package: splibrary(move)## Loading required package: geosphere## Loading required package: rgdal## rgdal: version: 1.4-3, (SVN revision 828)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Program Files/R/R-3.5.3/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Program Files/R/R-3.5.3/library/rgdal/proj
## Linking to sp version: 1.3-1library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(rgdal)Load State Boundaries
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
States = readOGR(dsn = "./Bounds",
layer = "StUS_GCS12",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\humph173\Desktop\UD_Example\Bounds", layer: "StUS_GCS12"
## with 48 features
## It has 7 fieldsStates = spTransform(States, LL84)Project States
MProj = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" #Albers Equal Area (meters)
Statesplt = spTransform(States, MProj)Load Full Dataset (all birds)
MoveBirds.df = read.csv("Birds_041019.csv", stringsAsFactors=FALSE)Convert to SpatialPoints
MoveBirds = SpatialPointsDataFrame(MoveBirds.df[,c("Longitude","Latitude")], MoveBirds.df)
proj4string(MoveBirds) = proj4string(States)
#Project birds
MoveBirds = spTransform(MoveBirds, MProj)Create Raster Template for UDs
tmp.ras = raster(res=0.71) #degrees
extent(tmp.ras) = extent(States)
UD.r = rasterize(States,
tmp.ras,
field = 0,
background = NA)
UD.r1 = projectRaster(UD.r, crs=CRS(proj4string(MoveBirds)))
UD.r = extend(UD.r1, c(200,200))Filter Birds to Spring Migration
Also, updating projected coordinates.
MoveBirds.df = MoveBirds@data %>%
mutate(Long = MoveBirds@coords[,1],
Lat = MoveBirds@coords[,2],
Type = "gps",
Error = Error_Radius,
LocT = as.POSIXct(LocT,
format = "%Y-%m-%d %H:%M:%S")) %>%
filter(DOY >= 76 & DOY <= 165) %>%
select(Animal_ID, Long, Lat, LocT, Type, Error)Remove Duplicates
Software package canโt handle them.
MoveBirds.df$Duplicated = duplicated(MoveBirds.df[, c("Animal_ID", "LocT")])
MoveBirds.df = MoveBirds.df %>% filter(Duplicated == FALSE)
#Order by Bird ID
MoveBirds.df = arrange(MoveBirds.df, Animal_ID)Run Model in Loop
Bird.lvls = levels(factor(MoveBirds.df$Animal_ID)) #Individual birds
for(i in 1:length(Bird.lvls)){
MoveBirds.df2 = MoveBirds.df %>% filter(Animal_ID == Bird.lvls[i]) #Take one bird at a time
Fixes = length(MoveBirds.df2$Animal_ID) #Count number of records for bird
if(Fixes < 31){next} #If fewer records than window size, skip the bird
else{
Duration.m = as.numeric(max(MoveBirds.df2$LocT) - min(MoveBirds.df2$LocT)) #Calculate fliht duration (hrs)
#Create "move object"
BirdM = move(x=MoveBirds.df2$Long, #Longitude
y=MoveBirds.df2$Lat, #Latitude
time=MoveBirds.df2$LocT, #Date/Timestamp
proj=CRS(proj4string(MoveBirds)), #Projection
data=MoveBirds.df2, #data frame
animal=factor(MoveBirds.df2$Animal_ID)) #Bird IDs
#Run model
dbbmm1 = brownian.bridge.dyn(object=BirdM, #Move object
raster = UD.r, #raster template
location.error= MoveBirds.df2$Error,#Argos error radius
time.step=15, #time interval
window.size=31, #Median sampling @ 48min. 30*48 = 1 day, but must be odd
margin=11) #1/2 of window, shifted down to be odd #Kranstauber et al 2012
#Save results as raster
writeRaster(dbbmm1, paste("./EUD_rasters/", Bird.lvls[i], ".tif", sep=""), overwrite=TRUE)
#re-read the just saved raster
Temp.rast = raster(paste("./EUD_rasters/", Bird.lvls[i], ".tif", sep=""))
#Weight raster by flight duration
Temp.rast = Temp.rast*Duration.m
#Save weighted raster
writeRaster(Temp.rast, paste("./EUD_rasters/", Bird.lvls[i], ".tif", sep=""), overwrite=TRUE)
}
}## Computational size: 2.6e+10## Computational size: 3.9e+09## Computational size: 5.2e+09## Computational size: 5.3e+09## Computational size: 5.1e+09
## Computational size: 5.1e+09## Computational size: 3.9e+09## Computational size: 2.2e+09## Computational size: 2.8e+09## Computational size: 5.1e+09## Computational size: 2.5e+09## Computational size: 2.9e+10## Computational size: 1.7e+08## Computational size: 3.0e+09## Computational size: 2.9e+10## Computational size: 2.6e+09## Computational size: 5.2e+09## Computational size: 4.1e+09## Computational size: 3.7e+09## Computational size: 5.0e+09## Computational size: 2.1e+09## Computational size: 3.2e+09## Computational size: 2.7e+10## Computational size: 3.0e+09## Computational size: 4.8e+09## Computational size: 1.4e+09## Computational size: 4.4e+09## Computational size: 7.7e+08## Computational size: 4.0e+09## Computational size: 3.6e+09## Computational size: 6.3e+08Sum Individual Rasters
UDfiles = list.files(path="./EUD_rasters", #Read individual raster names
pattern="*.tif", full.names=T, recursive=T)
UD.stk = stack(UDfiles) #Create raster stack of individuals
UD.sum = sum(UD.stk) #Sum all individuals (weighted inprevious code chunk)
UD.sum = crop(UD.sum, extent(UD.r1)) #Crop out estra extension added at top
UD.sum = UD.sum + UD.r1 #Set areas outside of US to NA
UD.sum = projectRaster(UD.sum, crs=CRS(proj4string(States))) #Convert back to Lat-Long
Scale = function(x){(x-min(x, na.rm=T))/(max(x, na.rm=T)-min(x, na.rm=T))}
values(UD.sum) = Scale(values(UD.sum)) #Scale 0 to 1Classify
ValVect = values(UD.sum)
ValVect[ValVect == 0] = NA
Quant1 = quantile(ValVect, probs = 0.5, na.rm=T)
Quant2 = quantile(ValVect, probs = 0.75, na.rm=T)
Quant3 = quantile(ValVect, probs = 0.99, na.rm=T)
values(UD.sum) = ifelse(values(UD.sum) < Quant1, 1,
ifelse(values(UD.sum) >= Quant1 & values(UD.sum) < Quant2, 2,
ifelse(values(UD.sum) >= Quant2 & values(UD.sum) < Quant3, 3,
ifelse(values(UD.sum) >= Quant3, 4, NA))))Plot
plot(UD.sum)