if (FALSE) {
# load data
sim_pu_raster <- get_sim_pu_raster()
sim_features <- get_sim_features()
sim_zones_pu_raster <- get_sim_zones_pu_raster()
sim_zones_features <- get_sim_zones_features()
# create minimal problem
p1 <-
problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.3) %>%
add_binary_decisions() %>%
add_default_solver(verbose = FALSE)
# create problem with contiguity constraints
p2 <- p1 %>% add_contiguity_constraints()
# create problem with constraints to represent features in contiguous
# units
p3 <- p1 %>% add_feature_contiguity_constraints()
# create problem with constraints to represent features in contiguous
# units that contain highly suitable habitat values
# (specifically in the top 5th percentile)
cm4 <- lapply(seq_len(terra::nlyr(sim_features)), function(i) {
# create connectivity matrix using the i'th feature's habitat data
m <- connectivity_matrix(sim_pu_raster, sim_features[[i]])
# convert matrix to 0/1 values denoting values in top 5th percentile
m <- round(m > quantile(as.vector(m), 1 - 0.05, names = FALSE))
# remove 0s from the sparse matrix
m <- Matrix::drop0(m)
# return matrix
m
})
p4 <- p1 %>% add_feature_contiguity_constraints(data = cm4)
# solve problems
s1 <- c(solve(p1), solve(p2), solve(p3), solve(p4))
names(s1) <- c(
"basic solution", "contiguity constraints",
"feature contiguity constraints",
"feature contiguity constraints with data"
)
# plot solutions
plot(s1, axes = FALSE)
# create minimal problem with multiple zones, and limit the solver to
# 30 seconds to obtain solutions in a feasible period of time
p5 <-
problem(sim_zones_pu_raster, sim_zones_features) %>%
add_min_set_objective() %>%
add_relative_targets(matrix(0.1, ncol = 3, nrow = 5)) %>%
add_binary_decisions() %>%
add_default_solver(time_limit = 30, verbose = FALSE)
# create problem with contiguity constraints that specify that the
# planning units used to conserve each feature in different management
# zones must form separate contiguous units
p6 <- p5 %>% add_feature_contiguity_constraints(diag(3))
# create problem with contiguity constraints that specify that the
# planning units used to conserve each feature must form a single
# contiguous unit if the planning units are allocated to zones 1 and 2
# and do not need to form a single contiguous unit if they are allocated
# to zone 3
zm7 <- matrix(0, ncol = 3, nrow = 3)
zm7[seq_len(2), seq_len(2)] <- 1
print(zm7)
p7 <- p5 %>% add_feature_contiguity_constraints(zm7)
# create problem with contiguity constraints that specify that all of
# the planning units in all three of the zones must conserve first feature
# in a single contiguous unit but the planning units used to conserve the
# remaining features do not need to be contiguous in any way
zm8 <- lapply(
seq_len(number_of_features(sim_zones_features)),
function(i) matrix(ifelse(i == 1, 1, 0), ncol = 3, nrow = 3)
)
print(zm8)
p8 <- p5 %>% add_feature_contiguity_constraints(zm8)
# solve problems
s2 <- lapply(list(p5, p6, p7, p8), solve)
s2 <- terra::rast(lapply(s2, category_layer))
names(s2) <- c("p5", "p6", "p7", "p8")
# plot solutions
plot(s2, axes = FALSE)
}
Run the code above in your browser using DataLab