The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
Load required packages
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.
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$WGS84To 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.
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)))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 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)))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)These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.