# R/utilityFunctions.R In MigConnectivity: Estimate Migratory Connectivity for Migratory Animals

#### Documented in distFromPos

```# Function for generating random transition probability matrix (expression of parametric uncertainty)
makePsiRand <- function(model, origin.set, target.set) {
vcv1 <- model\$results\$beta.vcv
beta1 <- model\$results\$beta\$estimate
new.beta <- MASS::mvrnorm(1, beta1, vcv1)
full.beta <- rep(NA, length(model\$results\$beta\$estimate))
new.real <- RMark::get.real(model, "Psi", new.beta,
design = model\$design.matrix, vcv=FALSE, se=FALSE)
mat <- RMark::TransitionMatrix(new.real)[origin.set, target.set]
return(mat)
}

# Calculates probability matrix based on exponential decline with distance
mlogitMat <- function(slope, dist) {
preMat <- exp(-slope/mean(dist)*dist)
diag(preMat) <- 0
nr <- nrow(dist)
nc <- ncol(dist)
outMat <- matrix(0, nr, nc)
for (b in 1:nr) {
outMat[b,] <- preMat[b,]/(1+sum(preMat[b, ]))
outMat[b,b] <- 1 - sum(outMat[b, ])
}
return(outMat)
}

# Calculates row of probability matrix based on exponential decline with distance
mlogitRow <- function(slope, dist, row.num) {
preMat <- exp(-slope/mean(dist)*dist)
diag(preMat) <- 0
out.row <- preMat[row.num,]/(1+sum(preMat[row.num, ]))
out.row[row.num] <- 1 - sum(out.row)
return(out.row)
}

# Crude optimizable function for developing MC pattern based on MC strength
mlogitMC <- function(slope, MC.in, origin.dist, target.dist, origin.rel.abund) {
nBreeding <- nrow(origin.dist)
nWintering <- nrow(target.dist)
psi <- mlogitMat(slope, origin.dist)
if (any(psi<0))
return(5*slope^2)
MC <- calcMC(origin.dist, target.dist, psi, origin.rel.abund)
return((MC.in - MC)^2)
}

# Samples from GL and/or GPS locations
# Called by estMantel and estMCGlGps
targetSample <- function(isGL,
geoBias,
geoVCov,
targetPoints,
animal.sample,
nSim = 1000,
targetSites = NULL,
targetAssignment = NULL,
resampleProjection = 'ESRI:53027',
maxTries = 300) {

# Target points should come in as sf object #
# length(targetPoints) works for sp objects #
# nrow(targetPoints) works for sf objects #
nAnimals <- nrow(targetPoints)

if (is.null(targetAssignment)) {
if (!is.null(targetSites)){
#Use the provided spatial layers to create targetAssignment
if(!sf::st_crs(targetSites)==sf::st_crs(targetPoints)){targetPoints <- sf::st_transform(targetPoints, crs = resampleProjection)}
targetAssignment <- suppressMessages(as.numeric(unclass(sf::st_intersects(x = targetPoints, y = targetSites, sparse = TRUE))))
#if targetAssignment is NA - convert to 0
targetAssignment[is.na(targetAssignment)] <- 0
}else
targetAssignment <- rep(0, nAnimals)
}

# Storage
target.sample <- rep(NA, nAnimals)
target.point.sample <- matrix(NA, nAnimals, 2)
# Fill in GPS values (no location uncertainty)
target.sample[which(!isGL[animal.sample])] <- targetAssignment[animal.sample[which(!isGL[animal.sample])]]
target.point.sample[which(!isGL[animal.sample]), ] <- sf::st_coordinates(targetPoints[animal.sample[which(!isGL[animal.sample])],])
# Which to sample still
toSample <- which(isGL[animal.sample])
# In this case, only have to sample each GL animal once, because no restrictions
if (is.null(targetSites) && any(isGL[animal.sample])) {
draws <- 0
while (length(toSample) > 0 && (is.null(maxTries) || draws <= maxTries)) {
draws <- draws + 1
geoBias2 <- array(rep(geoBias, length(toSample), each = nSim), c(nSim, 2, length(toSample)))

# generate sample and substract the bias
point.sample <- array(apply(sf::st_coordinates(targetPoints)[animal.sample[toSample], , drop = FALSE],
MARGIN = 1,
MASS::mvrnorm, n=nSim, Sigma=geoVCov),
# dimensions of the array
c(nSim, 2, length(toSample))) - geoBias2 #subtract the bias

# name so can
dimnames(point.sample)[[2]] <- c("x","y")

point.sample <- apply(point.sample, FUN = function(x){
sf::st_as_sf(data.frame(x), coords = c("x","y"),
crs = resampleProjection)}, MARGIN = 3)
good.sample <- apply(sapply(point.sample, sf::st_is_valid), 2,
function(x) min(which(x)))
if (any(!is.na(good.sample)))
target.point.sample[toSample[!is.na(good.sample)], ]<- t(mapply(x = good.sample[!is.na(good.sample)],
y = point.sample[!is.na(good.sample)],
FUN = function(x, y){sf::st_coordinates(y[x,])}))

toSample <- which(is.na(target.point.sample[,1]))
}
}else {
draws <- 0
while (length(toSample) > 0 && (is.null(maxTries) || draws <= maxTries)) {
draws <- draws + 1
# Create giant version of geoBias we can subtract from point.sample
geoBias2 <- array(rep(geoBias, length(toSample), each = nSim), c(nSim, 2, length(toSample)))
# Draw random points for GL animals from multivariate normal distribution
point.sample <- array(apply(sf::st_coordinates(targetPoints)[animal.sample[toSample], , drop = FALSE], 1,
MASS::mvrnorm, n=nSim, Sigma=geoVCov),
c(nSim, 2, length(toSample))) - geoBias2

dimnames(point.sample)[[2]] <- c("x","y")

# Convert those to sf

point.sample <- apply(point.sample, FUN = function(x){
sf::st_as_sf(data.frame(x), coords = c("x","y"),
crs = resampleProjection)}, MARGIN = 3)
#point.sample <- lapply(point.sample1, st_as_sf)
# Find out which sampled points are in a target site
if(!sf::st_crs(targetSites)==sf::st_crs(point.sample)){
targetSites <- sf::st_transform(targetSites, crs = resampleProjection)
}
target.sample0 <- sapply(point.sample, FUN = function(z){
inter <- sf::st_intersects(x = z, y = targetSites,
sparse = TRUE)
inter[lengths(inter)==0] <- NA
inter[lengths(inter)>1] <- sapply(inter[lengths(inter)>1],
function(x) x[1])
as.numeric(unlist(unclass(inter)))})
# Identify which animals have at least one valid sample point. good.sample
# will be NA for those that don't.  For those that do, it will location in
# point.sample first valid point can be found.
good.sample <- apply(target.sample0, 2, function(x) which(!is.na(x))[1])
# Fill in target site of valid sampled points
target.sample[toSample] <- apply(target.sample0, 2, function(x) x[!is.na(x)][1])
# Fill in target points of valid sampled points
if (any(!is.na(good.sample)))
target.point.sample[toSample[!is.na(good.sample)], ]<- t(mapply(x = good.sample[!is.na(good.sample)],
y = point.sample[!is.na(good.sample)],
FUN = function(x, y){sf::st_coordinates(y[x,])}))

toSample <- which(is.na(target.sample))
}
if (!is.null(maxTries) && draws > maxTries)
stop(paste0('maxTries (',maxTries,') reached during point resampling, exiting. Examine targetSites, geoBias, and geoVcov to determine why so few resampled points fall within targetSites.'))
}
return(list(target.sample = target.sample, target.point.sample = target.point.sample,
draws = draws))
}

# Function for sampling from geolocator, telemetry, and/or intrinsic location
# uncertainty on either origin or target side
locSample <- function(isGL,
isRaster,
isProb,
isTelemetry,
geoBias = NULL,
geoVCov = NULL,
points = NULL,
matvals = NULL,
matvals_crs = NULL,
singleCell = NULL,
pointsInSites = FALSE,
overlap1 = NULL,
assignment = NULL,
sites = NULL,
resampleProjection = 'ESRI:53027',
nSim = 1000,
maxTries = 300) {
#cat("Starting locSample\n")
# points should come in as sf object #
nAnimals <- length(isGL)
if (!is.null(singleCell)) {
# DETERMINE THE NUMBER OF RANDOM DRAWS FROM ISOSCAPE used in isoAssign

# How does targetIntrinsic make it into this function?
nrandomDraws <- dim(singleCell)[1]

# Stack 3D array into 2D, to make it easier to get different random point samples for each animal
intrinsic1 <- array(aperm(singleCell, c(1, 3, 2)),
dim = c(nAnimals * nrandomDraws, 2),
dimnames = list(NULL, c("x", "y")))

# project points to WGS #
intrinsic2 <- sf::st_as_sf(as.data.frame(intrinsic1), coords = c("x", "y"),
crs = 4326)
}

# If assignment isn't defined, create one from sites and points
if (is.null(assignment)) {
if (!is.null(sites) && !is.null(points)){
#Use the provided spatial layers to create assignment
if(!sf::st_crs(sites)==sf::st_crs(points)){
points <- sf::st_transform(points, crs = resampleProjection)
sites <- sf::st_transform(sites, crs = resampleProjection)
}
assignment <- suppressMessages(unclass(sf::st_intersects(x = points,
y = sites,
sparse = TRUE)))
#if assignment is not there - convert to 0
assignment[sapply(assignment, function(x) length(x)==0)] <- 0
assignment <- array(unlist(assignment), nAnimals)
}
else
assignment <- array(0, nAnimals)
}

# Storage
site.sample <- rep(NA, nAnimals)
point.sample <- matrix(NA, nAnimals, 2, dimnames = list(NULL, c('x', 'y')))
# Fill in telemetry/GPS values (no location uncertainty)
if (length(dim(assignment))==2){
site.sample[which(isTelemetry)] <- apply(assignment[isTelemetry, ,
drop = FALSE], 1,
which.max)
if (is.list(site.sample)) {
site.sample[lengths(site.sample)==0] <- 0
site.sample <- unlist(site.sample)
}

}else{
site.sample[which(isTelemetry)] <- assignment[which(isTelemetry)]
site.sample[which(isTelemetry & is.na(site.sample))] <- 0
}
if (!is.null(points))
point.sample[which(isTelemetry), ] <- sf::st_coordinates(points[which(isTelemetry),])
# Which to sample still
toSample <- which(!isTelemetry)
draws <- 0
# Loop until all animals have valid sample or have exceeded maxTries
while (length(toSample) > 0 && (is.null(maxTries) || draws < maxTries)) {
draws <- draws + 1
# Convert toSample (numbers) into T/F, so can combine with other T/F
toSampleBool <- 1:nAnimals %in% toSample
# cat(toSampleBool, '\n')
# cat(isGL[1:10], '\n')
# Sample geolocator points
if (any(isGL & toSampleBool)) {
#cat("*************", sum(isGL & toSampleBool), "*************\n")
# Create giant version of geoBias we can subtract from point.sample
geoBias2 <- array(rep(geoBias, sum(isGL & toSampleBool), each = nSim),
c(nSim, 2, sum(isGL & toSampleBool)))
# Draw random points for GL animals from multivariate normal distribution
point.sample0 <- array(apply(sf::st_coordinates(points)[toSampleBool & isGL,
, drop = FALSE], 1,
MASS::mvrnorm, n=nSim, Sigma=geoVCov),
c(nSim, 2, sum(isGL & toSampleBool))) - geoBias2

dimnames(point.sample0)[[2]] <- c("x","y")

# Convert those to SpatialPoints

point.sample0 <- apply(point.sample0,
FUN = function(x){sf::st_as_sf(data.frame(x),
coords = c("x","y"),
crs = resampleProjection)},
MARGIN = 3)
if (!is.null(sites)) {
# Find out which sampled points are in a target site
if(!sf::st_crs(sites)==sf::st_crs(point.sample0[[1]])){
sites <- sf::st_transform(sites, crs = resampleProjection)
}
site.sample0 <- sapply(point.sample0, FUN = function(z){
inter <- sf::st_intersects(x = z, y = sites,
sparse = TRUE)
inter[lengths(inter)==0] <- NA
inter[lengths(inter)>1] <- sapply(inter[lengths(inter)>1], function(x) x[1])
as.numeric(unlist(unclass(inter)))})

# Identify which animals have at least one valid sample point. good.sample0
# will be NA for those that don't.  For those that do, it will location in
# point.sample first valid point can be found.
good.sample0 <- apply(site.sample0, 2, function(x) which(!is.na(x))[1])
}
else {
site.sample0 <- array(0, c(nSim, length(point.sample0)))
good.sample0 <- rep(1, length(point.sample0))
}
# cat(good.sample0, "\n")
}
if (any(isProb & toSampleBool)) {
site.sample1 <- apply(assignment[isProb & toSampleBool, , drop = FALSE],
1, sample.int, n = ncol(assignment), size = nSim,
replace = TRUE)
if (is.null(dim(site.sample1)))
site.sample1 <- matrix(site.sample1, 1)
}
# walk through this section #
if (any(isRaster & toSampleBool)) {
# some raster -
if (pointsInSites && !any(isRaster & toSampleBool & (isGL | isProb))) {
# IF ALL THE POINTS ARE WITHIN SITES and no animals have two data types
samp <- sample.int(nrandomDraws,
size = sum(isRaster & toSampleBool),
replace = TRUE)
# GRAB LOCATIONS FROM RASTER #
samp2 <- samp + ((1:nAnimals)[toSampleBool & isRaster] - 1) * nrandomDraws

point.sample2 <- intrinsic2[samp2, ]

# Changed to make sure x,y coords stack correctly
# target.point.sample[toSample,1]<- point.sample@coords[,1]
# target.point.sample[toSample,2]<- point.sample@coords[,2]

# ANIMALS WITHIN samp2 that don't have raster
# We ULTIMATELY NEED POINTS BECAUSE THEY ONLY HAVE RASTER DATA
if (!is.null(overlap1))
# Grab the relevant sites
site.sample2 <- overlap1[samp, (1:nAnimals)[toSampleBool & isRaster], drop = FALSE]
good.sample2 <- seq(from = 1, by = nSim,
length.out = sum(isRaster & toSampleBool))
}
else if (is.null(singleCell)) {
multidraw <- apply(matvals[ , (1:nAnimals)[toSampleBool & isRaster] + 2, drop = FALSE],
2, function(x) (rmultinom(prob = x, n = nSim, size = 1)))
multidraw <- array(multidraw, c(nrow(matvals), nSim,
sum(toSampleBool & isRaster)))
#point.sample2 <- matvals[which(multidraw == 1, arr.ind = TRUE)[, 1], 1:2]
point.sample2 <- apply(multidraw, c(2:3), function(x) which(x==1))
point.sample2 <- array(as.matrix(matvals[c(point.sample2), 1:2]),
c(nSim * sum(toSampleBool & isRaster), 2),
list(NULL, c("x","y")))
#dimnames(point.sample2)[[2]] <- c("x","y")

# Convert those to sf
point.sample2 <- sf::st_as_sf(as.data.frame(point.sample2),
coords = c("x","y"),
#crs = 4326)
crs = sf::st_crs(matvals_crs))

# point.sample2 <- apply(point.sample2,
#                        FUN = function(x){
#                          sf::st_as_sf(data.frame(x),
#                                       coords = c("x","y"),
#                             crs = 4326)},
#                        MARGIN = 2)
#point.sample <- lapply(point.sample1, st_as_sf)
if (!is.null(sites)) {
# Find out which sampled points are in a target site
if(!sf::st_crs(sites)==sf::st_crs(point.sample2)){
sites <- sf::st_transform(sites, crs = resampleProjection)
point.sample2 <- sf::st_transform(point.sample2,
crs = resampleProjection)
}

# site.sample2 <- sapply(point.sample2, FUN = function(z){
#   suppressMessages(as.numeric(unclass(sf::st_intersects(x = z, y = sites,
#                                        sparse = TRUE))))})
# Check which points are in target sites
site.sample2 <- suppressMessages(sf::st_intersects(point.sample2,
y = sites))

# Give values without intersection an NA
len_intersect <- lengths(site.sample2)
# quick check to ensure that all the points fall exactly in one site #
if(any(len_intersect)>1){
warning("Overlapping targetSites or originSites may cause issues\n")
site.sample2 <- lapply(site.sample2, function (x) x[1])
}
site.sample2[lengths(site.sample2)==0] <- NA
site.sample2 <- unlist(as.numeric(site.sample2))

# Organize into matrix (separate by animal)
site.sample2 <- matrix(site.sample2, nSim, sum(isRaster & toSampleBool))

# Identify which animals have at least one valid sample point. good.sample2
# will be NA for those that don't.  For those that do, it will location in
# point.sample first valid point can be found.
good.sample2 <- apply(site.sample2, 2, function(x) which(!is.na(x))[1])+
seq(from = 0, by = nSim, length.out = sum(isRaster & toSampleBool))
#cat(good.sample2, "\n")
}
else {
site.sample2 <- array(0, c(nSim, sum(isRaster & toSampleBool)))
good.sample2 <- rep(1, sum(isRaster & toSampleBool))
}
}
else {
# Select nSim points for each animal still to be sampled
samp <- sample.int(nrandomDraws,
size = sum(isRaster & toSampleBool) * nSim,
replace = TRUE)

samp2 <- samp + rep((1:nAnimals)[toSampleBool & isRaster] - 1, each = nSim) * nrandomDraws

point.sample2 <- intrinsic2[samp2, ]
if (!is.null(sites)) {
if(!sf::st_crs(sites)==sf::st_crs(point.sample2)){
point.sample2 <- sf::st_transform(point.sample2,
crs = resampleProjection)
}
# Check which points are in target sites
site.sample2 <- suppressMessages(sf::st_intersects(point.sample2,
y = sites))

# Give values without intersection an NA
len_intersect <- lengths(site.sample2)
# quick check to ensure that all the points fall exactly in one site #
if(any(len_intersect)>1){
warning("Overlapping targetSites or originSites may cause issues\n")
site.sample2 <- lapply(site.sample2, function (x) x[1])
}
site.sample2[len_intersect==0] <- NA
site.sample2 <- unlist(as.numeric(site.sample2))

# Organize into matrix (separate by animal)
site.sample2 <- matrix(site.sample2, nSim, sum(isRaster & toSampleBool))
# Identify which animals have a valid point (inside a site).
# good.sample2, for those that don't, will be NA. For those that do, it
# will be location in point.sample where first valid point is.

good.sample2 <- apply(site.sample2, 2, function(x){which(!is.na(x))[1]})+
seq(from = 0, by = nSim, length.out = sum(isRaster & toSampleBool))
}
else {
site.sample2 <- array(0, c(nSim, sum(isRaster & toSampleBool)))
good.sample2 <- rep(1, sum(isRaster & toSampleBool))
}
}
}
if (any(isProb & !isGL & !isRaster & toSampleBool)){
site.sample[isProb & !isGL & !isRaster & toSampleBool] <-
site.sample1[1, !(which(isProb & toSampleBool) %in% which(isGL | isRaster)), drop = FALSE]
}
if (any(!isProb & isGL & !isRaster & toSampleBool)){
if (any(!is.na(good.sample0[!(which(isGL) %in% which(isProb | isRaster))]))){
site.sample[!isProb & isGL & !isRaster & toSampleBool] <-
apply(site.sample0[, !(which(isGL & toSampleBool) %in% which(isProb | isRaster)), drop = FALSE],
2,
function(x) x[!is.na(x)][1])
# Fill in target points of valid sampled points
point.sample[which(!isProb & isGL & !isRaster & toSampleBool)[which(!is.na(good.sample0[!(which(isGL & toSampleBool) %in% which(isProb | isRaster))]))], ]<-
t(mapply(x = good.sample0[!is.na(good.sample0) & !(which(isGL & toSampleBool) %in% which(isProb | isRaster))],
y = point.sample0[which(!(which(isGL & toSampleBool) %in% which(isProb | isRaster)) & (!is.na(good.sample0)))],
FUN = function(x, y) sf::st_coordinates(y[x,])))
}
}
if (any(isProb & isGL & !isRaster & toSampleBool)){
if (any(!is.na(good.sample0))){
compare <- site.sample0[, which(isGL & toSampleBool) %in% which(isProb & !isRaster),
drop = FALSE] ==
site.sample1[, which(isProb & toSampleBool) %in% which(isGL & !isRaster),
drop = FALSE]
good.sample01 <- apply(compare, 2, function(x) which(!is.na(x) & x)[1])
if (any(!is.na(good.sample01))){
for (col in 1:sum(isProb & isGL & !isRaster & toSampleBool))
site.sample[which(isProb & isGL & !isRaster & toSampleBool)[col]] <-
site.sample0[good.sample01[col],
which(which(isGL & toSampleBool) %in% which(isProb & !isRaster))[col]]
# Fill in target points of valid sampled points
# point.sample[which(isProb & isGL & !isRaster & toSampleBool)[which(!is.na(good.sample01))], ]<-
#   t(mapply(x = good.sample01[!is.na(good.sample01)],
#            y = point.sample0[which(which(isGL & toSampleBool) %in% which(isProb & !isRaster))[which(!is.na(good.sample01))]],
#            FUN = function(x, y) sf::st_coordinates(y[x,])))
}
}
}
if (any(!isProb & !isGL & isRaster & toSampleBool)){
if (any(!is.na(good.sample2[!(which(isRaster) %in% which(isProb | isGL))]))){
site.sample[which(!isProb & !isGL & isRaster & toSampleBool)] <-
apply(site.sample2[, !(which(isRaster & toSampleBool) %in% which(isProb | isGL)),
drop = FALSE],
2,
function(x) x[!is.na(x)][1])
rows.select <- good.sample2[!(which(isRaster & toSampleBool) %in%
which(isProb | isGL) |
is.na(good.sample2))]
# Fill in target points of valid sampled points
point.sample[which(!isProb & !isGL & isRaster & toSampleBool)[which(!is.na(good.sample2))], ] <-
sf::st_coordinates(point.sample2[rows.select, ])
#[which(!is.na(good.sample2))]
# t(mapply(x = good.sample2[!(which(isRaster & toSampleBool) %in% which(isProb | isGL))],
#          y = point.sample2[which(which(isRaster & toSampleBool) %in% which(!isProb & !isGL & !is.na(good.sample2))),,
#                            drop = FALSE],
#          FUN = function(x, y) sf::st_coordinates(y[x,])))
}
}
if (any(isProb & !isGL & isRaster & toSampleBool)){
if (any(!is.na(good.sample2))){
compare <- site.sample2[, which(isRaster & toSampleBool) %in% which(isProb & !isGL),
drop = FALSE] ==
site.sample1[, which(isProb & toSampleBool) %in% which(!isGL & isRaster),
drop = FALSE]
good.sample12 <- apply(compare, 2, function(x) which(!is.na(x) & x)[1])
if (any(!is.na(good.sample12))){
for (col in 1:sum(isProb & !isGL & isRaster & toSampleBool))
site.sample[which(isProb & !isGL & isRaster & toSampleBool)[col]] <-
site.sample2[good.sample12[col],
which(which(isRaster & toSampleBool) %in% which(isProb & !isGL))[col]]
# Fill in target points of valid sampled points
# Doesn't get used, not going to include prob with estMantel
# point.sample[which(isProb & !isGL & isRaster & toSampleBool)[which(!is.na(good.sample12))], ]<-
#   t(mapply(x = good.sample12[!is.na(good.sample12)],
#            y = point.sample2[which(which(isRaster & toSampleBool) %in% which(isProb & !isGL))[which(!is.na(good.sample12))], ],
#            FUN = function(x, y) sf::st_coordinates(y[x,])))
}
}
}
if (any(!isProb & isGL & isRaster & toSampleBool)){
if (any(!is.na(good.sample0)) && any(!is.na(good.sample2))){
compare <- site.sample2[, which(isRaster & toSampleBool) %in% which(!isProb & isGL),
drop = FALSE] ==
site.sample0[, which(isGL & toSampleBool) %in% which(!isProb & isRaster),
drop = FALSE]
good.sample02 <- apply(compare, 2, function(x) which(!is.na(x) & x)[1])
#cat(good.sample02, '\n')
if (any(!is.na(good.sample02))){
for (col in 1:sum(!isProb & isGL & isRaster & toSampleBool))
site.sample[which(!isProb & isGL & isRaster & toSampleBool)[col]] <-
site.sample2[good.sample02[col],
which(which(isRaster & toSampleBool) %in% which(!isProb & isGL))[col]]
# Just choosing GL point instead of fussing with raster point - generally
# those are a lot more precise so probably also more accurate.
point.sample[which(!isProb &
isGL &
isRaster &
toSampleBool)[which(!is.na(good.sample02))],] <-
t(mapply(
x = good.sample02[!is.na(good.sample02)],
y = point.sample0[which(which(isGL & toSampleBool) %in%
which(!isProb & isRaster))[which(!is.na(good.sample02))]],
FUN = function(x, y)
sf::st_coordinates(y[x, ])
))
}
}
}
if (any(isProb & isGL & isRaster & toSampleBool)){
if (any(!is.na(good.sample0)) && any(!is.na(good.sample2))){
compare <- (site.sample2[, which(isRaster & toSampleBool) %in% which(isProb & isGL),
drop = FALSE] ==
site.sample0[, which(isGL & toSampleBool) %in% which(isProb & isRaster),
drop = FALSE]) &
(site.sample1[, which(isProb & toSampleBool) %in% which(isRaster & isGL),
drop = FALSE] ==
site.sample0[, which(isGL & toSampleBool) %in% which(isProb & isRaster),
drop = FALSE])
good.sample012 <- apply(compare, 2, function(x) which(!is.na(x) & x)[1])
if (any(!is.na(good.sample012))){
for (col in 1:sum(isProb & isGL & isRaster & toSampleBool))
site.sample[which(isProb & isGL & isRaster & toSampleBool)[col]] <-
site.sample2[good.sample012[col],
which(which(isRaster & toSampleBool) %in% which(isProb & isGL))[col]]
# This one might not matter (not planning to include isProb in
# estMantel or look at points in estTransition).
# point.sample[which(isProb & isGL & isRaster & toSampleBool)[which(!is.na(good.sample012))]]<-
#   t(mapply(x = good.sample012[!is.na(good.sample012)],
#            y = point.sample0[which(which(isGL & toSampleBool) %in% which(isProb & isRaster))[which(!is.na(good.sample012))]],
#            FUN = function(x, y) sf::st_coordinates(y[x,])))
}
}
}
toSample <- which(is.na(site.sample))
}
if (length(toSample) > 0){
notfind <- data.frame(Animal = toSample, isGL = isGL[toSample],
isRaster = isRaster[toSample], isProb = isProb[toSample],
isTelemetry = isTelemetry[toSample])
return(list(site.sample = site.sample, point.sample = point.sample,
draws = draws, notfind = notfind))
}
return(list(site.sample = site.sample, point.sample = point.sample,
draws = draws))
}

# Samples from isotope target points
# Called by estMCisotope
targetSampleIsotope <- function(targetIntrinsic, animal.sample,
nSim = 10, targetSites = NULL,
resampleProjection = "+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs",
maxTries = 300, pointsAssigned = FALSE,
targCon = NULL, pointsInSites = FALSE) {
# Here is a catch to save Mike from himself
#if(is.null(nSim)){nSim<-1000}
if (pointsAssigned) {
nAnimals <- dim(targetIntrinsic\$SingleCell)[3]
# DETERMINE THE NUMBER OF RANDOM DRAWS FROM ISOSCAPE used in isoAssign
nrandomDraws <- dim(targetIntrinsic\$SingleCell)[1]

# Stack 3D array into 2D, to make it easier to get different random point samples for each animal
targetIntrinsic1 <- array(aperm(targetIntrinsic\$SingleCell, c(1, 3, 2)),
dim = c(nAnimals * nrandomDraws, 2),
dimnames = list(NULL, c("x","y")))

# project target points to WGS #
targetIntrinsic2 <- sf::st_as_sf(as.data.frame(targetIntrinsic1), coords = c("x", "y"),
crs = 4326)

}
else {
nAnimals <- dim(targetIntrinsic\$probassign)[3]
}
if (!pointsInSites){
if (!inherits(targetIntrinsic\$probassign, "SpatRaster"))
targetIntrinsic\$probassign <- terra::rast(targetIntrinsic\$probassign)
# converts raster to matrix of XY then probs
matvals <- terra::as.data.frame(targetIntrinsic\$probassign, xy = TRUE)
}
target.sample <- rep(NA, nAnimals)
target.point.sample <- matrix(NA, nAnimals, 2)
toSample <- 1:nAnimals
# In this case, only need to draw once, because already determined that points fall in targetSites
if (pointsInSites) {
draws <- 1
samp <- sample.int(nrandomDraws, size = length(toSample), replace = TRUE)
samp2 <- samp + (animal.sample[toSample] - 1) * nrandomDraws
point.sample <- targetIntrinsic2[samp2, ]
# Changed to make sure x,y coords stack correctly
#target.point.sample[toSample,1]<- point.sample@coords[,1]
#target.point.sample[toSample,2]<- point.sample@coords[,2]
target.point.sample[toSample,1] <- sf::st_coordinates(point.sample)[,1]
target.point.sample[toSample,2] <- sf::st_coordinates(point.sample)[,2]
if (!is.null(targCon))
# Grab the relevant target sites
target.sample[toSample] <- targCon[samp2]
}
# In this case, only need to draw once, as no limits on where points can be
else if (is.null(targetSites)) {
draws <- 1
# This draws samples nSamples per animal (faster than looping over nSamples) and fills the xysim with x,y coords
for(i in 3:ncol(matvals)) {
animals <- which(animal.sample==i - 2)
if (length(animals) > 0) {
multidraw <- rmultinom(n = length(animals), size = 1, prob = matvals[,i])
target.point.sample[animals,1] <- matvals[which(multidraw == 1, arr.ind = TRUE)[,1],1]
target.point.sample[animals,2] <- matvals[which(multidraw == 1, arr.ind = TRUE)[,1],2]
}
}
}
else {
draws <- 0
# Make sure targetSites are WGS84
targetSites <- sf::st_transform(targetSites, 4326)

while (length(toSample) > 0 && (is.null(maxTries) || draws <= maxTries)) {
draws <- draws + 1
if (!pointsAssigned) {
# This draws samples nSamples per animal (faster than looping over nSamples) and fills the xysim with x,y coords
for(i in 3:ncol(matvals)) {
animals <- which(animal.sample[toSample]==i - 2)
if (length(animals) > 0) {
multidraw <- rmultinom(n = length(animals)*nSim, size = 1, prob = matvals[,i])

point.sample <- matvals[which(multidraw == 1, arr.ind = TRUE)[, 1], 1:2]

# covert to sf
point.sample1 <- sf::st_as_sf(as.data.frame(point.sample),
coords = c("x","y"),
crs = 4326)

# Check which points are in target sites

target.sample0 <- suppressMessages(as.numeric(unclass(sf::st_intersects(x = point.sample1,
y = targetSites,
sparse = TRUE))))

good.sample <- which(!is.na(target.sample0))

nToFill <- min(length(good.sample), length(animals))
if (nToFill > 0) {
target.sample[animals[1:nToFill]] <- target.sample0[good.sample[1:nToFill]]
target.point.sample[animals[1:nToFill], 1] <- point.sample[good.sample[1:nToFill], 1]
target.point.sample[animals[1:nToFill], 2] <- point.sample[good.sample[1:nToFill], 2]
}
}
}
}
else {
# Select nSim points for each animal still to be sampled
samp <- sample.int(nrandomDraws, size = length(toSample) * nSim,
replace = TRUE)
samp2 <- samp + rep(animal.sample[toSample] - 1, each = nSim) * nrandomDraws
point.sample <- targetIntrinsic2[samp2,]
# Check which points are in target sites
target.sample0 <- suppressMessages(as.numeric(unclass(sf::st_intersects(x = point.sample,
y = targetSites,
sparse = TRUE))))
# Organize into matrix (separate by animal)
target.sample1 <- matrix(target.sample0, nSim, length(toSample))
# Identify which animals have a valid point (inside a target site).
# good.sample, for those that don't, will be NA. For those that do, it
# will be location in point.sample where first valid point is.
good.sample <- apply(target.sample1, 2, function(x) which(!is.na(x))[1]) +
seq(from = 0, by = nSim, length.out = length(toSample))
# Put in target sites of animals with valid points sampled
target.sample[toSample] <- apply(target.sample1, 2, function(x) x[!is.na(x)][1])
# Put in target points of animals with valid points sampled
if (any(!is.na(good.sample))) {
target.point.sample[toSample[!is.na(good.sample)],1]<- sf::st_coordinates(point.sample[good.sample[!is.na(good.sample)],])[,1]
target.point.sample[toSample[!is.na(good.sample)],2]<- sf::st_coordinates(point.sample[good.sample[!is.na(good.sample)],])[,2]
}
}
toSample <- which(is.na(target.sample))
}
if (!is.null(maxTries) && draws > maxTries)
stop(paste0('maxTries (',maxTries,') reached during point resampling, exiting. Examine targetSites, geoBias, and geoVcov to determine why so few resampled points fall within targetSites.'))
}
return(list(target.sample = target.sample, target.point.sample = target.point.sample,
draws = draws))
}

#' Distance matrix from position matrix
#'
#' @param pos Number of sites by 2 matrix with positions of each site.  If
#'    \code{surface} is 'ellipsoid' or 'sphere', then column 1 should be
#'    longitude and column 2 should be latitude. If \code{surface} is 'plane',
#'    column 1 can be x-position and column 2 y-position.
#' @param surface Surface to calculate distances on.  Either 'ellipsoid'
#'    (default), 'sphere', or 'plane'.
#' @param units Units of return distance matrix. If \code{surface} is 'plane',
#'    then this argument is ignored and the return units will be the same as the
#'    \code{pos} units. Options are 'km' (kilometers, default), 'm' (meters),
#'    'miles', and 'nautical miles'.
#'
#' @return Square matrix of distances between sites. If \code{surface} is
#'    'ellipsoid' or 'sphere', then argument \code{units} will determine units;
#'    if \code{surface} is 'plane', the units will be the same as the \code{pos}
#'    units.
#'
#' @note In version 0.4.3 we switched package dependencies from \code{geosphere}
#'    to \code{geodist}. As a result, spherical distances (and possibly
#'    ellipsoid distances) may differ slightly from those calculated with earlier
#'    versions of our package.
#'
#' @export
#'
#' @examples
#' nBreeding <- 100
#' nWintering <- 100
#' breedingPos <- matrix(c(rep(seq(-99, -81, 2), each = sqrt(nBreeding)),
#'                         rep(seq(49, 31, -2), sqrt(nBreeding))),
#'                       nBreeding, 2)
#' winteringPos <- matrix(c(rep(seq(-79, -61, 2), each = sqrt(nWintering)),
#'                          rep(seq(9, -9, -2), sqrt(nWintering))),
#'                        nWintering, 2)
#' tail(breedingPos)
#' tail(winteringPos)
#'
#' breedDist <- distFromPos(breedingPos, 'ellipsoid')
#' nonbreedDist <- distFromPos(winteringPos, 'ellipsoid')
#' breedDist[1:12, 1:12]
#' breedDist[1:12, c(1,91,100)]
distFromPos <- function(pos, surface = 'ellipsoid',
units = c("km", "m", "miles", "nautical miles")) {
if (!is.matrix(pos) || dim(pos)[2] != 2)
stop("pos should be a matrix with two column (e.g., longitude and latitude)")
units <- match.arg(units)
if (!(surface %in% c('plane', 'sphere', 'ellipsoid')))
stop('surface must be "plane", "sphere", or "ellipsoid".')
if (is.null(colnames(pos)) && surface != 'plane')
colnames(pos) <- c("lon", "lat")
div <- ifelse(units == "km", 1000,
ifelse(units == "m", 1,
ifelse(units == "miles", 1609.34, 1852)))
if (surface == 'ellipsoid')
dist <- geodist::geodist(pos, measure = "geodesic") / div
else if (surface == 'sphere')
dist <- geodist::geodist(pos, measure = "vincenty") / div
else {
nSites <- nrow(pos)
dist <- matrix(0, nSites, nSites)
for (b1 in 2:nSites) {
for (b2 in 1:(b1-1)) {
dist[b1, b2] <- sqrt((pos[b1, 1] - pos[b2, 1]) ^ 2 +
(pos[b1, 2] - pos[b2, 2]) ^ 2)
dist[b2, b1] <- dist[b1, b2]
}
}
}
return(dist)
}

reassignInds <- function(dataOverlapSetting = "none",
originPoints = NULL, targetPoints = NULL,
originAssignment = NULL, targetAssignment = NULL,
isGL = FALSE, isTelemetry = FALSE,
isRaster = FALSE, isProb = FALSE,
captured = "origin",
#originRaster = NULL,
originRasterXYZ = NULL,
originRasterXYZcrs = NULL,
originSingleCell = NULL,
#targetRaster = NULL,
targetRasterXYZ = NULL,
targetRasterXYZcrs = NULL,
targetSingleCell = NULL,
originSites = NULL,
targetSites = NULL) {
if (dataOverlapSetting != "none") {
stop('dataOverlapSetting "named" not set up yet')
}
if (is.null(targetPoints))
nTargetPoints <- 0
else
nTargetPoints <- nrow(targetPoints)
if (is.null(originPoints))
nOriginPoints <- 0
else
nOriginPoints <- nrow(originPoints)
if (any(isGL) || any(isTelemetry)) {
nPoints <- max(nTargetPoints, nOriginPoints)
}
else {
nPoints <- 0
}
# capturePointsGL <- T
# if (any(isGL))
#   if (sum(captured[isGL]=="origin")>nOriginPoints)
if (any(isRaster)) {
if (!is.null(originRasterXYZ)) {
dimORast <- dim(originRasterXYZ)
nORast <- dimORast[2] - 2
}
else
nORast <- 0
if (!is.null(targetRasterXYZ)) {
dimTRast <- dim(targetRasterXYZ)
nTRast <- dimTRast[2] - 2
}
else
nTRast <- 0
}
else {
nORast <- nTRast <- 0
}
if (any(isProb)) {
if (is.null(originAssignment)){
nOAss <- 0
dimOAss <- NULL
}
else {
dimOAss <- dim(originAssignment)
if (is.null(dimOAss)) {
nOAss <- length(originAssignment)
tt2 <- matrix(0, nOAss, max(length(unique(originAssignment)),
nrow(originSites)))
for (i in 1:nOAss)
tt2[i, originAssignment[i]] <- 1
originAssignment <- tt2
dimOAss <- dim(originAssignment)
}
else
nOAss <- dimOAss[1]
}
if (is.null(targetAssignment)){
nTAss <- 0
dimTAss <- NULL
}
else {
dimTAss <- dim(targetAssignment)
if (is.null(dimTAss)){
nTAss <- length(targetAssignment)
tt2 <- matrix(0, nTAss, max(length(unique(targetAssignment)),
nrow(targetSites)))
for (i in 1:nTAss)
tt2[i, targetAssignment[i]] <- 1
targetAssignment <- tt2
dimTAss <- dim(targetAssignment)
}
else
nTAss <- dimTAss[1]
}
}
nTotal <- max(length(isGL), length(isTelemetry), length(isRaster),
length(isProb))
if (nTotal < 2) {
nTotal <- length(captured)
if (nTotal < 2 &&
((any(isGL) + any(isTelemetry) + any(isRaster) + any(isProb)) > 1)) {
warning("For multiple datasets, we recommend having an entry in captured, isGL, isTelemetry, etc. for each animal")
nTotal <- nPoints + nORast + nTRast + max(nOAss, nTAss)
}
}
if (length(isGL)==1)
isGL <- rep(isGL, nTotal)
if (length(isTelemetry)==1)
isTelemetry <- rep(isTelemetry, nTotal)
if (length(isRaster)==1)
isRaster <- rep(isRaster, nTotal)
if (length(isProb)==1)
isProb <- rep(isProb, nTotal)
if (length(captured)==1)
captured <- rep(captured, nTotal)
capturePoint <- rep(FALSE, nTotal)
nTPremain <- nTargetPoints
nOPremain <- nOriginPoints
if (any(isGL)){
nTPremain <- nTPremain - sum(captured[isGL]!="target")
nOPremain <- nOPremain - sum(captured[isGL]!="origin")
if (nTPremain > 0) {
capturePoint[captured=="target" & isGL] <- TRUE
nTPremain <- nTPremain - sum(captured[isGL]=="target")
}
if (nOPremain > 0) {
capturePoint[captured=="origin" & isGL] <- TRUE
nOPremain <- nOPremain - sum(captured[isGL]=="origin")
}
}
if ((nTPremain > 0 || nOPremain > 0) && any(isTelemetry)) {
nTPremain <- nTPremain - sum(captured[isTelemetry]!="target")
nOPremain <- nOPremain - sum(captured[isTelemetry]!="origin")
if (nTPremain > 0) {
capturePoint[captured=="target" & isTelemetry] <- TRUE
nTPremain <- nTPremain - sum(captured[isTelemetry]=="target")
}
if (nOPremain > 0) {
capturePoint[captured=="origin" & isTelemetry] <- TRUE
nOPremain <- nOPremain - sum(captured[isTelemetry]=="origin")
}
}
if ((nTPremain > 0 || nOPremain > 0) && any(isRaster)) {
if (nTPremain > 0) {
capturePoint[captured=="target" & isRaster] <- TRUE
nTPremain <- nTPremain - sum(captured[isRaster]=="target")
}
if (nOPremain > 0) {
capturePoint[captured=="origin" & isRaster] <- TRUE
nOPremain <- nOPremain - sum(captured[isRaster]=="origin")
}
}
if ((nTPremain > 0 || nOPremain > 0) && any(isProb)) {
if (nTPremain > 0) {
capturePoint[captured=="target" & isProb] <- TRUE
nTPremain <- nTPremain - sum(captured[isProb]=="target")
}
if (nOPremain > 0) {
capturePoint[captured=="origin" & isProb] <- TRUE
nOPremain <- nOPremain - sum(captured[isProb]=="origin")
}
}
if (nTargetPoints > 0 || nOriginPoints > 0) {
if (is.null(targetPoints)) {
targetPoints2 <- NULL
}
else {
if (nTargetPoints < nTotal) {
dummyPoint <- targetPoints[1, ]
targetPoints2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isGL[i] || isTelemetry[i] || (captured[i] == "target" &&
capturePoint[i])) {
place <- place + 1
targetPoints2 <- rbind(targetPoints2,
targetPoints[place, ])
}
else
targetPoints2 <- rbind(targetPoints2,
dummyPoint)
}
}
else
targetPoints2 <- targetPoints
}
if (is.null(originPoints)) {
originPoints2 <- NULL
}
else {
if (nOriginPoints < nTotal) {
dummyPoint <- originPoints[1, ]
originPoints2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isGL[i] || isTelemetry[i] || captured[i] == "origin" &&
capturePoint[i]) {
place <- place + 1
originPoints2 <- rbind(originPoints2,
originPoints[place, ])
}
else
originPoints2 <- rbind(originPoints2,
dummyPoint)
}
}
else
originPoints2 <- originPoints
}
}
else {
originPoints2 <- originPoints; targetPoints2 <- targetPoints
}
if (any(isRaster)) {
if (!is.null(originRasterXYZ)) {
if (all(isRaster & captured != "origin")) {
originRasterXYZ2 <- originRasterXYZ
}
else {
originRasterXYZ2 <- originRasterXYZ[, 1:2]
column <- 2
for (i in 1:nTotal) {
if (isRaster[i] && captured[i] != "origin") {
column <- column + 1
originRasterXYZ2 <- cbind(originRasterXYZ2,
originRasterXYZ[, column])
}
else {
originRasterXYZ2 <- cbind(originRasterXYZ2,
array(1/dimORast[1], c(dimORast[1], 1)))
}
}
if (!is.null(originSingleCell)) {
dimOCell <- dim(originSingleCell)
dummyVals <- c(originSingleCell[ , , 1])
originSingleCell2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isRaster[i] && captured[i] != "origin") {
place <- place + 1
originSingleCell2 <- c(originSingleCell2, originSingleCell[,,place])
}
else
originSingleCell2 <- c(originSingleCell2, dummyVals)

}
originSingleCell <- array(originSingleCell2,
c(dimOCell[1], dimOCell[2], nTotal),
list(NULL, c("Longitude", "Latitude"), NULL))
}
}
}
else{
originRasterXYZ2 <- NULL
}
if (!is.null(targetRasterXYZ)) {
if (all(isRaster & captured != "target")) {
targetRasterXYZ2 <- targetRasterXYZ
}
else {
xvalues <- range(targetRasterXYZ[,1])
yvalues <- range(targetRasterXYZ[,2])
targetRasterXYZ2 <- targetRasterXYZ[, 1:2]
column <- 2
for (i in 1:nTotal) {
if (isRaster[i] && captured[i] != "target") {
column <- column + 1
targetRasterXYZ2 <- cbind(targetRasterXYZ2,
targetRasterXYZ[, column])
}
else{
targetRasterXYZ2 <- cbind(targetRasterXYZ2,
array(1/dimTRast[1], c(dimTRast[1], 1)))
}
}
if (!is.null(targetSingleCell)) {
dimTCell <- dim(targetSingleCell)
dummyVals <- c(targetSingleCell[ , , 1])
targetSingleCell2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isRaster[i] && captured[i] != "target") {
place <- place + 1
targetSingleCell2 <- c(targetSingleCell2, targetSingleCell[,,place])
}
else
targetSingleCell2 <- c(targetSingleCell2, dummyVals)

}
targetSingleCell <- array(targetSingleCell2,
c(dimTCell[1], dimTCell[2], nTotal))
}
}
}
else {
targetRasterXYZ2 <- NULL
}
}
else {
originRasterXYZ2 <- originRasterXYZ; targetRasterXYZ2 <- targetRasterXYZ
}
if (any(isProb)) {
if (nTAss == 0) {
targetAssignment2 <- NULL
}
else {
if (nTAss < nTotal) {
dummyAss <- array(1/dimTAss[2], c(1, dimTAss[2]))
targetAssignment2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isProb[i] && (captured[i] != "target" || !capturePoint[i])) {
place <- place + 1
targetAssignment2 <- rbind(targetAssignment2,
targetAssignment[place, ])
}
else {
if (captured[i] == "target" && capturePoint[i]){
ass <- sf::st_intersects(x = targetPoints[i,], y = targetSites,
sparse = FALSE)
if (sum(ass)<1){
warning("Not all target capture locations are within targetSites. Assigning to closest site.\n",
"Affects animal: ", i)
ass <- matrix(0, 1, dimTAss[2])
ass[sf::st_nearest_feature(x = targetPoints[i,],
y = targetSites)] <- 1
}
else if (sum(ass) > 1){
warning("Overlapping targetSites may cause issues\n")
ass0 <- sf::st_intersects(x = targetPoints[i,], y = targetSites,
sparse = TRUE)
ass <- matrix(0, 1, dimTAss[2])
ass[ass0[1]] <- 1
}
targetAssignment2 <- rbind(targetAssignment2, ass)
}

else
targetAssignment2 <- rbind(targetAssignment2, dummyAss)
}
}
}
else
targetAssignment2 <- targetAssignment
}
if (is.null(dimOAss)) {
originAssignment2 <- NULL
}
else {
if (dimOAss[1] < nTotal) {
dummyAss <- array(1/dimOAss[2], c(1, dimOAss[2]))
originAssignment2 <- NULL
place <- 0
for (i in 1:nTotal) {
if (isProb[i] && (captured[i] != "origin" || !capturePoint[i])) {
place <- place + 1
originAssignment2 <- rbind(originAssignment2,
originAssignment[place, ])
}
else {
if (captured[i] == "origin" && capturePoint[i]){
ass <- sf::st_intersects(x = originPoints[i,], y = originSites,
sparse = FALSE)
if (sum(ass)<1){
warning("Not all origin capture locations are within originSites. Assigning to closest site.\n",
"Affects animal: ", i)
ass <- matrix(0, 1, dimOAss[2])
ass[sf::st_nearest_feature(x = originPoints[i,],
y = originSites)] <- 1

}
else if (sum(ass) > 1){
warning("Overlapping originSites may cause issues\n")
ass0 <- sf::st_intersects(x = originPoints[i,], y = originSites,
sparse = TRUE)
ass <- matrix(0, 1, dimOAss[2])
ass[ass0[1]] <- 1
}
originAssignment2 <- rbind(originAssignment2, ass)
}
else
originAssignment2 <- rbind(originAssignment2,
dummyAss)
}

}
}
else
originAssignment2 <- originAssignment
}
}
else {
originAssignment2 <- originAssignment
targetAssignment2 <- targetAssignment
}
return(list(originPoints = originPoints2, targetPoints = targetPoints2,
originAssignment = originAssignment2,
targetAssignment = targetAssignment2,
isGL = isGL, isTelemetry = isTelemetry, isRaster = isRaster,
isProb = isProb, captured = captured,
# originRaster = originRaster2,
originRasterXYZ = originRasterXYZ2,
originSingleCell = originSingleCell,
# targetRaster = targetRaster2,
targetRasterXYZ = targetRasterXYZ2,
targetSingleCell = targetSingleCell))
}

# Return randomly sampled points from sites, based on assignment
randomPoints <- function(sites, assignment, geomName = "geometry") {
nSites <- nrow(sites)
assignmentSum <- c(table(factor(assignment, levels = 1:nSites)))
points <- sf::st_sample(sites, assignmentSum)
points <- sf::st_as_sf(points)
points <- points[rank(assignment, ties.method = "first"), ]
names(points)[1] <- geomName
sf::st_geometry(points) <- geomName
return(points)
}

# Summary functions
summarize3D <- function(array3D, nSites1, nSites2, names1,
names2, point = NULL, alpha = 0.05) {
meanA <- apply(array3D, 2:3, mean)
medianA <- apply(array3D, 2:3, median)
seA <- apply(array3D, 2:3, sd)
simpleCIA <- apply(array3D, 2:3, quantile, probs = c(alpha/2, 1-alpha/2),
na.rm=TRUE)
dimnames(simpleCIA)[[1]] <- c("lower", "upper")
matrixA <- array(c(array3D), c(dim(array3D)[1], nSites1 * nSites2),
list(NULL, paste(rep(names1, nSites2),
rep(names2, each = nSites1),
sep = "#")))
mcmcA <- coda::as.mcmc(matrixA)
hpdCI <- coda::HPDinterval(mcmcA, 1-alpha)
hpdCI <- array(hpdCI, c(nSites1, nSites2, 2),
list(names1, names2, c("lower", "upper")))
hpdCI <- aperm(hpdCI, c(3, 1, 2))
bcCIA <- array(NA, dim = c(2, nSites1, nSites2),
dimnames = list(c("lower", "upper"), names1, names2))
for (i in 1:nSites1) {
for (j in 1:nSites2) {
psi.z0 <- qnorm(sum(array3D[, i, j] < meanA[i, j], na.rm = TRUE) /
length(which(!is.na(array3D[, i, j]))))
bcCIA[ , i, j] <- quantile(array3D[, i, j],
pnorm(2 * psi.z0 + qnorm(c(alpha/2, 1-alpha/2))),
na.rm=TRUE, names = FALSE)
}
}
return(list(sample = array3D, mean = meanA, se = seA,
simpleCI = simpleCIA, bcCI = bcCIA, hpdCI = hpdCI,
median = medianA, point = point))
}

summarizeAbund <- function(abund, nSites1, names1, pointAbund = NULL,
alpha = 0.05) {
meanAbund <- apply(abund, 2, mean, na.rm=TRUE)
medianAbund <- apply(abund, 2, median, na.rm=TRUE)
seAbund <- apply(abund, 2, sd, na.rm=TRUE)
# Calculate confidence intervals using quantiles of sampled MC
simpleCIAbund <- apply(abund, 2, quantile,
probs = c(alpha/2, 1-alpha/2),
na.rm=TRUE)
dimnames(simpleCIAbund)[[1]] <- c("lower", "upper")
bcCIAbund <- array(NA, dim = c(2, nSites1),
dimnames = list(c("lower", "upper"), names1))
for (i in 1:nSites1) {
abund.z0 <- qnorm(sum(abund[, i] < meanAbund[i], na.rm = TRUE) /
length(which(!is.na(abund[, i]))))
bcCIAbund[ , i] <- quantile(abund[, i],
pnorm(2 * abund.z0 + qnorm(c(alpha/2, 1-alpha/2))),
na.rm=TRUE, names = FALSE)
}
return(list(sample = abund, mean = meanAbund, se = seAbund,
simpleCI = simpleCIAbund, bcCI = bcCIAbund, median = medianAbund,
point = pointAbund))
}

# Currently called through reverseTransition (in calcConnectivity.R)
reverseEstTransition <- function(psi, originRelAbund, pi = NULL,
originSites=NULL, targetSites=NULL,
originNames = NULL, targetNames = NULL,
nSamples = 1000, row0 = 0, alpha = 0.05) {
if (is.null(psi)) {
if (is.array(pi)) {
if (length(dim(pi))!=3)
stop('Pi should either be 2-(for fixed migratory combination probabilities) or 3-dimensional array')
piIn <- pi
}
else if (inherits(pi, "estPi")) {
if (is.null(originNames)) {
originNames <- pi\$input\$originNames
}
if (is.null(targetNames)) {
targetNames <- pi\$input\$targetNames
}
piIn <- pi
pi <- pi\$pi\$sample
}
nOriginSites <- dim(pi)[2]
nTargetSites <- dim(pi)[3]
if (is.null(originNames)) {
originNames <- dimnames(pi)[[2]]
}
if (is.null(targetNames)) {
targetNames <- dimnames(pi)[[3]]
}
piBase <- apply(pi, 2:3, mean)
if (dim(pi)[1]>=nSamples)
piSamples <- round(seq(from = 1, to = dim(pi)[1],
length.out = nSamples))
else
piSamples <- sample.int(dim(pi)[1], nSamples, replace = TRUE)
pointRev <- reverseTransition(pi = piBase)
pointGamma <- pointRev\$gamma
pointTargetAbund <- pointRev\$targetRelAbund
pointPsi <- pointRev\$psi
pointOriginAbund <- pointRev\$originRelAbund
gamma <- array(NA, c(nSamples, nTargetSites, nOriginSites),
dimnames = list(NULL, targetNames, originNames))
targetRelAbund <- array(NA, c(nSamples, nTargetSites),
dimnames = list(NULL, targetNames))
psi <- array(NA, c(nSamples, nOriginSites, nTargetSites),
dimnames = list(NULL, originNames, targetNames))
originRelAbund <- array(NA, c(nSamples, nOriginSites),
dimnames = list(NULL, originNames))
for (i in 1:nSamples) {
piNew <- pi[piSamples[i],,]
# Reverse for new pi
sampleRev <- reverseTransition(pi = piNew)
gamma[i, , ] <- sampleRev\$gamma
targetRelAbund[i, ] <- sampleRev\$targetRelAbund
psi[i, , ] <- sampleRev\$psi
originRelAbund[i, ] <- sampleRev\$originRelAbund
}
gammaReturn <- summarize3D(gamma, nTargetSites, nOriginSites, targetNames,
originNames, pointGamma, alpha)
targetAbundReturn <- summarizeAbund(targetRelAbund, nTargetSites,
targetNames, pointTargetAbund, alpha)
psiReturn <- summarize3D(psi, nOriginSites, nTargetSites, originNames,
targetNames, pointPsi, alpha)
originAbundReturn <- summarizeAbund(originRelAbund, nOriginSites,
originNames, pointOriginAbund, alpha)
rev <- list(gamma = gammaReturn,
targetRelAbund = targetAbundReturn,
psi = psiReturn, originRekAbund = originAbundReturn,
input = list(pi = piIn,
originSites = originSites, targetSites = targetSites,
originNames = originNames, targetNames = targetNames,
nSamples = nSamples, row0 = row0, alpha = alpha))
class(rev) <- c("estGamma", "estTargetRelAbund", "estPsi",
"estOriginRelAbund", "estMigConnectivity")
return(rev)
}
if (is.matrix(psi)) {
psiFixed <- TRUE
psiVCV <- NULL
nOriginSites <- dim(psi)[1]
nTargetSites <- dim(psi)[2]
if (is.null(originNames)) {
originNames <- dimnames(psi)[[1]]
}
if (is.null(targetNames)) {
targetNames <- dimnames(psi)[[2]]
}
psiBase <- psi
psiIn <- psi
}
else if (inherits(psi, "mark")) {
psiFixed <- FALSE
if (!is.numeric(originSites) || !is.numeric(targetSites))
stop('Must specify which RMark Psi parameters represent origin and target sites')
nOriginSites <- length(originSites)
nTargetSites <- length(targetSites)
psiVCV <- psi\$results\$beta.vcv
psiBase <- RMark::TransitionMatrix(RMark::get.real(psi, "Psi",
se=TRUE))[originSites,
targetSites]
if (any(diag(psi\$results\$beta.vcv) < 0))
stop("Can't sample model, negative beta variances")
psiIn <- psi
}
else if (is.array(psi)) {
if (length(dim(psi))!=3)
stop('Psi should either be 2-(for fixed transition probabilities) or 3-dimensional array')
psiFixed <- FALSE
nOriginSites <- dim(psi)[2]
nTargetSites <- dim(psi)[3]
if (is.null(originNames)) {
originNames <- dimnames(psi)[[2]]
}
if (is.null(targetNames)) {
targetNames <- dimnames(psi)[[3]]
}
psiBase <- apply(psi, 2:3, mean)
psiVCV <- NULL
if (dim(psi)[1]>=nSamples)
psiSamples <- round(seq(from = 1, to = dim(psi)[1],
length.out = nSamples))
else
psiSamples <- sample.int(dim(psi)[1], nSamples, replace = TRUE)
psiIn <- psi
}
else if (inherits(psi, "estPsi") || inherits(psi, "estMC")) {
psiFixed <- FALSE
if (is.null(originNames)) {
originNames <- psi\$input\$originNames
}
if (is.null(targetNames)) {
targetNames <- psi\$input\$targetNames
}
psiIn <- psi
psi <- psi\$psi\$sample
nOriginSites <- dim(psi)[2]
nTargetSites <- dim(psi)[3]
psiBase <- apply(psi, 2:3, mean)
psiVCV <- NULL
if (dim(psi)[1]>=nSamples)
psiSamples <- round(seq(from = 1, to = dim(psi)[1],
length.out = nSamples))
else
psiSamples <- sample.int(dim(psi)[1], nSamples, replace = TRUE)
}
originRelAbundIn <- originRelAbund
if (coda::is.mcmc(originRelAbund) || coda::is.mcmc.list(originRelAbund)) {
originRelAbund <- as.matrix(originRelAbund)
}
if (is.matrix(originRelAbund) && dim(originRelAbund)[1]>1) {
abundFixed <- FALSE
if (dim(originRelAbund)[2]>nOriginSites)
abundParams <- paste('relN[', 1:nOriginSites, ']', sep='')
else if (dim(originRelAbund)[2]==nOriginSites)
abundParams <- 1:nOriginSites
else
stop('Number of origin sites must be constant between psi and abundance')
if (dim(originRelAbund)[1] >= nSamples)
abundRows <- round(seq(from = row0 + 1, to = dim(originRelAbund)[1],
length.out = nSamples))
else
abundRows <- sample((row0 + 1):dim(originRelAbund)[1], nSamples,
replace = TRUE)
originRelAbund <- as.matrix(originRelAbund[abundRows, abundParams])
abundBase <- colMeans(originRelAbund)
}
else {
abundFixed <- TRUE
if (length(originRelAbund)!=nOriginSites)
stop('Number of origin sites must be constant between psi and abundance')
abundBase <- originRelAbund
}
pointRev <- reversePsiRelAbund(psi = psiBase, originRelAbund = abundBase)
pointGamma <- pointRev\$gamma
pointAbund <- pointRev\$targetRelAbund
pointPi <- pointRev\$pi
gamma <- array(NA, c(nSamples, nTargetSites, nOriginSites),
dimnames = list(NULL, targetNames, originNames))
targetRelAbund <- array(NA, c(nSamples, nTargetSites),
dimnames = list(NULL, targetNames))
pi <- array(NA, c(nSamples, nOriginSites, nTargetSites),
dimnames = list(NULL, originNames, targetNames))
for (i in 1:nSamples) {
# Generate random transition probability matrices
if (psiFixed)
psiNew <- psiBase
else if (is.null(psiVCV))
psiNew <- psi[psiSamples[i],,]
else
psiNew <- makePsiRand(psi, originSites, targetSites)
# Point estimates of breeding densities
if (abundFixed)
abundNew <- abundBase
else
abundNew <- originRelAbund[i, abundParams]
# Reverse for new psis and relative breeding densities
sampleRev <- reversePsiRelAbund(psi = psiNew, originRelAbund = abundNew)
gamma[i, , ] <- sampleRev\$gamma
targetRelAbund[i, ] <- sampleRev\$targetRelAbund
pi[i, , ] <- sampleRev\$pi
}
gammaReturn <- summarize3D(gamma, nTargetSites, nOriginSites, targetNames,
originNames, pointGamma, alpha)
targetAbundReturn <- summarizeAbund(targetRelAbund, nTargetSites,
targetNames, pointAbund, alpha)
piReturn <- summarize3D(pi, nOriginSites, nTargetSites, originNames,
targetNames, pointPi, alpha)
rev <- list(gamma = gammaReturn,
targetRelAbund = targetAbundReturn,
pi = piReturn,
input = list(psi = psiIn, originRelAbund = originRelAbundIn,
originSites = originSites, targetSites = targetSites,
originNames = originNames, targetNames = targetNames,
nSamples = nSamples, row0 = row0, alpha = alpha))
class(rev) <- c("estGamma", "estTargetRelAbund", "estPi", "estMigConnectivity")
return(rev)

}

assignmentToMatrix <- function(assignment, nSites) {
nInds <- length(assignment)
assignMatrix <- array(0, c(nInds, nSites))
for (i in 1:nInds)
if (!is.na(assignment[i]))
assignMatrix[i, assignment[i]] <- 1
return(assignMatrix)
}

rasterToProb <- function(originSites = NULL, targetSites = NULL,
originAssignment = NULL, targetAssignment = NULL,
isRaster = FALSE, isProb = FALSE,
captured = "origin",
originRaster = NULL, targetRaster = NULL) {
failTarget <- failOrigin <- FALSE
nInds <- length(captured)
if (is.null(originSites))
failOrigin <- TRUE
else
nOriginSites <- nrow(originSites)
if (is.null(targetSites))
failTarget <- TRUE
else
nTargetSites <- nrow(targetSites)
inds <- isRaster & !isProb
capTarget <- captured=="target"
capOrigin <- captured=="origin"
whichTarget <- which(inds & !capTarget)
whichOrigin <- which(inds & !capOrigin)
if (length(whichTarget>0) && !is.null(targetSites)) {
targetSites <- sf::st_transform(targetSites, terra::crs(targetRaster))
targetSum <- terra::extract(targetRaster,
terra::vect(targetSites), fun = sum,
na.rm = TRUE, ID = FALSE)
targetSum <- t(targetSum[, whichTarget, drop = FALSE])
targetSum <- sweep(targetSum, 1, rowSums(targetSum), "/")
if (is.null(targetAssignment)){
targetAssignment <- array(0, c(nInds, nTargetSites))
}
else if (!is.array(targetAssignment) || length(dim(targetAssignment))<2)
targetAssignment <- assignmentToMatrix(targetAssignment, nTargetSites)
targetAssignment[whichTarget, ] <- targetSum
isRaster[whichTarget] <- FALSE
isProb[whichTarget] <- TRUE
}
if (length(whichOrigin>0) && !is.null(originSites)) {
originSites <- sf::st_transform(originSites, terra::crs(originRaster))
originSum <- terra::extract(originRaster,
terra::vect(originSites), fun = sum,
na.rm = TRUE, ID = FALSE)
originSum <- t(originSum[, whichOrigin, drop = FALSE])
originSum <- sweep(originSum, 1, rowSums(originSum), "/")
if (is.null(originAssignment)){
originAssignment <- array(0, c(nInds, nOriginSites))
}
else if (!is.array(originAssignment) || length(dim(originAssignment))<2)
originAssignment <- assignmentToMatrix(originAssignment, nOriginSites)
originAssignment[whichOrigin, ] <- originSum
isRaster[whichOrigin] <- FALSE
isProb[whichOrigin] <- TRUE
}
return(list(originAssignment = originAssignment,
targetAssignment = targetAssignment,
failOrigin = failOrigin, failTarget = failTarget,
isRaster = isRaster, isProb = isProb))
}

assignRasterStats <- function(theRaster) {
if (is.null(theRaster)){
PointsAssigned <- FALSE
SingleCell <- NULL
RasterXYZ <- NULL
RasterXYZcrs <- NULL
outRaster <- NULL
}
else {
if (inherits(theRaster, "isoAssign")) {
if (!inherits(theRaster\$probassign, "SpatRaster"))
theRaster\$probassign <- terra::rast(theRaster\$probassign)
PointsAssigned <- !is.null(theRaster\$SingleCell) &&
!is.null(dim(theRaster\$SingleCell))
RasterXYZ <- terra::as.data.frame(theRaster\$probassign, xy = TRUE)
RasterXYZcrs <- terra::crs(theRaster\$probassign)
SingleCell <- theRaster\$SingleCell
outRaster <- theRaster\$probassign
}
else {
if (!inherits(theRaster, "SpatRaster")) {
theRaster <- terra::rast(theRaster)
}
if (is.na(terra::crs(theRaster))){
stop("Please provide a projection (crs) for each raster\n")
}
RasterXYZ <- terra::as.data.frame(theRaster, xy = TRUE)
RasterXYZcrs <- terra::crs(theRaster)
SingleCell <- NULL
PointsAssigned <- FALSE
outRaster <- theRaster
}
}
return(list(Raster = outRaster, PointsAssigned = PointsAssigned,
SingleCell = SingleCell, RasterXYZ = RasterXYZ,
RasterXYZcrs = RasterXYZcrs))
}

# function to make an odds ratio (likely vs unlikely) assignment
oddsFun <- function(x, odds = odds){
predict(smooth.spline(x = cumsum(sort(x)),
sort(x),
spar = 0.1),(1-odds))\$y
}

propSpatRaster <- function(assignments) {
if (terra::nlyr(assignments)<2){
assign2prob  <- assignments /
unlist(c(terra::global(assignments, fun = "sum", na.rm = T)))
return(assign2prob)
}
summy <- terra::global(assignments, "sum", na.rm = TRUE)
test <- lapply(summy[,1],
function(x) terra::rast(nrows = terra::nrow(assignments),
ncols = terra::ncol(assignments),
xmin = terra::xmin(assignments),
xmax = terra::xmax(assignments),
ymin = terra::ymin(assignments),
ymax = terra::ymax(assignments),
vals = x,
crs = sf::st_crs(4326)\$wkt))
test <- terra::rast(test)
terra::crs(test) <- terra::crs(assignments)
assign2prob <- assignments / test
terra::crs(assign2prob) <- sf::st_crs(4326)\$wkt
return(assign2prob)
}
```

## Try the MigConnectivity package in your browser

Any scripts or data that you put into this service are public.

MigConnectivity documentation built on May 29, 2024, 8:37 a.m.