In SMBC-NZP/MigConnectivity: Estimate Migratory Connectivity for Migratory Animals

knitr::opts_chunk$set(collapse = TRUE, 
                      comment = "#>", 
                      tidy = FALSE,
                      message = FALSE,
                      warning = FALSE)

Load required packages

library(terra)
library(MigConnectivity)
library(sf)

The following is a simulation that tests the how the spatial arrangement of target sites influences MC from stable-hydrogen isotopes. The following simulation is run using data generated within the code but we use the Ovenbird as an example species.

Ovenbird distribution

Read in the Ovenbird distribution and create a species distribution map from the abundance data.

# read in raster layer
download.file(
"https://raw.githubusercontent.com/SMBC-NZP/MigConnectivity/master/data-raw/Spatial_Layers/bbsoven.txt",
              destfile = "bbsoven.txt")

OVENabund <- terra::rast("bbsoven.txt")

OVENdist <- OVENabund
OVENdist[OVENdist>0]<-1
OVENdist[OVENdist==0]<-NA

OVEN_single_poly <- terra::as.polygons(OVENdist)
terra::crs(OVEN_single_poly) <- MigConnectivity::projections$WGS84
terra::crs(OVENabund) <- MigConnectivity::projections$WGS84

To complete the simulation we need a template to ensure the raster resolution is the same as the assignment raster. To do this, we use the isotope data as our template. We grab the isotope base-map using the getIsoMap function.

### get isotope raster for template
rasterTemplate <- getIsoMap()

Crop template to distribution and mask with distribution

terra::crs(rasterTemplate) <- MigConnectivity::projections$WGS84
rasterTemplate <- terra::crop(terra::mask(rasterTemplate, 
                                          OVEN_single_poly),
                              OVEN_single_poly)

# rasterize the distribution for relative abundance so that raster
# dimensions and resolution match the isotope layer
relativeAbund <- terra::project(OVENabund,rasterTemplate)
relativeAbund  <- relativeAbund / 
  unlist(c(terra::global(relativeAbund, fun = "sum", na.rm = T)))

Generate target sites

The simulation is focused on the effect of target sites on MC when using stable-hydrogen isotopes. The following code generates various target site layers used in the simulation.

# generate target sites
targetRanges <- vector('list',5)
# 3' latitude
targetRanges[[1]] <- terra::as.polygons(terra::rast(xmin = -180,
                                         xmax = -40,
                                         ymin = 25,
                                         ymax = 85,
                                         resolution = c(140,3)))

# 5' latitude
targetRanges[[2]] <- terra::as.polygons(terra::rast(xmin = -180,
                                             xmax = -40,
                                             ymin = 25,
                                             ymax = 85,
                                             resolution = c(140,5)))

# 10' latitude
targetRanges[[3]] <- terra::as.polygons(terra::rast(xmin = -180,
                                         xmax = -40,
                                         ymin = 25,
                                         ymax = 85,
                                         resolution = c(140,10)))

# 12 isotope units
featherIso <- (0.95*getIsoMap(period = "GrowingSeason")+(-17.57))
iso <- terra::crop(featherIso, terra::ext(c(-180,-40,25,85)))
isocut <- terra::classify(iso, 
                          rcl = seq(terra::global(iso, fun = "min", 
                                                  na.rm = TRUE)$min,
                                    terra::global(iso, fun = "max", 
                                                  na.rm = TRUE)$max,12))
targetRanges[[4]] <- terra::as.polygons(isocut)


# 12*2 isotope units
isocut <- terra::classify(iso, 
                          rcl = seq(terra::global(iso, fun = "min", 
                                                  na.rm = TRUE)$min, 
                                    terra::global(iso,fun = "max",
                                                  na.rm = TRUE)$max, 24))
targetRanges[[5]] <- terra::as.polygons(isocut)


TR <- lapply(targetRanges, sf::st_as_sf)
vTR <- lapply(TR, sf::st_make_valid)
vTR <- lapply(vTR, FUN = function(x){sf::st_set_crs(x, 4326)})

sf_use_s2(FALSE)
OVEN_single_poly <- sf::st_as_sf(OVEN_single_poly) %>% sf::st_make_valid()


# Keep only the targetSites that intersect with the OVEN polygon
targetRanges <- lapply(vTR, FUN = function(x){sf::st_intersection(x,OVEN_single_poly)})

targetRanges <- lapply(targetRanges, FUN = function(x){y <- sf::st_as_sf(x)
                                                       z <- sf::st_transform(y,4326)
                                                       return(z)})

for(i in 1:5){
  targetRanges[[i]]$Target <- 1:nrow(targetRanges[[i]])
}
targetRanges <- lapply(targetRanges, sf::st_make_valid)
sf_use_s2(TRUE)

Generate simulated data

#Generate random breeding locations using the 10' target sites
Site1 <- st_as_sf(sf::st_sample(targetRanges[[3]][1:2,], size =100, type = "random"))
Site2 <- st_as_sf(sf::st_sample(targetRanges[[3]][2:3,], size =100, type = "random"))
Site3 <- st_as_sf(sf::st_sample(targetRanges[[3]][2:4,], size =100, type = "random"))

# Capture coordinates
capCoords <- array(NA,c(3,2))
capCoords[1,] <- cbind(-98.17,28.76)
capCoords[2,] <- cbind(-93.70,29.77)
capCoords[3,] <- cbind(-85.000,29.836)

featherIso <- (0.95*getIsoMap(period = "GrowingSeason")+(-17.57))

# Extract simulated data
iso_dat <- data.frame(Site = rep(1:3, each = 100),
                      xcoords = c(sf::st_coordinates(Site1)[,1], 
                                  sf::st_coordinates(Site2)[,1], 
                                  sf::st_coordinates(Site3)[,1]),
                      ycoords = c(sf::st_coordinates(Site1)[,2], 
                                  sf::st_coordinates(Site2)[,2],
                                  sf::st_coordinates(Site3)[,2]),
                      targetSite = unlist(unclass(sf::st_intersects(rbind(Site1,
                                                                          Site2,
                                                                          Site3),
                                        targetRanges[[3]]))),
                      featherIso = terra::extract(featherIso,
                                                  terra::vect(rbind(Site1,Site2,
                                                                    Site3))))

iso_dat <- iso_dat[complete.cases(iso_dat),]

# generate transition data from simulation
sim_psi <- table(iso_dat$Site, iso_dat$targetSite)

for(i in 1:nrow(sim_psi)){
  sim_psi[i,]<-sim_psi[i,]/sum(sim_psi[i,])
}

states <- rnaturalearth::ne_states(country = "United States of America", 
                                   returnclass = "sf")
originSites <- states[(states$woe_name %in% c("Texas","Louisiana","Florida")),]
originSites <- sf::st_transform(originSites, 4326)

originDist <- distFromPos(st_coordinates(sf::st_centroid(originSites,
                                                            byid = TRUE)))
targetDistMC <- distFromPos(sf::st_coordinates(sf::st_centroid(targetRanges[[3]], 
                                                              byid = TRUE)))

Start of simulations

warning this takes a long time to run.

originRelAbund <- rep(1/3, 3)
nTargetSetups <- 5
nSims <- 2            #SET LOW FOR EXAMPLE
# nSims <- 200        
nOriginSites = 3

targetPoints0 <- matrix(NA, 300, 2, dimnames = list(NULL, c("x", "y")))


a <- Sys.time()
output.sims <- vector("list", nTargetSetups)
for(target in 1:nTargetSetups){
  sim.output <- data.frame(targetSetup = rep(NA,nSims),
                           sim = rep(NA,nSims),
                           MC.generated = rep(NA,nSims),
                           MC.realized = rep(NA,nSims),
                           MC.est = rep(NA,nSims),
                           MC.low = rep(NA,nSims),
                           MC.high = rep(NA,nSims),
                           rM.realized = rep(NA,nSims),
                           rM.est = rep(NA,nSims),
                           rM.low = rep(NA,nSims),
                           rM.high = rep(NA,nSims))
  set.seed(9001)
  targetSites <- targetRanges[[target]]
  sf_use_s2(FALSE)
  targetSites <- sf::st_make_valid(targetSites)
  targetDist <- distFromPos(sf::st_coordinates(sf::st_centroid(targetSites, 
                                                               byid = TRUE)))

  # Extract simulated data
  iso_dat <- data.frame(Site = rep(1:3,each = 100),
                        xcoords = c(sf::st_coordinates(Site1)[,1],
                                    sf::st_coordinates(Site2)[,1],
                                    sf::st_coordinates(Site3)[,1]),
                        ycoords = c(sf::st_coordinates(Site1)[,2],
                                    sf::st_coordinates(Site2)[,2],
                                    sf::st_coordinates(Site3)[,2]),
                        targetSite = unlist(unclass(sf::st_intersects(rbind(Site1,
                                                                            Site2,
                                                                            Site3),
                                        targetSites))),
                        featherIso = terra::extract(featherIso,
                                                    terra::vect(rbind(Site1,
                                                                      Site2,
                                                                      Site3))))

  iso_dat <- iso_dat[complete.cases(iso_dat),]

  # generate transition data from simulation
  sim_psi <- prop.table(table(iso_dat$Site,factor(iso_dat$targetSite,
                                                  1:nrow(targetSites))), 1)

  sf_use_s2(TRUE)
  MC.generated <- calcMC(originDist = originDist,
                          targetDist = targetDist,
                          originRelAbund = originRelAbund,
                          psi = sim_psi)

  for (sim in 1:nSims) {
    cat("Simulation Run", sim, "of", nSims, "for target",target,"at", date(), "\n")
    sim_move <- simMove(rep(100, nOriginSites), originDist, targetDist, sim_psi, 1, 1)
    originAssignment <- sim_move$animalLoc[,1,1,1]
    targetAssignment <- sim_move$animalLoc[,2,1,1]
    sf_use_s2(FALSE)
    for (i in 1:300) {
      targetPoint1 <- sf::st_sample(targetSites[targetAssignment[i],], 
                               size = 1, type = "random", iter = 25)
      targetPoints0[i,] <- sf::st_coordinates(targetPoint1)
    }
    targetPoints <- sf::st_as_sf(as.data.frame(targetPoints0), crs = 4326,
                                 coords = c("x", "y"))

    # Extract simulated data
    iso_dat <- data.frame(Site = originAssignment,
                          xcoords = targetPoints0[,1],
                          ycoords = targetPoints0[,2],
                          targetSite = unlist(unclass(sf::st_intersects(targetPoints,
                                                                        targetSites))),
                          featherIso = terra::extract(featherIso,
                                                      terra::vect(targetPoints)))


    iso_dat <- iso_dat[complete.cases(iso_dat),]

    # generate transition data from simulation
    sim_psi0 <- table(iso_dat$Site,
                      factor(iso_dat$targetSite,
                             min(targetSites$Target):max(targetSites$Target)))
    sim_psi_realized <- prop.table(sim_psi0, 1)

    # get points ready for analysis
    nSite1 <- table(iso_dat$Site)[1]
    nSite2 <- table(iso_dat$Site)[2]
    nSite3 <- table(iso_dat$Site)[3]

    nTotal <- nSite1+nSite2+nSite3

    originCap <- array(NA, c(nTotal,2))

    wherecaught <- rep(paste0("Site", 1:3), c(nSite1, nSite2, nSite3))

    originCap[which(pmatch(wherecaught,"Site1",dup = TRUE)==1),1] <- capCoords[1,1]
    originCap[which(pmatch(wherecaught,"Site1",dup = TRUE)==1),2] <- capCoords[1,2]

    originCap[which(pmatch(wherecaught,"Site2",dup = TRUE)==1),1] <- capCoords[2,1]
    originCap[which(pmatch(wherecaught,"Site2",dup = TRUE)==1),2] <- capCoords[2,2]

    originCap[which(pmatch(wherecaught,"Site3",dup = TRUE)==1),1] <- capCoords[3,1]
    originCap[which(pmatch(wherecaught,"Site3",dup = TRUE)==1),2] <- capCoords[3,2]

    originPoints <- sf::st_as_sf(as.data.frame(originCap), crs = 4326, coords = 1:2)
    sf_use_s2(TRUE)

    MC.realized <- calcMC(originDist = originDist,
                          targetDist = targetDist,
                          originRelAbund = originRelAbund,
                          psi = sim_psi_realized,
                          sampleSize=nTotal)



    originPointDists <- distFromPos(originCap)
    targetPointDists <- distFromPos(cbind(iso_dat$xcoords, iso_dat$ycoords))


    simAssign <- isoAssign(isovalues = iso_dat$featherIso.GrowingSeasonD,
                           isoSTD = 12,
                           intercept = -17.57,
                           slope = 0.95,
                           odds = 0.67,
                           restrict2Likely = FALSE,
                           nSamples = 500,
                           sppShapefile = OVEN_single_poly,
                           relAbund = relativeAbund,
                           verbose = 2,
                           isoWeight = -0.7,
                           abundWeight = 0,
                           assignExtent = c(-179,-60,15,89),
                           element = "Hydrogen",
                           period = "GrowingSeason")
    psi <- estTransition(targetRaster = simAssign,
                        targetSites = targetSites,
                        originPoints = originPoints,
                        originSites = originSites, isRaster = TRUE,
                        nSamples = 100, nSim = 5, verbose = 0)
    rM <- estMantel(targetRaster = simAssign,
                    targetSites = targetSites,
                    originPoints = originPoints, 
                    isGL = FALSE, isTelemetry = FALSE,
                    originSites = originSites, isRaster = TRUE,
                    nBoot = 100, nSim = 5, verbose = 0, captured = "origin")
    simEst <- estStrength(originRelAbund = originRelAbund,
                          targetDist = targetDist,
                          psi = psi,
                          originDist = originDist,
                          nSamples = 100,
                          verbose = 0,
                          alpha = 0.05)

    sim.output$targetSetup[sim] <- target
    sim.output$sim[sim]<-sim
    sim.output$MC.generated[sim] <- MC.generated
    sim.output$MC.realized[sim] <- MC.realized
    sim.output$MC.est[sim] <- simEst$MC$mean
    sim.output$MC.low[sim] <- simEst$MC$bcCI[1]
    sim.output$MC.high[sim] <- simEst$MC$bcCI[2]
    sim.output$rM.realized[sim] <- calcMantel(originDist = originPointDists, 
                                              targetDist = targetPointDists)$pointCorr
    sim.output$rM.est[sim] <- rM$corr$mean
    sim.output$rM.low[sim] <- rM$corr$bcCI[1]
    sim.output$rM.high[sim] <- rM$corr$bcCI[2]

  }
  output.sims[[target]] <- sim.output
}
Sys.time()-a

output.sims <- do.call("rbind", output.sims)
#

Summarize the output

sim.output <- transform(output.sims,
                        MC.generated.cover = as.integer((MC.low <= MC.generated & 
                                                           MC.high >= MC.generated)),
                        MC.realized.cover = as.integer((MC.low <= MC.realized & 
                                                          MC.high >= MC.realized)),
                        MC.generated.error = MC.est - MC.generated, 
                        MC.realized.error = MC.est - MC.realized,
                        rM.cover = as.integer((rM.low <= rM.realized & 
                                                 rM.high >= rM.realized)),
                        rM.error = rM.est - rM.realized)

summary(sim.output)
# Examine results
aggregate(MC.generated.error ~ targetSetup, sim.output, mean)
aggregate(MC.generated.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(MC.generated.cover ~ targetSetup, sim.output, mean)
aggregate(MC.realized.error ~ targetSetup, sim.output, mean)
aggregate(MC.realized.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(MC.realized.cover ~ targetSetup, sim.output, mean)
aggregate(I(MC.realized - MC.generated) ~ targetSetup, sim.output, mean)
aggregate(rM.error ~ targetSetup, sim.output, mean)
aggregate(rM.error ~ targetSetup, sim.output, function (x) mean(abs(x)))
aggregate(rM.cover ~ targetSetup, sim.output, mean)
aggregate(MC.est ~ targetSetup, sim.output, mean)