docs/simulateIsoMC.R
In SMBC-NZP/MigConnectivity: Estimate Migratory Connectivity for Migratory Animals
## ----echo = FALSE----------------------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE,
comment = "#>",
tidy = FALSE,
message = FALSE,
warning = FALSE)
## ----eval = FALSE----------------------------------------------------------------------------
# library(terra)
# library(MigConnectivity)
# library(sf)
## ----eval = FALSE----------------------------------------------------------------------------
# # 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
## ----eval = FALSE----------------------------------------------------------------------------
# ### get isotope raster for template
# rasterTemplate <- getIsoMap()
## ----eval = FALSE----------------------------------------------------------------------------
# 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)))
## ----eval = FALSE----------------------------------------------------------------------------
# # 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)
## ----eval = FALSE----------------------------------------------------------------------------
# #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)))
## ----eval = FALSE----------------------------------------------------------------------------
# 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)
# #
## ----eval = FALSE----------------------------------------------------------------------------
# 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)
SMBC-NZP/MigConnectivity documentation built on March 26, 2024, 4:22 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----echo = FALSE----------------------------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE,
comment = "#>",
tidy = FALSE,
message = FALSE,
warning = FALSE)
## ----eval = FALSE----------------------------------------------------------------------------
# library(terra)
# library(MigConnectivity)
# library(sf)
## ----eval = FALSE----------------------------------------------------------------------------
# # 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
## ----eval = FALSE----------------------------------------------------------------------------
# ### get isotope raster for template
# rasterTemplate <- getIsoMap()
## ----eval = FALSE----------------------------------------------------------------------------
# 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)))
## ----eval = FALSE----------------------------------------------------------------------------
# # 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)
## ----eval = FALSE----------------------------------------------------------------------------
# #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)))
## ----eval = FALSE----------------------------------------------------------------------------
# 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)
# #
## ----eval = FALSE----------------------------------------------------------------------------
# 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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.