diff --git a/.Rbuildignore b/.Rbuildignore index c0537c2..2850a5b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,3 +7,6 @@ ^docs$ ^pkgdown$ ^\.github$ +^doc$ +^Meta$ +figure$ diff --git a/.gitignore b/.gitignore index f47ffab..e77a603 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ .httr-oauth .DS_Store inst/doc +/doc/ +/Meta/ diff --git a/DESCRIPTION b/DESCRIPTION index 72da618..0cecf9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: patchwise Title: Create Patches of Feature Groups for Use in Spatial Prioritizations -Version: 0.1.0 +Version: 0.2.0 Authors@R: c( person(given = "Echelle S.", family = "Burns", email = "echelle_burns@ucsb.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3902-1410")), person(given = "Jason", family = "Flower", email = "jflower@ucsb.edu", role = c("aut"), comment = c(ORCID = "0000-0002-6731-8182")) @@ -20,18 +20,13 @@ Imports: tidyselect, rlang Roxygen: list(markdown = TRUE) -Remotes: - github::emlab-ucsb/spatialgridr, - github::emlab-ucsb/oceandatr -RoxygenNote: 7.3.1 Suggests: knitr, rmarkdown, codetools, testthat (>= 3.0.0), - oceandatr, - spatialgridr, prioritizr Config/testthat/edition: 3 VignetteBuilder: knitr URL: https://emlab-ucsb.github.io/patchwise/ +Config/roxygen2/version: 8.1.0 diff --git a/R/constraints_targets.R b/R/constraints_targets.R index 37cc4e2..c965ab2 100644 --- a/R/constraints_targets.R +++ b/R/constraints_targets.R @@ -1,41 +1,42 @@ -#' Create a dataframe for manual targets for all layers, including the features and the constraints +#' Create a dataframe for manual targets for all layers, including the features +#' and the constraints #' -#' @description This function creates a dataframe with the comibination of feature targets, generated by `features_targets()` and the targets for constraints to be used in prioritizr +#' @description This function creates a dataframe with the combination of +#' feature targets, generated by `features_targets()` and the targets for +#' constraints to be used in prioritizr #' -#' @param feature_targets a dataframe of feature targets generated by `features_targets()` -#' @param patch_df a dataframe generated by `create_patch_df()` that includes constraints for each patch and grid cell combination +#' @param feature_targets a dataframe of feature targets generated by +#' `features_targets()` +#' @param patch_df a dataframe generated by `create_patch_df()` that includes +#' constraints for each patch and grid cell combination #' #' @return A dataframe to be used to specify all manual targets for prioritizr #' @export #' #' @examples -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) -#'# Create patch dataframe -#'patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, -#' patches = patches_raster, costs = cost_raster) -#'# Create boundary matrix for prioritizr -#'boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, -#' patches = patches_raster, patch_df = patches_raster_df) -#'# Create target features - using just 20% for every feature -#'features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), -#' features = features_raster, pre_patches = seamounts_raster) -#'# Create constraint targets -#'constraint_targets <- constraints_targets(feature_targets = features_targets, -#' patch_df = patches_raster_df) +#'# Import some planning data +#' +#'#import planning units/ cost raster +#'pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + +#'#import species distributions - feature data for planning +#'species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) +#' +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#'# Create seamount patches as multi-layer raster +#'patches_raster <- create_patches(seamounts) +#' +#'patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +#'patches = patches_raster, costs = pu_raster) +#' +#'# Create targets for protection - use 20% for each feature (including 20% of whole seamounts) in this example +#'targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) +#' +#'# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +#'constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) + constraints_targets <- function(feature_targets, patch_df) { diff --git a/R/convert_solution.R b/R/convert_solution.R index 7da12c0..f8f353f 100644 --- a/R/convert_solution.R +++ b/R/convert_solution.R @@ -1,56 +1,61 @@ #' Converts prioritization solution into a more digestible output #' -#' @description This function converts the prioritization solution that uses patches to a raster or sf object (depending on the format of your spatial_grid) that clearly marks areas suggested for protection (1) and areas not suggested for protection (0) +#' @description This function converts the prioritization solution that uses +#' patches to a raster or sf object (depending on the format of your +#' spatial_grid) that clearly marks areas suggested for protection (1) and +#' areas not suggested for protection (0) #' -#' @param solution the solution that results from using `prioritizr::solve()` on `prioritizr::problem()` -#' @param patch_df a dataframe generated by `create_patch_df()` that includes constraints for each patch and grid cell combination -#' @param spatial_grid a raster or sf template with the desired resolution and coordinate reference system generated by `spatialgridr::get_grid()`; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object) +#' @param solution the solution that results from using `prioritizr::solve()` on +#' `prioritizr::problem()` +#' @param patch_df a dataframe generated by `create_patch_df()` that includes +#' constraints for each patch and grid cell combination +#' @param spatial_grid a raster or sf template with the desired resolution and +#' coordinate reference system; values in areas of interest are 1, while all +#' other values are NA (only required if feature is a sf object) #' -#' @return A raster or sf object that denotes areas suggested for protection (value = 1) +#' @return A raster or sf object that denotes areas suggested for protection +#' (value = 1) #' @export #' #' @examples #' \dontrun{ -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) -#'# Create patch dataframe -#'patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, -#' patches = patches_raster, costs = cost_raster) -#'# Create boundary matrix for prioritizr -#'boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, patches = patches_raster, -#' patch_df = patches_raster_df) -#'# Create target features - using just 20% for every feature -#'features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), -#' features = features_raster, pre_patches = seamounts_raster) -#'# Create constraint targets -#'constraint_targets <- constraints_targets(feature_targets = features_targets, -#' patch_df = patches_raster_df) -#'# Create prioritization problem -#'problem_raster <- prioritizr::problem(x = patches_raster_df, -#' features = constraint_targets$feature, cost_column = "cost") %>% -#' prioritizr::add_min_set_objective() %>% -#' prioritizr::add_manual_targets(constraint_targets) %>% -#' prioritizr::add_binary_decisions() %>% -#' prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix) %>% -#' prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) -#'# Solve problem -#'solution <- solve(problem_raster) -#'# Convert to a more digestible format -#'suggested_protection <- convert_solution(solution = solution, patch_df = patches_raster_df, -#' spatial_grid = planning_raster) +#'#import planning units/ cost raster +#'pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + +#'#import species distributions - feature data for planning +#'species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) +#' +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#'# Create seamount patches as multi-layer raster +#'patches_raster <- create_patches(seamounts) +#' +#'patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +#'patches = patches_raster, costs = pu_raster) +#' +#'# Create targets for protection - use 20% for each feature (including 20% of whole seamounts) in this example +#'targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) +#' +#'# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +#'constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) +#' +#'library(prioritizr) +#' +#'# Create the prioritization problem +#'problem_patches <- problem(x = patches_df, features = constraints_df$feature, cost_column = "cost") |> +#' add_min_set_objective() |> +#' add_manual_targets(constraints_df) |> +#' add_binary_decisions() |> +#' add_default_solver() +#' +#'# Solve the problem +#'solution_patches_tbl <- solve(problem_patches) +#' +#'# Convert the solution into a raster using `convert_solution()` +#'sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df, spatial_grid = pu_raster) +#' +#'terra::plot(sol_rast_patches) #' } convert_solution <- function(solution, patch_df, spatial_grid) { @@ -61,7 +66,7 @@ convert_solution <- function(solution, patch_df, spatial_grid) { if(!check_df(patch_df)) { stop("patch_df must be a dataframe object")} - planning_unit_id <- unique(unlist(patch_df$id[which(solution$solution_1 > 0.5)])) + planning_unit_id <- unique(unlist(patch_df$idx[which(solution$solution_1 > 0.5)])) if(class(spatial_grid)[1] %in% c("RasterLayer", "SpatRaster")){ solution <- spatial_grid*0 diff --git a/R/create_boundary_matrix.R b/R/create_boundary_matrix.R index 45a2c11..f18f24b 100644 --- a/R/create_boundary_matrix.R +++ b/R/create_boundary_matrix.R @@ -1,36 +1,42 @@ #' Create boundary matrix for prioritizr if you are using patches #' -#' @description This function creates a boundary matrix to plug into prioritizr +#' @description This function creates a boundary matrix which is required if you +#' want to use a boundary penalty in prioritizr #' -#' @param spatial_grid a raster or sf template with the desired resolution and coordinate reference system generated by `spatialgridr::get_grid()`; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object) -#' @param patches a raster or sf object generated by `create_patches()`; a single feature split into spatially distinct groups; each layer of the raster or each column of the sf object identifies the location of each patch -#' @param patch_df a dataframe generated by `create_patch_df()` that includes constraints for each patch and grid cell combination +#' @param spatial_grid a raster or sf template with the desired resolution and +#' coordinate reference system; values in areas of interest are 1, while all +#' other values are NA (only required if feature is a sf object) +#' @param patches a raster or sf object generated by `create_patches()`; a +#' single feature split into spatially distinct groups; each layer of the +#' raster or each column of the sf object identifies the location of each +#' patch +#' @param patch_df a dataframe generated by `create_patch_df()` that includes +#' constraints for each patch and grid cell combination #' -#' @return A boundary matrix to plug into prioritizr +#' @return A sparse matrix #' @export #' #' @examples -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) -#'# Create patch dataframe -#'patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, -#' patches = patches_raster, costs = cost_raster) -#'# Create boundary matrix for prioritizr -#'boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, -#' patches = patches_raster, patch_df = patches_raster_df) +#'# Import some planning data +#' +#'#import planning units/ cost raster +#'pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + +#'#import species distributions - feature data for planning +#'species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) +#' +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#'# Create seamount patches as multi-layer raster +#'patches_raster <- create_patches(seamounts) +#' +#'patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +#'patches = patches_raster, costs = pu_raster) +#' +#'#Create boundary matrix if you want to use a boundary penatly function in prioritizr +#'boundary_matrix <- create_boundary_matrix(spatial_grid = pu_raster, patches = patches_raster, patch_df = patches_df) + create_boundary_matrix <- function(spatial_grid, patches, patch_df){ diff --git a/R/create_patch_df.R b/R/create_patch_df.R index 3dce719..f9848c6 100644 --- a/R/create_patch_df.R +++ b/R/create_patch_df.R @@ -1,37 +1,52 @@ -#' Create a dataframe that counts each patch independently for the prioritizr input +#' Create a dataframe that counts each patch independently for the prioritizr +#' input #' -#' @description This function takes all of the inputs of prioritizr and creates a dataframe precursor for prioritizr. This specific function includes the use of "patches" which have been already generated by `create_patches()`. "Patches" are spatial groups of a particular feature, such as seamounts. Prioritizr will protect entire patches to meet the protection goal for that feature. The output of this function should be passed to `create_boundary_matrix()` before being passed to prioritizr. +#' @description This function takes all of the inputs of prioritizr and creates +#' a dataframe precursor for prioritizr. This specific function includes the +#' use of "patches" which have been already generated by `create_patches()`. +#' "Patches" are spatial groups of a particular feature, such as seamounts. +#' Prioritizr will protect entire patches to meet the protection goal for that +#' feature. The output of this function should be passed to +#' `create_boundary_matrix()` before being passed to prioritizr. #' -#' @param spatial_grid a raster or sf template with the desired resolution and coordinate reference system generated by `spatialgridr::get_grid()`; values in areas of interest are 1, while all other values are NA -#' @param features a raster or sf object that includes all relevant features to be used in the prioritization; each layer of the raster or each column of the sf object identifies the location of each feature -#' @param patches a raster or sf object generated by `create_patches()`; a single feature split into spatially distinct groups; each layer of the raster or each column of the sf object identifies the location of each patch -#' @param costs a raster or sf object of the costs for protecting each cell in the planning grid -#' @param locked_out a raster or sf object for areas to be locked out (absolutely not protected) in the prioritization -#' @param locked_in a raster or sf object for areas to be locked in (absolutely protected) in the prioritization +#' @param spatial_grid a raster or sf template with the desired resolution and +#' coordinate reference system; values in areas of interest are 1, while all +#' other values are NA +#' @param features a raster or sf object that includes all relevant features to +#' be used in the prioritization; each layer of the raster or each column of +#' the sf object identifies the location of each feature +#' @param patches a raster or sf object generated by `create_patches()`; a +#' single feature split into spatially distinct groups; each layer of the +#' raster or each column of the sf object identifies the location of each +#' patch +#' @param costs a raster or sf object of the costs for protecting each cell in +#' the planning grid +#' @param locked_out a raster or sf object for areas to be locked out +#' (absolutely not protected) in the prioritization +#' @param locked_in a raster or sf object for areas to be locked in (absolutely +#' protected) in the prioritization #' #' @return A data frame to be used as an input for `create_boundary_matrix()` #' @importFrom rlang := #' @export #' #' @examples -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) -#'# Create patch dataframe -#'patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, -#' patches = patches_raster, costs = cost_raster) +#'# Import some planning data +#' +#'#import planning units/ cost raster +#'pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + +#'#import species distributions - feature data for planning +#'species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) +#' +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#'# Create seamount patches as multi-layer raster +#'patches_raster <- create_patches(seamounts) +#' +#'patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +#'patches = patches_raster, costs = pu_raster) create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locked_out = NULL, locked_in = NULL){ @@ -58,7 +73,7 @@ create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locke if (class(spatial_grid)[1] %in% c("RasterLayer", "SpatRaster")) { # Initialize - pu_grid_data <- tibble::tibble(id = as.list(seq_len(nrow(terra::as.data.frame(spatial_grid, na.rm = FALSE))))) %>% + pu_grid_data <- tibble::tibble(idx = as.list(seq_len(nrow(terra::as.data.frame(spatial_grid, na.rm = FALSE))))) %>% dplyr::bind_cols(tibble::tibble(terra::as.data.frame(features, na.rm = FALSE))) %>% dplyr::mutate(patch = 0) @@ -73,7 +88,7 @@ create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locke pu_sm_data <- lapply(names(patches), function(i) { curr_sm_pu <- tibble::tibble( - id = list(as.numeric(row.names(terra::as.data.frame(patches[[i]], na.rm = FALSE))[which(terra::as.data.frame(patches[[i]], na.rm = FALSE) > 0.5)]))) %>% + idx = list(as.numeric(row.names(terra::as.data.frame(patches[[i]], na.rm = FALSE))[which(terra::as.data.frame(patches[[i]], na.rm = FALSE) > 0.5)]))) %>% dplyr::bind_cols( terra::as.data.frame(features * patches[[i]], na.rm = FALSE) %>% stats::setNames(names(features)) %>% @@ -127,7 +142,7 @@ create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locke } # Initialize - pu_grid_data <- tibble::tibble(id = as.list(seq_len(nrow(spatial_grid)))) %>% + pu_grid_data <- tibble::tibble(idx = as.list(seq_len(nrow(spatial_grid)))) %>% dplyr::bind_cols(features) %>% dplyr::mutate(patch = 0) @@ -144,7 +159,7 @@ create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locke pu_sm_data <- lapply(names(patches), function(i) { curr_sm_pu <- tibble::tibble( - id = list(as.numeric(row.names(patches[which(patches[,i] > 0.5),])))) %>% + idx = list(as.numeric(row.names(patches[which(patches[,i] > 0.5),])))) %>% dplyr::bind_cols(data.frame(t(features %>% dplyr::mutate_all(., ~(.*patches[,i])) %>% colSums(., na.rm = T)))) @@ -182,14 +197,17 @@ create_patch_df <- function(spatial_grid, features, patches, costs = NULL, locke # planning units constraints <- data.frame(test = rep(0,nrow(pu_data))) + # add in id column for planning unit data + pu_data$id <- seq_len(nrow(pu_data)) + index <- 1 #Create a vector with a '1' for each combination of patch-level # and grid-cell level planning unit that overlaps with that patch. for (i in seq_len(nrow(pu_sm_data))) { - print(paste0("Processing patch ", i, " of ", nrow(pu_sm_data))) + message(paste0("Processing patch ", i, " of ", nrow(pu_sm_data))) - for (j in pu_sm_data$id[[i]]) { + for (j in pu_sm_data$idx[[i]]) { v <- rep(0, nrow(pu_data)) # initialize with zeros v[nrow(pu_grid_data) + i] <- 1 # specify patch-level planning unit v[as.numeric(j)] <- 1 # specify grid cell-level planning unit diff --git a/R/create_patches.R b/R/create_patches.R index 0fd40af..e48401b 100644 --- a/R/create_patches.R +++ b/R/create_patches.R @@ -1,29 +1,30 @@ -#' Create patches for a specific feature +#'Create patches for a specific feature #' -#' @description This function takes a feature and splits it into multiple "patches" of the feature. This is ideal for features that you want to protect x% of, but want to make sure that entire patches (as opposed to pieces of patches) are protected to meet that target. +#'@description This function takes a feature and splits it into multiple +#' "patches" of the feature. This is ideal for features that you want to +#' protect x% of, but want to make sure that entire patches (as opposed to +#' pieces of patches) are protected to meet that target. #' -#' @param feature a raster or sf object with the feature of interest present (value = 1) or absent (value = NA) -#' @param spatial_grid a raster or sf template with the desired resolution and coordinate reference system generated by `spatialgridr::get_grid()`; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object) +#'@param feature a raster or sf object with the feature of interest present +#' (value = 1) or absent (value = NA) +#'@param spatial_grid a raster or sf template with the desired resolution and +#' coordinate reference system generated; values in areas of interest are 1, +#' while all other values are NA (only required if feature is a sf object) #' -#' @return A raster or sf object with independent layers (raster)/columns (sf) designating the location of each patch -#' @export +#'@return A raster or sf object with independent layers (raster)/columns (sf) +#' designating the location of each patch +#'@export #' #' @examples -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) +#'# Import some planning data +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#'terra::plot(seamounts) +#' +#'# Create seamount patches as multi-layer raster +#'patches_raster <- create_patches(seamounts) +#'terra::plot(patches_raster) create_patches <- function(feature, spatial_grid = NULL) { diff --git a/R/features_targets.R b/R/features_targets.R index 965abd7..7f7b737 100644 --- a/R/features_targets.R +++ b/R/features_targets.R @@ -1,41 +1,39 @@ #' Create a dataframe for the relative and absolute targets for features #' -#' @description This function creates a dataframe of targets for each of the features to be used in the prioritization +#' @description This function creates a dataframe of targets for each of the +#' features to be used in the prioritization #' -#' @param targets a vector of targets for protection (range between 0 and 1); must be the same length as the number of features + the pre-patches variable -#' @param features a raster or sf object that includes all relevant features to be used in the prioritization; each layer of the raster or each column of the sf object identifies the location of each feature -#' @param pre_patches a raster or sf object that includes the feature that is to be split into patches (not the layer that is already split into patches) -#' @param locked_out a raster or sf object for areas to be locked out (absolutely not protected) in the prioritization -#' @param locked_in a raster or sf object for areas to be locked in (absolutely protected) in the prioritization +#' @param targets a vector of targets for protection (range between 0 and 1); +#' must be the same length as the number of features + the pre-patches +#' variable +#' @param features a raster or sf object that includes all relevant features to +#' be used in the prioritization; each layer of the raster or each column of +#' the sf object identifies the location of each feature +#' @param pre_patches a raster or sf object that includes the feature that is to +#' be split into patches (not the layer that is already split into patches) +#' @param locked_out a raster or sf object for areas to be locked out +#' (absolutely not protected) in the prioritization +#' @param locked_in a raster or sf object for areas to be locked in (absolutely +#' protected) in the prioritization #' -#' @return A data frame to be used to specify targets for features. Will need to be plugged into `constraint_targets()` before being implemented into prioritizr +#' @return A data frame to be used to specify targets for features. Will need to +#' be plugged into `constraint_targets()` before being implemented into +#' prioritizr #' @export #' #' @examples -#'# Start with a little housekeeping to get the data from oceandatr -#'# Choose area of interest (Bermuda EEZ) -#'area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -#'projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -#'# Create a planning grid -#'planning_raster <- spatialgridr::get_grid(area, projection = projection) -#'# Grab all relevant data -#'features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -#'# Separate seamount data - we want to protect entire patches -#'seamounts_raster <- features_raster[["seamounts"]] -#'features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -#'# Create a "cost" to protecting a cell - just a uniform cost for this example -#'cost_raster <- stats::setNames(planning_raster, "cost") -#'# Create patches from layer -#'patches_raster <- create_patches(seamounts_raster) -#'# Create patch dataframe -#'patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, -#' patches = patches_raster, costs = cost_raster) -#'# Create boundary matrix for prioritizr -#'boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, patches = patches_raster, -#' patch_df = patches_raster_df) -#'# Create target features - using just 20% for every feature -#'features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), -#' features = features_raster, pre_patches = seamounts_raster) +#'# Import some planning data +#' +#'#import species distributions - feature data for planning +#'species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) +#' +#'#import seamounts +#'seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) +#' +#' +#'# Create targets for protection - use 20% for each feature (including 20% of whole seamounts) in this example +#'targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) +#'head(targets_df) features_targets <- function(targets, features, pre_patches, locked_out = NULL, locked_in = NULL){ diff --git a/README.Rmd b/README.Rmd new file mode 100644 index 0000000..7baef07 --- /dev/null +++ b/README.Rmd @@ -0,0 +1,256 @@ +--- +output: github_document +--- + + + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "man/figures/README-", + out.width = "100%" +) + +devtools::load_all() +library(prioritizr) +library(terra) +library(sf) +``` + +# patchwise patchwise website + + + + +`patchwise` is intended to be used as a supplementary package to `prioritizr` for instances in which users wish to protect entire contiguous patches of features rather than portions of many features. For example, consider conservation planning for an area of ocean where seamounts are one of the biodiversity features that are targeted for protection. If the seamounts span multiple planning units and representation target of say 20% is used, portions of many seamounts could be protected, but it might be better to protect the entirety (contiguous patches) of a smaller number of seamounts. `patchwise` provides this option, ensuring that representation targets are met by representing whole features. + +## Installation + +You can install the development version of `patchwise` from [GitHub](https://github.com/) with: + +```{r eval=FALSE} +# install.packages("pak") +pak::pak("emlab-ucsb/patchwise") +``` + +## Example using raster data + +First load the libraries needed. Apart from `patchwise`, the `terra` package is used for raster data manipulation, and `prioritizr` is used for the spatial prioritization. + +```{r libraries, eval=FALSE} +library(patchwise) +library(prioritizr) +library(terra) +``` + +We will import some basic planning data to demonstrate how `patchwise` can be used. For this example we will use a a 40 x 40 raster planning grid which has a cost value of 1 for each cell, and the following features that will be targeted in the prioritization: + +- 4 random, binary (0 or 1) data layers, representing species distributions; for this example we are calling them fish +- 3 "patches" of contiguous features; for this example we are calling them seamounts + +```{r data-setup} +#import planning units/ cost raster +pu_raster <- rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + +#import fish distributions +fish_distributions <- rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- rast(system.file("extdata/seamounts.tif", package = "patchwise")) +``` + +Let's look at our fish distributions and our seamounts + +```{r visualize-data} +plot(c(fish_distributions, seamounts)) +``` + +Now we can use `patchwise` to do some pre-processing of the seamounts data so the seamounts can be prioritized as whole patches in the following prioritization. + +```{r patchwise-prep} +# Create seamount patches +patches_rast <- create_patches(seamounts) + +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_rast <- create_patch_df(spatial_grid = pu_raster, features = fish_distributions, patches = patches_rast, costs = pu_raster) +``` + +With that pre-processing done, we can now use `patchwise` to create protection targets for our features, including seamounts. In this example, we will use 20%, including 20% of whole seamounts. + +```{r targets-constraints} +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_rast <- features_targets(targets = rep(0.2, (nlyr(fish_distributions) + 1)), features = fish_distributions, pre_patches = seamounts) + +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) +``` + +With all the data preparation now done, we can run a prioritization using `prioritizr`: + +```{r prioritization, include=FALSE} +# Create the prioritization problem +problem_patches <- problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_rast) |> + add_binary_decisions() |> + add_default_solver() + +# Solve the prioritization +solution_patches_tbl <- solve(problem_patches) +``` + +The solution object is a tibble. To convert this into a raster object for plotting, we use the `convert_solution()` function from `patchwise` + +```{r convert-solution} +# Convert the solution into a raster using the patchwise function `convert_solution()` +sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df_rast, spatial_grid = pu_raster) |> + setNames("With patchwise") +``` + +For comparison, we will run a prioritization without using `patchwise` + +```{r prioritization-no-patches, include=FALSE} +features_rast <- c(fish_distributions, seamounts) + +problem_no_patches <- problem(x = pu_raster, features = features_rast) |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() + +# Solve the prioritization +solution_no_patches <- solve(problem_no_patches) |> + setNames("Without patchwise") +``` + +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). + +```{r prioritization-plot} +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(as.polygons(seamounts), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) +``` + +If you want to use a boundary penalty in the prioritization, we need to manually create a boundary matrix using the `patchwise` function `create_boundary_matrix()` + +```{r boundary-matrix} +boundary_matrix_rast <- create_boundary_matrix(spatial_grid = pu_raster, patches = patches_rast, patch_df = patches_df_rast) +``` + +We can now re-run the prioritization with a boundary penalty + +```{r prioritization-boundary, include=FALSE} +problem_boundary <- problem_patches |> + add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) + +solution_boundary_tbl <- solve(problem_boundary) + +solution_rast_boundary <- convert_solution(solution = solution_boundary_tbl, patch_df = patches_df_rast, spatial_grid = pu_raster) +``` + +We now get a solution with planning units more clustered together + +```{r prioritization-boundary-plot} +plot(solution_rast_boundary, plg = list(legend = c("Not selected", "Selected"))) +lines(as.polygons(seamounts), col = "red") +``` + +## Example with sf data + +The previous example used raster input data for the prioritization, but `patchwise` also handles `sf` data. + +```{r eval=FALSE} +library(sf) +``` + + +First we need to create suitable `sf` data inputs. We can do this by polygonizing the raster data: + +```{r} +#create the planning grid, which is also the same as the cost grid since we are using planning units all with cost = 1 +pu_sf <- as.polygons(pu_raster, aggregate = FALSE) |> + st_as_sf() + +features_sf <- as.polygons(fish_distributions, aggregate = FALSE) |> + st_as_sf() + +seamounts_sf <- as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + st_as_sf() + +#replace NAs with zeroes +seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 +``` + +Let's check our features and seamounts look ok: + +```{r} +plot(cbind(features_sf, st_drop_geometry(seamounts_sf))) +``` + +Now we can run the same `patchwise` functions as we did with raster data to prepare the data for prioritization + +```{r patchwise-prep-sf} + +# Create seamount patches +patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) + +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf) + +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf) + +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +``` + +With all the data preparation now done, we can run a prioritization using `prioritizr`: + +```{r prioritization_sf, include=FALSE} +# Create the prioritization problem +problem_patches_sf <- problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_sf) |> + add_binary_decisions() |> + add_default_solver() + +# Solve the prioritization +solution_patches_sf_tbl <- solve(problem_patches_sf) +``` + +The solution object is a tibble. To convert this into a raster object for plotting, we use the `convert_solution()` function from `patchwise` + +```{r convert-solution-sf} +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) +``` + +For comparison, we will run a prioritization without using `patchwise` + +```{r prioritization-no-patches-sf, include=FALSE} +sf_planning_data <- cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(pu_sf)) + +problem_no_patches_sf <- problem(x = sf_planning_data, + features = names(st_drop_geometry(sf_planning_data)[-6]), + cost_column = "cost") |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() + +# Solve the prioritization +solution_no_patches_sf <- solve(problem_no_patches_sf) +``` + +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). + +```{r prioritization-plot-sf} +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(as.polygons(seamounts, aggregate = TRUE), col = "red")) +``` diff --git a/README.md b/README.md index 312b841..459e705 100644 --- a/README.md +++ b/README.md @@ -1,214 +1,236 @@ + + + # patchwise patchwise website -`patchwise` is intended to be used as a supplementary package to `oceandatr` (and `spatialgridr`) for instances in which users wish to protect entire "chunks" of areas using `prioritizr`. -One example is when a user wishes to include seamounts as a feature to protect in `prioritizr` with a target of 20% protection. Instead of protecting a little bit of each seamount until the 20% is reached, `patchwise` makes it easy to ensure that entire seamounts are protected sequentially to meet the protection target. + -## Installation -You can install the development version of `patchwise` from GitHub with: + -``` -if (!require(devtools)) install.packages("devtools") -devtools::install_github("emlab-ucsb/patchwise") -``` +`patchwise` is intended to be used as a supplementary package to +`prioritizr` for instances in which users wish to protect entire +contiguous patches of features rather than portions of many features. +For example, consider conservation planning for an area of ocean where +seamounts are one of the biodiversity features that are targeted for +protection. If the seamounts span multiple planning units and +representation target of say 20% is used, portions of many seamounts +could be protected, but it might be better to protect the entirety +(contiguous patches) of a smaller number of seamounts. `patchwise` +provides this option, ensuring that representation targets are met by +representing whole features. + +## Installation -You can install `oceandatr` and `spatialgridr` from GitHub with: +You can install the development version of `patchwise` from +[GitHub](https://github.com/) with: +``` r +# install.packages("pak") +pak::pak("emlab-ucsb/patchwise") ``` -if (!require(devtools)) install.packages("devtools") -devtools::install_github("emlab-ucsb/oceandatr") -devtools::install_github("emlab-ucsb/spatialgridr") -``` - -## Examples of usage -### Using `raster` objects as inputs +## Example using raster data -Since this package is intended to be used in combination with `oceandatr`, there are several housekeeping steps that need to be completed first. +First load the libraries needed. Apart from `patchwise`, the `terra` +package is used for raster data manipulation, and `prioritizr` is used +for the spatial prioritization. -``` -# Load libraries +``` r library(patchwise) +library(prioritizr) +library(terra) +``` -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <- '+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' - -# Create a planning grid -planning_rast <- spatialgridr::get_grid(area, projection = projection) +We will import some basic planning data to demonstrate how `patchwise` +can be used. For this example we will use a a 40 x 40 raster planning +grid which has a cost value of 1 for each cell, and the following +features that will be targeted in the prioritization: -# Grab all relevant data -features_rast <- oceandatr::get_features(spatial_grid = planning_rast) +- 4 random, binary (0 or 1) data layers, representing species + distributions; for this example we are calling them fish +- 3 “patches” of contiguous features; for this example we are calling + them seamounts -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_rast <- stats::setNames(planning_rast, "cost") +``` r +#import planning units/ cost raster +pu_raster <- rast(system.file("extdata/pu_raster.tif", package = "patchwise")) -# Separate seamount data - we want to protect entire patches -seamounts_rast <- features_rast[["seamounts"]] -features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] +#import fish distributions +fish_distributions <- rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) -# Show what seamounts look like... -terra::plot(seamounts_rast) # there are 7 seamount areas (seamounts that are touching) +#import seamounts +seamounts <- rast(system.file("extdata/seamounts.tif", package = "patchwise")) ``` -![](https://github.com/echelleburns/patchwise/assets/40546424/7e717d95-4673-4dac-8f59-dbb4865398ea) +Let’s look at our fish distributions and our seamounts +``` r +plot(c(fish_distributions, seamounts)) ``` -# Create seamount patches - seamount areas that touch are considered the same patch -patches_rast <- patchwise::create_patches(seamounts_rast) -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast) + -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +Now we can use `patchwise` to do some pre-processing of the seamounts +data so the seamounts can be prioritized as whole patches in the +following prioritization. -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast) +``` r +# Create seamount patches +patches_rast <- create_patches(seamounts) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) - -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_rast <- create_patch_df(spatial_grid = pu_raster, features = fish_distributions, patches = patches_rast, costs = pu_raster) +#> Processing patch 1 of 3 +#> Processing patch 2 of 3 +#> Processing patch 3 of 3 +``` -# Solve the prioritization -solution_rast <- solve(problem_rast) +With that pre-processing done, we can now use `patchwise` to create +protection targets for our features, including seamounts. In this +example, we will use 20%, including 20% of whole seamounts. -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +``` r +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_rast <- features_targets(targets = rep(0.2, (nlyr(fish_distributions) + 1)), features = fish_distributions, pre_patches = seamounts) -# Show the results -terra::plot(result_rast) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) ``` -![](https://github.com/emlab-ucsb/patchwise/assets/40546424/82d88030-e915-4fde-a3ec-a99b62ff596d) - +With all the data preparation now done, we can run a prioritization +using `prioritizr`: -Areas in green were identified by `prioritizr` as areas worth protecting. We can see that entire seamounts were used to meet the target objective of protecting 20% of seamounts. We can compare this result to a `prioritizr` run that does not protect whole seamounts: +The solution object is a tibble. To convert this into a raster object +for plotting, we use the `convert_solution()` function from `patchwise` +``` r +# Convert the solution into a raster using the patchwise function `convert_solution()` +sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df_rast, spatial_grid = pu_raster) |> + setNames("With patchwise") ``` -# Grab all relevant data -features_rast_nopatch <- oceandatr::get_features(spatial_grid = planning_rast) -# Run the prioritization -problem_rast_nopatch <- prioritizr::problem(x = cost_rast, features = features_rast_nopatch) %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_relative_targets(rep(0.2, terra::nlyr(features_rast_nopatch))) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +For comparison, we will run a prioritization without using `patchwise` -# Solve the prioritization -solution_nopatch <- solve(problem_rast_nopatch) +We can now plot the solutions, and overlaying the outlines of the +seamounts (in red), we can see that one entire seamount is included in +the solution that used `patchwise`, whereas the solution without +`patchwise` selects a few planning units in each seamount. Note that in +the solution with `patchwise`, planning units that overlap seamounts +that are not entirely selected have also been selected to meet targets +for other features (fish distributions). -# Show the results -terra::plot(solution_nopatch) +``` r +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(as.polygons(seamounts), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) ``` -![](https://github.com/emlab-ucsb/patchwise/assets/40546424/0d27f6df-32e4-49ff-abbc-720b930c26d4) + -Only portions of seamount units are protected here. +If you want to use a boundary penalty in the prioritization, we need to +manually create a boundary matrix using the `patchwise` function +`create_boundary_matrix()` -### Using `sf` objects as inputs +``` r +boundary_matrix_rast <- create_boundary_matrix(spatial_grid = pu_raster, patches = patches_rast, patch_df = patches_df_rast) +``` + +We can now re-run the prioritization with a boundary penalty -Since this package is intended to be used in combination with `oceandatr`, there are several housekeeping steps that need to be completed first. +We now get a solution with planning units more clustered together +``` r +plot(solution_rast_boundary, plg = list(legend = c("Not selected", "Selected"))) +lines(as.polygons(seamounts), col = "red") ``` -# Load libraries -library(patchwise) -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <- '+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' + -# Create a planning grid -planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") +## Example with sf data -# Grab all relevant data -features_sf <- oceandatr::get_features(spatial_grid = planning_sf) +The previous example used raster input data for the prioritization, but +`patchwise` also handles `sf` data. -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) +``` r +library(sf) +``` -# Separate seamount data - we want to protect entire patches -seamounts_sf <- features_sf %>% - dplyr::select(seamounts) +First we need to create suitable `sf` data inputs. We can do this by +polygonizing the raster data: -features_sf <- features_sf %>% - dplyr::select(-seamounts) +``` r +#create the planning grid, which is also the same as the cost grid since we are using planning units all with cost = 1 +pu_sf <- as.polygons(pu_raster, aggregate = FALSE) |> + st_as_sf() -# Show what seamounts look like... -plot(seamounts_sf, border = F) # there are 7 seamount areas (seamounts that are touching) -``` +features_sf <- as.polygons(fish_distributions, aggregate = FALSE) |> + st_as_sf() -![](https://github.com/emlab-ucsb/patchwise/assets/40546424/892e0ee5-6791-4fb7-acac-6e3fcb8718b0) +seamounts_sf <- as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + st_as_sf() +#replace NAs with zeroes +seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 ``` -# Create seamount patches - seamount areas that touch are considered the same patch -patches_sf <- patchwise::create_patches(seamounts_sf, spatial_grid = planning_sf) -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf) +Let’s check our features and seamounts look ok: -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) +``` r +plot(cbind(features_sf, st_drop_geometry(seamounts_sf))) +``` -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf) + -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +Now we can run the same `patchwise` functions as we did with raster data +to prepare the data for prioritization -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +``` r -# Solve the prioritization -solution_sf <- solve(problem_sf) +# Create seamount patches +patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf) +#> Processing patch 1 of 3 +#> Processing patch 2 of 3 +#> Processing patch 3 of 3 -# Show the results -plot(result_sf, border = F) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf) + +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) ``` -![](https://github.com/emlab-ucsb/patchwise/assets/40546424/4fe72641-1693-43f5-a317-b86d30e7c54c) +With all the data preparation now done, we can run a prioritization +using `prioritizr`: -Areas in yellow were identified by `prioritizr` as areas worth protecting. We can see that entire seamounts were used to meet the target objective of protecting 20% of seamounts. We can compare this result to a `prioritizr` run that does not protect whole seamounts: +The solution object is a tibble. To convert this into a raster object +for plotting, we use the `convert_solution()` function from `patchwise` +``` r +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) ``` -# Grab all relevant data -features_sf_nopatch <- oceandatr::get_features(spatial_grid = planning_sf) %>% - dplyr::mutate(cost = 1) %>% # create a cost column - dplyr::relocate(cost, .before = x) # make sure cost column is before geometry column -# Run the prioritization -problem_sf_nopatch <- prioritizr::problem(x = features_sf_nopatch, features = names(features_sf_nopatch)[1:(ncol(features_sf_nopatch)-1)], cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_relative_targets(rep(0.2, ncol(features_sf_nopatch)-1)) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +For comparison, we will run a prioritization without using `patchwise` -# Solve the prioritization -solution_sf_nopatch <- solve(problem_sf_nopatch) +We can now plot the solutions, and overlaying the outlines of the +seamounts (in red), we can see that one entire seamount is included in +the solution that used `patchwise`, whereas the solution without +`patchwise` selects a few planning units in each seamount. Note that in +the solution with `patchwise`, planning units that overlap seamounts +that are not entirely selected have also been selected to meet targets +for other features (fish distributions). -# Show the results -plot(solution_sf_nopatch %>% dplyr::select(solution_1), border = F) +``` r +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(as.polygons(seamounts, aggregate = TRUE), col = "red")) ``` -![](https://github.com/emlab-ucsb/patchwise/assets/40546424/a4fee456-4206-4668-bd61-76e22b8b0222) - -Only portions of seamount units are protected here. + diff --git a/data-raw/example-data.R b/data-raw/example-data.R new file mode 100644 index 0000000..46763af --- /dev/null +++ b/data-raw/example-data.R @@ -0,0 +1,32 @@ +# Creation of a raster planning grid and some features for use in package examples + +# create a 40 x 40 planning grid, with cost values all equal to 1 + +max_x_and_y <- 40 +pu_raster <- terra::rast(nrows = max_x_and_y, ncols = max_x_and_y, xmin= 0, xmax= max_x_and_y, ymin= 0, ymax= max_x_and_y, resolution = 1 , vals = 1) |> + setNames("cost") + +#create a 4 layer raster of random, binary (0 or 1) features +# these could represent distributions of species; for this example they are named fish + +number_species <- 4 + +fish_distributions <- lapply(seq_len(number_species), function(x) { + terra::setValues(pu_raster, round(runif(terra::ncell(pu_raster)) > 0.5)) +}) |> + setNames(paste0("fish_", seq_len(number_species))) |> + terra::rast() + +# create patches; for this example, they are named seamounts +seamounts <- data.frame(x = c(5, 22, 32), + y = c(5, 18, 29)) |> + terra::vect() |> + terra::buffer(4) |> + terra::rasterize(pu_raster, field = 1) |> + setNames("Seamounts") + +terra::writeRaster(pu_raster, "inst/extdata/pu_raster.tif") + +terra::writeRaster(fish_distributions, "inst/extdata/spp_distributions.tif") + +terra::writeRaster(seamounts, "inst/extdata/seamounts.tif") diff --git a/figure/unnamed-chunk-10-1.png b/figure/unnamed-chunk-10-1.png deleted file mode 100644 index 42b235d..0000000 Binary files a/figure/unnamed-chunk-10-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-11-1.png b/figure/unnamed-chunk-11-1.png deleted file mode 100644 index d96a92e..0000000 Binary files a/figure/unnamed-chunk-11-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-12-1.png b/figure/unnamed-chunk-12-1.png deleted file mode 100644 index a4c1264..0000000 Binary files a/figure/unnamed-chunk-12-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-13-1.png b/figure/unnamed-chunk-13-1.png deleted file mode 100644 index 1b36b74..0000000 Binary files a/figure/unnamed-chunk-13-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-3-1.png b/figure/unnamed-chunk-3-1.png deleted file mode 100644 index 5b99861..0000000 Binary files a/figure/unnamed-chunk-3-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-4-1.png b/figure/unnamed-chunk-4-1.png deleted file mode 100644 index 97cc144..0000000 Binary files a/figure/unnamed-chunk-4-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-5-1.png b/figure/unnamed-chunk-5-1.png deleted file mode 100644 index d7795b1..0000000 Binary files a/figure/unnamed-chunk-5-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-6-1.png b/figure/unnamed-chunk-6-1.png deleted file mode 100644 index ef16ed0..0000000 Binary files a/figure/unnamed-chunk-6-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-7-1.png b/figure/unnamed-chunk-7-1.png deleted file mode 100644 index efa0531..0000000 Binary files a/figure/unnamed-chunk-7-1.png and /dev/null differ diff --git a/figure/unnamed-chunk-9-1.png b/figure/unnamed-chunk-9-1.png deleted file mode 100644 index a98ce2b..0000000 Binary files a/figure/unnamed-chunk-9-1.png and /dev/null differ diff --git a/inst/extdata/pu_raster.tif b/inst/extdata/pu_raster.tif new file mode 100644 index 0000000..ed7f4b3 Binary files /dev/null and b/inst/extdata/pu_raster.tif differ diff --git a/inst/extdata/seamounts.tif b/inst/extdata/seamounts.tif new file mode 100644 index 0000000..1760ba0 Binary files /dev/null and b/inst/extdata/seamounts.tif differ diff --git a/inst/extdata/spp_distributions.tif b/inst/extdata/spp_distributions.tif new file mode 100644 index 0000000..3b399fa Binary files /dev/null and b/inst/extdata/spp_distributions.tif differ diff --git a/man/constraints_targets.Rd b/man/constraints_targets.Rd index d081d82..17b0e17 100644 --- a/man/constraints_targets.Rd +++ b/man/constraints_targets.Rd @@ -2,47 +2,46 @@ % Please edit documentation in R/constraints_targets.R \name{constraints_targets} \alias{constraints_targets} -\title{Create a dataframe for manual targets for all layers, including the features and the constraints} +\title{Create a dataframe for manual targets for all layers, including the features +and the constraints} \usage{ constraints_targets(feature_targets, patch_df) } \arguments{ -\item{feature_targets}{a dataframe of feature targets generated by \code{features_targets()}} +\item{feature_targets}{a dataframe of feature targets generated by +\code{features_targets()}} -\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes constraints for each patch and grid cell combination} +\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes +constraints for each patch and grid cell combination} } \value{ A dataframe to be used to specify all manual targets for prioritizr } \description{ -This function creates a dataframe with the comibination of feature targets, generated by \code{features_targets()} and the targets for constraints to be used in prioritizr +This function creates a dataframe with the combination of +feature targets, generated by \code{features_targets()} and the targets for +constraints to be used in prioritizr } \examples{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) -# Create patch dataframe -patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, - patches = patches_raster, costs = cost_raster) -# Create boundary matrix for prioritizr -boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, - patches = patches_raster, patch_df = patches_raster_df) -# Create target features - using just 20\% for every feature -features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), - features = features_raster, pre_patches = seamounts_raster) -# Create constraint targets -constraint_targets <- constraints_targets(feature_targets = features_targets, - patch_df = patches_raster_df) +# Import some planning data + +#import planning units/ cost raster +pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) +#import species distributions - feature data for planning +species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + +# Create seamount patches as multi-layer raster +patches_raster <- create_patches(seamounts) + +patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +patches = patches_raster, costs = pu_raster) + +# Create targets for protection - use 20\% for each feature (including 20\% of whole seamounts) in this example +targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) + +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) } diff --git a/man/convert_solution.Rd b/man/convert_solution.Rd index ecd4499..2394f33 100644 --- a/man/convert_solution.Rd +++ b/man/convert_solution.Rd @@ -7,59 +7,63 @@ convert_solution(solution, patch_df, spatial_grid) } \arguments{ -\item{solution}{the solution that results from using \code{prioritizr::solve()} on \code{prioritizr::problem()}} +\item{solution}{the solution that results from using \code{prioritizr::solve()} on +\code{prioritizr::problem()}} -\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes constraints for each patch and grid cell combination} +\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes +constraints for each patch and grid cell combination} -\item{spatial_grid}{a raster or sf template with the desired resolution and coordinate reference system generated by \code{spatialgridr::get_grid()}; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object)} +\item{spatial_grid}{a raster or sf template with the desired resolution and +coordinate reference system; values in areas of interest are 1, while all +other values are NA (only required if feature is a sf object)} } \value{ -A raster or sf object that denotes areas suggested for protection (value = 1) +A raster or sf object that denotes areas suggested for protection +(value = 1) } \description{ -This function converts the prioritization solution that uses patches to a raster or sf object (depending on the format of your spatial_grid) that clearly marks areas suggested for protection (1) and areas not suggested for protection (0) +This function converts the prioritization solution that uses +patches to a raster or sf object (depending on the format of your +spatial_grid) that clearly marks areas suggested for protection (1) and +areas not suggested for protection (0) } \examples{ \dontrun{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) -# Create patch dataframe -patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, - patches = patches_raster, costs = cost_raster) -# Create boundary matrix for prioritizr -boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, patches = patches_raster, - patch_df = patches_raster_df) -# Create target features - using just 20\% for every feature -features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), - features = features_raster, pre_patches = seamounts_raster) -# Create constraint targets -constraint_targets <- constraints_targets(feature_targets = features_targets, - patch_df = patches_raster_df) -# Create prioritization problem -problem_raster <- prioritizr::problem(x = patches_raster_df, - features = constraint_targets$feature, cost_column = "cost") \%>\% - prioritizr::add_min_set_objective() \%>\% - prioritizr::add_manual_targets(constraint_targets) \%>\% - prioritizr::add_binary_decisions() \%>\% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix) \%>\% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) -# Solve problem -solution <- solve(problem_raster) -# Convert to a more digestible format -suggested_protection <- convert_solution(solution = solution, patch_df = patches_raster_df, - spatial_grid = planning_raster) +#import planning units/ cost raster +pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) +#import species distributions - feature data for planning +species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + +# Create seamount patches as multi-layer raster +patches_raster <- create_patches(seamounts) + +patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +patches = patches_raster, costs = pu_raster) + +# Create targets for protection - use 20\% for each feature (including 20\% of whole seamounts) in this example +targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) + +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) + +library(prioritizr) + +# Create the prioritization problem +problem_patches <- problem(x = patches_df, features = constraints_df$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_df) |> + add_binary_decisions() |> + add_default_solver() + +# Solve the problem +solution_patches_tbl <- solve(problem_patches) + +# Convert the solution into a raster using `convert_solution()` +sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df, spatial_grid = pu_raster) + +terra::plot(sol_rast_patches) } } diff --git a/man/create_boundary_matrix.Rd b/man/create_boundary_matrix.Rd index d63c90f..11f2080 100644 --- a/man/create_boundary_matrix.Rd +++ b/man/create_boundary_matrix.Rd @@ -7,38 +7,42 @@ create_boundary_matrix(spatial_grid, patches, patch_df) } \arguments{ -\item{spatial_grid}{a raster or sf template with the desired resolution and coordinate reference system generated by \code{spatialgridr::get_grid()}; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object)} +\item{spatial_grid}{a raster or sf template with the desired resolution and +coordinate reference system; values in areas of interest are 1, while all +other values are NA (only required if feature is a sf object)} -\item{patches}{a raster or sf object generated by \code{create_patches()}; a single feature split into spatially distinct groups; each layer of the raster or each column of the sf object identifies the location of each patch} +\item{patches}{a raster or sf object generated by \code{create_patches()}; a +single feature split into spatially distinct groups; each layer of the +raster or each column of the sf object identifies the location of each +patch} -\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes constraints for each patch and grid cell combination} +\item{patch_df}{a dataframe generated by \code{create_patch_df()} that includes +constraints for each patch and grid cell combination} } \value{ -A boundary matrix to plug into prioritizr +A sparse matrix } \description{ -This function creates a boundary matrix to plug into prioritizr +This function creates a boundary matrix which is required if you +want to use a boundary penalty in prioritizr } \examples{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) -# Create patch dataframe -patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, - patches = patches_raster, costs = cost_raster) -# Create boundary matrix for prioritizr -boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, - patches = patches_raster, patch_df = patches_raster_df) +# Import some planning data + +#import planning units/ cost raster +pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) +#import species distributions - feature data for planning +species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + +# Create seamount patches as multi-layer raster +patches_raster <- create_patches(seamounts) + +patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +patches = patches_raster, costs = pu_raster) + +#Create boundary matrix if you want to use a boundary penatly function in prioritizr +boundary_matrix <- create_boundary_matrix(spatial_grid = pu_raster, patches = patches_raster, patch_df = patches_df) } diff --git a/man/create_patch_df.Rd b/man/create_patch_df.Rd index e5768cf..e48b276 100644 --- a/man/create_patch_df.Rd +++ b/man/create_patch_df.Rd @@ -2,7 +2,8 @@ % Please edit documentation in R/create_patch_df.R \name{create_patch_df} \alias{create_patch_df} -\title{Create a dataframe that counts each patch independently for the prioritizr input} +\title{Create a dataframe that counts each patch independently for the prioritizr +input} \usage{ create_patch_df( spatial_grid, @@ -14,41 +15,54 @@ create_patch_df( ) } \arguments{ -\item{spatial_grid}{a raster or sf template with the desired resolution and coordinate reference system generated by \code{spatialgridr::get_grid()}; values in areas of interest are 1, while all other values are NA} +\item{spatial_grid}{a raster or sf template with the desired resolution and +coordinate reference system; values in areas of interest are 1, while all +other values are NA} -\item{features}{a raster or sf object that includes all relevant features to be used in the prioritization; each layer of the raster or each column of the sf object identifies the location of each feature} +\item{features}{a raster or sf object that includes all relevant features to +be used in the prioritization; each layer of the raster or each column of +the sf object identifies the location of each feature} -\item{patches}{a raster or sf object generated by \code{create_patches()}; a single feature split into spatially distinct groups; each layer of the raster or each column of the sf object identifies the location of each patch} +\item{patches}{a raster or sf object generated by \code{create_patches()}; a +single feature split into spatially distinct groups; each layer of the +raster or each column of the sf object identifies the location of each +patch} -\item{costs}{a raster or sf object of the costs for protecting each cell in the planning grid} +\item{costs}{a raster or sf object of the costs for protecting each cell in +the planning grid} -\item{locked_out}{a raster or sf object for areas to be locked out (absolutely not protected) in the prioritization} +\item{locked_out}{a raster or sf object for areas to be locked out +(absolutely not protected) in the prioritization} -\item{locked_in}{a raster or sf object for areas to be locked in (absolutely protected) in the prioritization} +\item{locked_in}{a raster or sf object for areas to be locked in (absolutely +protected) in the prioritization} } \value{ A data frame to be used as an input for \code{create_boundary_matrix()} } \description{ -This function takes all of the inputs of prioritizr and creates a dataframe precursor for prioritizr. This specific function includes the use of "patches" which have been already generated by \code{create_patches()}. "Patches" are spatial groups of a particular feature, such as seamounts. Prioritizr will protect entire patches to meet the protection goal for that feature. The output of this function should be passed to \code{create_boundary_matrix()} before being passed to prioritizr. +This function takes all of the inputs of prioritizr and creates +a dataframe precursor for prioritizr. This specific function includes the +use of "patches" which have been already generated by \code{create_patches()}. +"Patches" are spatial groups of a particular feature, such as seamounts. +Prioritizr will protect entire patches to meet the protection goal for that +feature. The output of this function should be passed to +\code{create_boundary_matrix()} before being passed to prioritizr. } \examples{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) -# Create patch dataframe -patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, - patches = patches_raster, costs = cost_raster) +# Import some planning data + +#import planning units/ cost raster +pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) +#import species distributions - feature data for planning +species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + +# Create seamount patches as multi-layer raster +patches_raster <- create_patches(seamounts) + +patches_df <- create_patch_df(spatial_grid = pu_raster, features = species_distributions, +patches = patches_raster, costs = pu_raster) } diff --git a/man/create_patches.Rd b/man/create_patches.Rd index 508a35a..791ac89 100644 --- a/man/create_patches.Rd +++ b/man/create_patches.Rd @@ -7,30 +7,31 @@ create_patches(feature, spatial_grid = NULL) } \arguments{ -\item{feature}{a raster or sf object with the feature of interest present (value = 1) or absent (value = NA)} +\item{feature}{a raster or sf object with the feature of interest present +(value = 1) or absent (value = NA)} -\item{spatial_grid}{a raster or sf template with the desired resolution and coordinate reference system generated by \code{spatialgridr::get_grid()}; values in areas of interest are 1, while all other values are NA (only required if feature is a sf object)} +\item{spatial_grid}{a raster or sf template with the desired resolution and +coordinate reference system generated; values in areas of interest are 1, +while all other values are NA (only required if feature is a sf object)} } \value{ -A raster or sf object with independent layers (raster)/columns (sf) designating the location of each patch +A raster or sf object with independent layers (raster)/columns (sf) +designating the location of each patch } \description{ -This function takes a feature and splits it into multiple "patches" of the feature. This is ideal for features that you want to protect x\% of, but want to make sure that entire patches (as opposed to pieces of patches) are protected to meet that target. +This function takes a feature and splits it into multiple +"patches" of the feature. This is ideal for features that you want to +protect x\% of, but want to make sure that entire patches (as opposed to +pieces of patches) are protected to meet that target. } \examples{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) +# Import some planning data +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + +terra::plot(seamounts) + +# Create seamount patches as multi-layer raster +patches_raster <- create_patches(seamounts) +terra::plot(patches_raster) } diff --git a/man/features_targets.Rd b/man/features_targets.Rd index 853e0df..ce5d8aa 100644 --- a/man/features_targets.Rd +++ b/man/features_targets.Rd @@ -13,45 +13,43 @@ features_targets( ) } \arguments{ -\item{targets}{a vector of targets for protection (range between 0 and 1); must be the same length as the number of features + the pre-patches variable} +\item{targets}{a vector of targets for protection (range between 0 and 1); +must be the same length as the number of features + the pre-patches +variable} -\item{features}{a raster or sf object that includes all relevant features to be used in the prioritization; each layer of the raster or each column of the sf object identifies the location of each feature} +\item{features}{a raster or sf object that includes all relevant features to +be used in the prioritization; each layer of the raster or each column of +the sf object identifies the location of each feature} -\item{pre_patches}{a raster or sf object that includes the feature that is to be split into patches (not the layer that is already split into patches)} +\item{pre_patches}{a raster or sf object that includes the feature that is to +be split into patches (not the layer that is already split into patches)} -\item{locked_out}{a raster or sf object for areas to be locked out (absolutely not protected) in the prioritization} +\item{locked_out}{a raster or sf object for areas to be locked out +(absolutely not protected) in the prioritization} -\item{locked_in}{a raster or sf object for areas to be locked in (absolutely protected) in the prioritization} +\item{locked_in}{a raster or sf object for areas to be locked in (absolutely +protected) in the prioritization} } \value{ -A data frame to be used to specify targets for features. Will need to be plugged into \code{constraint_targets()} before being implemented into prioritizr +A data frame to be used to specify targets for features. Will need to +be plugged into \code{constraint_targets()} before being implemented into +prioritizr } \description{ -This function creates a dataframe of targets for each of the features to be used in the prioritization +This function creates a dataframe of targets for each of the +features to be used in the prioritization } \examples{ -# Start with a little housekeeping to get the data from oceandatr -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <-'+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' -# Create a planning grid -planning_raster <- spatialgridr::get_grid(area, projection = projection) -# Grab all relevant data -features_raster <- oceandatr::get_features(spatial_grid = planning_raster) -# Separate seamount data - we want to protect entire patches -seamounts_raster <- features_raster[["seamounts"]] -features_raster <- features_raster[[names(features_raster)[names(features_raster) != "seamounts"]]] -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_raster <- stats::setNames(planning_raster, "cost") -# Create patches from layer -patches_raster <- create_patches(seamounts_raster) -# Create patch dataframe -patches_raster_df <- create_patch_df(spatial_grid = planning_raster, features = features_raster, - patches = patches_raster, costs = cost_raster) -# Create boundary matrix for prioritizr -boundary_matrix <- create_boundary_matrix(spatial_grid = planning_raster, patches = patches_raster, - patch_df = patches_raster_df) -# Create target features - using just 20\% for every feature -features_targets <- features_targets(targets = rep(0.2, (terra::nlyr(features_raster)) + 1), - features = features_raster, pre_patches = seamounts_raster) +# Import some planning data + +#import species distributions - feature data for planning +species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + +#import seamounts +seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + + +# Create targets for protection - use 20\% for each feature (including 20\% of whole seamounts) in this example +targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), features = species_distributions, pre_patches = seamounts) +head(targets_df) } diff --git a/man/figures/README-prioritization-boundary-plot-1.png b/man/figures/README-prioritization-boundary-plot-1.png new file mode 100644 index 0000000..3928e9e Binary files /dev/null and b/man/figures/README-prioritization-boundary-plot-1.png differ diff --git a/man/figures/README-prioritization-plot-1.png b/man/figures/README-prioritization-plot-1.png new file mode 100644 index 0000000..8df690a Binary files /dev/null and b/man/figures/README-prioritization-plot-1.png differ diff --git a/man/figures/README-prioritization-plot-sf-1.png b/man/figures/README-prioritization-plot-sf-1.png new file mode 100644 index 0000000..2cc13e5 Binary files /dev/null and b/man/figures/README-prioritization-plot-sf-1.png differ diff --git a/man/figures/README-unnamed-chunk-5-1.png b/man/figures/README-unnamed-chunk-5-1.png new file mode 100644 index 0000000..be16b98 Binary files /dev/null and b/man/figures/README-unnamed-chunk-5-1.png differ diff --git a/man/figures/README-unnamed-chunk-6-1.png b/man/figures/README-unnamed-chunk-6-1.png new file mode 100644 index 0000000..be16b98 Binary files /dev/null and b/man/figures/README-unnamed-chunk-6-1.png differ diff --git a/man/figures/README-visualize-data-1.png b/man/figures/README-visualize-data-1.png new file mode 100644 index 0000000..f8e5983 Binary files /dev/null and b/man/figures/README-visualize-data-1.png differ diff --git a/tests/testthat/test-constraints_targets.R b/tests/testthat/test-constraints_targets.R index 2a63f63..3091080 100644 --- a/tests/testthat/test-constraints_targets.R +++ b/tests/testthat/test-constraints_targets.R @@ -1,83 +1,49 @@ test_that("create constraints - raster", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] + patches_raster <- create_patches(seamounts) - # Create seamount patches - seamount areas that touch are considered the same patch - patches_rast <- create_patches(seamounts_rast) + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_raster, features = species_distributions, + patches = patches_raster, costs = pu_raster)) - # Create a "cost" to protecting a cell - just a uniform cost for this example - cost_rast <- stats::setNames(planning_rast, "cost") + targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), + features = species_distributions, pre_patches = seamounts) - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_rast <- create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast) - - # Create boundary matrix - bnd_mat_rast <- create_boundary_matrix(spatial_grid = planning_rast, - patches = patches_rast, - patch_df = patches_df_rast) - - # Create targets - targets_rast <- features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), - features = features_rast, - pre_patches = seamounts_rast) - - # Create constraints - constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) - }, + constraints_targets(feature_targets = targets_df, patch_df = patches_df)}, class = "tbl_df") }) test_that("create constraints - sf", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) + pu_sf <- terra::as.polygons(pu_raster, aggregate = FALSE) |> + sf::st_as_sf() - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) - features_sf <- features_sf %>% dplyr::select(-seamounts) + features_sf <- terra::as.polygons(species_distributions, aggregate = FALSE) |> + sf::st_as_sf() - # Create seamount patches - seamount areas that touch are considered the same patch - patches_sf <- create_patches(seamounts_sf, spatial_grid = planning_sf) + seamounts_sf <- terra::as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() - # Add cost layer - cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_square <- create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf) + patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) - # Create boundary matrix - bnd_mat_square <- create_boundary_matrix(spatial_grid = planning_sf, - patches = patches_sf, - patch_df = patches_df_square) + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_sf, features = features_sf, + patches = patches_sf, costs = pu_sf)) - # Create features targets - targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), - features = features_sf, - pre_patches = seamounts_sf) + targets_df <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf) - # Create constraints - constraints_targets(feature_targets = targets_sf, patch_df = patches_df_square) - }, + constraints_targets(feature_targets = targets_df, patch_df = patches_df)}, class = "tbl_df") }) diff --git a/tests/testthat/test-convert_solution.R b/tests/testthat/test-convert_solution.R deleted file mode 100644 index 30bf524..0000000 --- a/tests/testthat/test-convert_solution.R +++ /dev/null @@ -1,111 +0,0 @@ -test_that("convert solution - raster", { - expect_s4_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' - - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) - - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) - - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] - - # Create seamount patches - seamount areas that touch are considered the same patch - patches_rast <- create_patches(seamounts_rast) - - # Create a "cost" to protecting a cell - just a uniform cost for this example - cost_rast <- stats::setNames(planning_rast, "cost") - - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_rast <- create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast) - - # Create boundary matrix - bnd_mat_rast <- create_boundary_matrix(spatial_grid = planning_rast, - patches = patches_rast, - patch_df = patches_df_rast) - - # Create targets - targets_rast <- features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), - features = features_rast, - pre_patches = seamounts_rast) - - # Create constraints - constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) - - # Run prioritization - p_rast <- suppressWarnings(prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = bnd_mat_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1)) - - # Solve - raw_solution <- solve(p_rast) - - # Convert solution into more digestible format - convert_solution(raw_solution, patches_df_rast, planning_rast) - }, - class = "SpatRaster") -}) - -test_that("create constraints - sf", { - expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' - - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") - - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) - - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) - features_sf <- features_sf %>% dplyr::select(-seamounts) - - # Create seamount patches - seamount areas that touch are considered the same patch - patches_sf <- create_patches(seamounts_sf, spatial_grid = planning_sf) - - # Add cost layer - cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) - - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_square <- create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf) - - # Create boundary matrix - bnd_mat_square <- create_boundary_matrix(spatial_grid = planning_sf, - patches = patches_sf, - patch_df = patches_df_square) - - # Create features targets - targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), - features = features_sf, - pre_patches = seamounts_sf) - - # Create constraints - constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_square) - - # Run prioritization - p_square <- suppressWarnings(prioritizr::problem(x = patches_df_square, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = bnd_mat_square) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1)) - - # Solve - raw_solution <- solve(p_square) - - # Convert solution - convert_solution(raw_solution, patches_df_square, planning_sf) - }, - class = "sf") -}) diff --git a/tests/testthat/test-convert_solution.R.notworking b/tests/testthat/test-convert_solution.R.notworking new file mode 100644 index 0000000..6f69003 --- /dev/null +++ b/tests/testthat/test-convert_solution.R.notworking @@ -0,0 +1,71 @@ +test_that("convert solution - raster", { + expect_s4_class( + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + + patches_raster <- create_patches(seamounts) + + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_raster, features = species_distributions, + patches = patches_raster, costs = pu_raster)) + + targets_df <- features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), + features = species_distributions, pre_patches = seamounts) + + constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) + + problem_patches <- prioritizr::problem(x = patches_df, features = constraints_df$feature, cost_column = "cost") |> + prioritizr::add_min_set_objective() |> + prioritizr::add_manual_targets(constraints_df) |> + prioritizr::add_binary_decisions() |> + prioritizr::add_gurobi_solver() + + solution_patches_tbl <- suppressMessages(solve(problem_patches)) + + convert_solution(solution = solution_patches_tbl, patch_df = patches_df, spatial_grid = pu_raster)}, + class = "SpatRaster") +}) + +test_that("create constraints - sf", { + expect_s3_class( + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) + + pu_sf <- terra::as.polygons(pu_raster, aggregate = FALSE) |> + sf::st_as_sf() + + features_sf <- terra::as.polygons(species_distributions, aggregate = FALSE) |> + sf::st_as_sf() + + seamounts_sf <- terra::as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() + + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 + + patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) + + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_sf, features = features_sf, + patches = patches_sf, costs = pu_sf)) + + targets_df <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf) + + constraints_df <- constraints_targets(feature_targets = targets_df, patch_df = patches_df) + + problem_patches <- prioritizr::problem(x = patches_df, features = constraints_df$feature, cost_column = "cost") |> + prioritizr::add_min_set_objective() |> + prioritizr::add_manual_targets(constraints_df) |> + prioritizr::add_binary_decisions() |> + prioritizr::add_gurobi_solver() + + solution_patches_tbl <- suppressMessages(solve(problem_patches)) + + convert_solution(solution = solution_patches_tbl, patch_df = patches_df, spatial_grid = pu_sf) + + }, + class = "sf") +}) diff --git a/tests/testthat/test-create_boundary_matrix.R b/tests/testthat/test-create_boundary_matrix.R index 3b2c429..2dcb7fc 100644 --- a/tests/testthat/test-create_boundary_matrix.R +++ b/tests/testthat/test-create_boundary_matrix.R @@ -1,68 +1,44 @@ test_that("create boundary matrix - raster", { expect_s4_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] + patches_raster <- create_patches(seamounts) - # Create seamount patches - seamount areas that touch are considered the same patch - patches_rast <- create_patches(seamounts_rast) + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_raster, features = species_distributions, + patches = patches_raster, costs = pu_raster)) - # Create a "cost" to protecting a cell - just a uniform cost for this example - cost_rast <- stats::setNames(planning_rast, "cost") - - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_rast <- create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast) - - # Create boundary matrix - create_boundary_matrix(spatial_grid = planning_rast, - patches = patches_rast, - patch_df = patches_df_rast) - - }, + create_boundary_matrix(spatial_grid = pu_raster, patches = patches_raster, patch_df = patches_df)}, class = "Matrix") }) test_that("create boundary matrix - sf", { expect_s4_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") + pu_sf <- terra::as.polygons(pu_raster, aggregate = FALSE) |> + sf::st_as_sf() - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) + features_sf <- terra::as.polygons(species_distributions, aggregate = FALSE) |> + sf::st_as_sf() - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) - features_sf <- features_sf %>% dplyr::select(-seamounts) + seamounts_sf <- terra::as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() - # Create seamount patches - seamount areas that touch are considered the same patch - patches_sf <- create_patches(seamounts_sf, spatial_grid = planning_sf) + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 - # Add cost layer - cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) + patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_square <- create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf) + patches_df <- suppressMessages(create_patch_df(spatial_grid = pu_sf, features = features_sf, + patches = patches_sf, costs = pu_sf)) - # Create boundary matrix - create_boundary_matrix(spatial_grid = planning_sf, - patches = patches_sf, - patch_df = patches_df_square) - }, + create_boundary_matrix(spatial_grid = pu_sf, patches = patches_sf, patch_df = patches_df)}, class = "Matrix") }) diff --git a/tests/testthat/test-create_patch_df.R b/tests/testthat/test-create_patch_df.R index 6b41240..1378b93 100644 --- a/tests/testthat/test-create_patch_df.R +++ b/tests/testthat/test-create_patch_df.R @@ -1,57 +1,41 @@ test_that("convert seamount data to df - raster", { expect_s3_class( # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] + patches_raster <- create_patches(seamounts) - # Create seamount patches - seamount areas that touch are considered the same patch - patches_rast <- create_patches(seamounts_rast) - - # Create a "cost" to protecting a cell - just a uniform cost for this example - cost_rast <- stats::setNames(planning_rast, "cost") - - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - suppressMessages(create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast)) - }, + suppressMessages(create_patch_df(spatial_grid = pu_raster, features = species_distributions, + patches = patches_raster, costs = pu_raster))}, class = "tbl_df") }) test_that("convert seamount data to df - sf", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {pu_raster <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise")) + + species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) + + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") + pu_sf <- terra::as.polygons(pu_raster, aggregate = FALSE) |> + sf::st_as_sf() - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) + features_sf <- terra::as.polygons(species_distributions, aggregate = FALSE) |> + sf::st_as_sf() - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) - features_sf <- features_sf %>% dplyr::select(-seamounts) + seamounts_sf <- terra::as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() - # Create seamount patches - seamount areas that touch are considered the same patch - patches_sf <- create_patches(seamounts_sf, spatial_grid = planning_sf) + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 - # Add cost layer - cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) + patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - suppressMessages(create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf)) - }, + suppressMessages(create_patch_df(spatial_grid = pu_sf, features = features_sf, + patches = patches_sf, costs = pu_sf))}, class = "tbl_df") }) diff --git a/tests/testthat/test-create_patches.R b/tests/testthat/test-create_patches.R index bd2319d..dd3665a 100644 --- a/tests/testthat/test-create_patches.R +++ b/tests/testthat/test-create_patches.R @@ -1,39 +1,23 @@ test_that("create patches from seamount data - raster", { expect_s4_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' - - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) - - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) - - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - - suppressMessages(create_patches(seamounts_rast)) - }, + {terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) |> + create_patches(seamounts)}, class = "SpatRaster") }) test_that("create patches from seamount data - sf", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' - - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") + {pu_sf <- terra::rast(system.file("extdata/pu_raster.tif", package = "patchwise"))|> + terra::as.polygons(aggregate = FALSE) |> + sf::st_as_sf() - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) + seamounts_sf <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) |> + terra::as.polygons(aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) + #replace NAs with zeroes + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 - suppressMessages(create_patches(seamounts_sf, spatial_grid = planning_sf)) - }, + create_patches(seamounts_sf, spatial_grid = pu_sf)}, class = "sf") }) diff --git a/tests/testthat/test-features_targets.R b/tests/testthat/test-features_targets.R index 950ed01..959d6da 100644 --- a/tests/testthat/test-features_targets.R +++ b/tests/testthat/test-features_targets.R @@ -1,79 +1,29 @@ test_that("create features targets - raster", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' + {species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Create a planning grid - planning_rast <- spatialgridr::get_grid(area, projection = projection) + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Grab all relevant data - features_rast <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_rast))) - - # Separate seamount data - we want to protect entire patches - seamounts_rast <- features_rast[["seamounts"]] - features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] - - # Create seamount patches - seamount areas that touch are considered the same patch - patches_rast <- create_patches(seamounts_rast) - - # Create a "cost" to protecting a cell - just a uniform cost for this example - cost_rast <- stats::setNames(planning_rast, "cost") - - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_rast <- create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast) - - # Create boundary matrix - bnd_mat_rast <- create_boundary_matrix(spatial_grid = planning_rast, - patches = patches_rast, - patch_df = patches_df_rast) - - # Create targets - features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), - features = features_rast, - pre_patches = seamounts_rast) - - }, + features_targets(targets = rep(0.2, (terra::nlyr(species_distributions) + 1)), + features = species_distributions, pre_patches = seamounts)}, class = "tbl_df") }) test_that("create features targets - sf", { expect_s3_class( - # Choose area of interest (Bermuda EEZ) - {area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") - projection <- 'PROJCS["ProjWiz_Custom_Lambert_Azimuthal", GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0], UNIT["Degree",0.0174532925199433]], PROJECTION["Lambert_Azimuthal_Equal_Area"], PARAMETER["False_Easting",0.0], PARAMETER["False_Northing",0.0], PARAMETER["Central_Meridian",-64.5], PARAMETER["Latitude_Of_Origin",32], UNIT["Meter",1.0]]' - - # Create a planning grid - planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") - - # Grab all relevant data - features_sf <- suppressMessages(suppressWarnings(oceandatr::get_features(spatial_grid = planning_sf))) - - # Separate seamount data - we want to protect entire patches - seamounts_sf <- features_sf %>% dplyr::select(seamounts) - features_sf <- features_sf %>% dplyr::select(-seamounts) + {species_distributions <- terra::rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) - # Create seamount patches - seamount areas that touch are considered the same patch - patches_sf <- create_patches(seamounts_sf, spatial_grid = planning_sf) + seamounts <- terra::rast(system.file("extdata/seamounts.tif", package = "patchwise")) - # Add cost layer - cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) + features_sf <- terra::as.polygons(species_distributions, aggregate = FALSE) |> + sf::st_as_sf() - # Create patches dataframe - this creates several constraints so that entire seamount units are protected together - patches_df_square <- create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf) + seamounts_sf <- terra::as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + sf::st_as_sf() - # Create boundary matrix - bnd_mat_square <- create_boundary_matrix(spatial_grid = planning_sf, - patches = patches_sf, - patch_df = patches_df_square) + seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 - # Create features targets - features_targets(targets = rep(0.2, ncol(features_sf)), - features = features_sf, - pre_patches = seamounts_sf) - }, + features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf)}, class = "tbl_df") }) diff --git a/vignettes/figure/prioritization-plot-1.png b/vignettes/figure/prioritization-plot-1.png new file mode 100644 index 0000000..eb435b7 Binary files /dev/null and b/vignettes/figure/prioritization-plot-1.png differ diff --git a/vignettes/figure/prioritization-plot-sf-1.png b/vignettes/figure/prioritization-plot-sf-1.png new file mode 100644 index 0000000..6b8cc13 Binary files /dev/null and b/vignettes/figure/prioritization-plot-sf-1.png differ diff --git a/vignettes/figure/prioritization-plot-sf2-1.png b/vignettes/figure/prioritization-plot-sf2-1.png new file mode 100644 index 0000000..25e2acf Binary files /dev/null and b/vignettes/figure/prioritization-plot-sf2-1.png differ diff --git a/vignettes/figure/prioritization-plot2-1.png b/vignettes/figure/prioritization-plot2-1.png new file mode 100644 index 0000000..9107484 Binary files /dev/null and b/vignettes/figure/prioritization-plot2-1.png differ diff --git a/vignettes/figure/unnamed-chunk-6-1.png b/vignettes/figure/unnamed-chunk-6-1.png index ef16ed0..8712454 100644 Binary files a/vignettes/figure/unnamed-chunk-6-1.png and b/vignettes/figure/unnamed-chunk-6-1.png differ diff --git a/vignettes/figure/unnamed-chunk-8-1.png b/vignettes/figure/unnamed-chunk-8-1.png new file mode 100644 index 0000000..d6bca76 Binary files /dev/null and b/vignettes/figure/unnamed-chunk-8-1.png differ diff --git a/vignettes/figure/visualize-data-1.png b/vignettes/figure/visualize-data-1.png new file mode 100644 index 0000000..2b8e69b Binary files /dev/null and b/vignettes/figure/visualize-data-1.png differ diff --git a/vignettes/locked_in_out.Rmd b/vignettes/locked_in_out.Rmd index 76b0208..585015e 100644 --- a/vignettes/locked_in_out.Rmd +++ b/vignettes/locked_in_out.Rmd @@ -27,775 +27,354 @@ A quick note about `prioritizr` - sometimes there are several cells that are equ To generate our data for this example, we will use data from Bermuda generated via `oceandatr`. -```r +``` r # Load libraries library(patchwise) +library(terra) +library(sf) +library(prioritizr) +``` -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <- '+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' +We will import some basic planning data, same as in the `patchwise` readme. For this example we will use a a 40 x 40 raster planning grid which has a cost value of 1 for each cell, and the following features that will be targeted in the prioritization: -# Create a planning grid -planning_rast <- spatialgridr::get_grid(area, projection = projection) +- 4 random, binary (0 or 1) data layers, representing species distributions; for this example we are calling them fish +- 3 "patches" of contiguous features; for this example we are calling them seamounts -# Grab all relevant data -features_rast <- oceandatr::get_features(spatial_grid = planning_rast) -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_rast <- stats::setNames(planning_rast, "cost") +``` r +#import planning units/ cost raster +pu_raster <- rast(system.file("extdata/pu_raster.tif", package = "patchwise")) -# Separate seamount data - we want to protect entire patches -seamounts_rast <- features_rast[["seamounts"]] -features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] -``` +#import fish distributions +fish_distributions <- rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) -We can generate a locked-in layer that would represent an existing MPA. - - -```r -# Create a function to produce a square area -create_area <- function(spatial_grid, mpa_size, edge_offset) { - if(class(spatial_grid)[1] %in% c("RasterLayer", "SpatRaster")) { - terra::rast(spatial_grid, - vals = c(rep(NA, (terra::ncol(spatial_grid)*edge_offset)), - rep(c(rep(NA, edge_offset), - rep(1, mpa_size), - rep(NA, (terra::ncol(spatial_grid)-(edge_offset+mpa_size)))), - mpa_size), - rep(NA, (terra::ncol(spatial_grid)*(terra::nrow(spatial_grid)-edge_offset-mpa_size))))) - } else if(class(spatial_grid)[1] == "sf") { - sf_to_grid <- sf::st_make_grid(spatial_grid, cellsize = sqrt(sf::st_area(spatial_grid[1,]))) - sf_centroid <- sf::st_centroid(sf_to_grid) %>% - sf::st_coordinates(.) %>% - as.data.frame() %>% - dplyr::mutate(X = round(X, digits = 4), - Y = round(Y, digits = 4)) - - # sf::st_make_grid changes the orientation of how cells are read - change them back - sf_to_grid <- sf_to_grid %>% - as.data.frame() %>% - dplyr::bind_cols(sf_centroid) %>% - dplyr::arrange(dplyr::desc(Y), X) %>% - dplyr::select(-X, -Y) - - cols_sf_to_grid <- length(unique(sf_centroid$X)) - rows_sf_to_grid <- length(unique(sf_centroid$Y)) - - sf_to_grid <- sf_to_grid %>% - sf::st_sf() %>% - dplyr::mutate(layer = c(rep(NA, (cols_sf_to_grid*edge_offset)), - rep(c(rep(NA, edge_offset), - rep(1, mpa_size), - rep(NA, (cols_sf_to_grid-(edge_offset+mpa_size)))), - mpa_size), - rep(NA, (cols_sf_to_grid*(rows_sf_to_grid-edge_offset-mpa_size))))) - - # the planning sf only includes polygons for the planning area (we need to remove them) so the number of cells matches the planning grid - sf_to_grid <- sf_to_grid[sf_to_grid$geometry %in% spatial_grid[[1]], ] - rownames(sf_to_grid) <- NULL - sf_to_grid %>% - dplyr::relocate(layer, 1) - } -} +#import seamounts +seamounts <- rast(system.file("extdata/seamounts.tif", package = "patchwise")) ``` +We will generate a square within our planning area that will be used as our locked in area -```r -# Create a fake MPA -mpa_location <- create_area(planning_rast, 20, 50) -# Plot fake MPA -terra::plot(mpa_location) +``` r +locked_area <- as.polygons(ext(5, 15, 25, 35)) |> + rasterize(pu_raster, background = 0) |> + setNames("Locked-in") ``` -![Map of potential locked-in areas.](figure/unnamed-chunk-3-1.png) +Let's look at our fish distributions, seamounts and locked-in area: -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). +``` r +plot(c(fish_distributions, seamounts, locked_area)) +``` -```r -# Create seamount patches - seamount areas that touch are considered the same patch -patches_rast <- patchwise::create_patches(seamounts_rast) +![plot of chunk visualize-data](figure/visualize-data-1.png) -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_in = mpa_location) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts) - this is the same process as shown in the main `patchwise` readme, but with the locked-in area. -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_in = mpa_location) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) +``` r +# Create seamount patches +patches_rast <- create_patches(seamounts) -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_rast <- create_patch_df(spatial_grid = pu_raster, features = fish_distributions, patches = patches_rast, costs = pu_raster, locked_in = locked_area) -# Solve the prioritization -solution_rast <- solve(problem_rast) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76296 rows, 56073 columns and 214386 nonzeros -#> Model fingerprint: 0xc25c3da2 -#> Variable types: 0 continuous, 56073 integer (56073 binary) -#> Coefficient statistics: -#> Matrix range [1e+00, 5e+02] -#> Objective range [2e-02, 5e+02] -#> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3740.9750000 -#> Presolve removed 2915 rows and 1615 columns -#> Presolve time: 0.64s -#> Presolved: 73381 rows, 54458 columns, 207121 nonzeros -#> Variable types: 0 continuous, 54458 integer (54458 binary) -#> Root relaxation presolved: 73381 rows, 54458 columns, 207121 nonzeros -#> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... -#> -#> Warning: Markowitz tolerance tightened to 0.25 -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 75523 3.7292125e+03 0.000000e+00 4.676925e+03 5s -#> 88193 3.7247051e+03 0.000000e+00 1.152036e+03 10s -#> 92763 3.7233916e+03 0.000000e+00 6.880988e+03 15s -#> 97523 3.7230699e+03 0.000000e+00 1.977410e+03 20s -#> 103193 3.7230205e+03 0.000000e+00 4.106123e+04 25s -#> 109387 3.7228462e+03 0.000000e+00 3.799363e+04 30s -#> 113716 3.7228405e+03 0.000000e+00 5.063170e+02 35s -#> 117511 3.7228325e+03 0.000000e+00 5.532806e+03 40s -#> 119951 3.7228309e+03 0.000000e+00 1.711754e+03 45s -#> 122461 3.7228286e+03 0.000000e+00 1.202599e+03 50s -#> Concurrent spin time: 0.04s -#> -#> Solved with dual simplex -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 69647 3.7229317e+03 0.000000e+00 0.000000e+00 51s -#> -#> Root relaxation: objective 3.722932e+03, 69647 iterations, 50.41 seconds (151.93 work units) -#> -#> Nodes | Current Node | Objective Bounds | Work -#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time -#> -#> 0 0 3722.93174 0 50601 3740.97500 3722.93174 0.48% - 51s -#> -#> Explored 1 nodes (69735 simplex iterations) in 51.50 seconds (153.73 work units) -#> Thread count was 9 (of 10 available processors) -#> -#> Solution count 1: 3740.98 -#> -#> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.740975000000e+03, best bound 3.722931735863e+03, gap 0.4823% - -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_rast <- features_targets(targets = rep(0.2, (nlyr(fish_distributions) + 1)), features = fish_distributions, pre_patches = seamounts, locked_in = locked_area) -# Show the results -terra::plot(result_rast) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) ``` -![Map of proposed solution with locked-in areas.](figure/unnamed-chunk-4-1.png) +With all the data preparation now done, we can run a prioritization using `prioritizr`: -The areas in green are those selected by the prioritization for protection. We can clearly see the locked-in area (the green square) was included! -# Implement locked-out areas using a `raster` input -What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude this area from the prioritization results. +For comparison, we will run a prioritization without using `patchwise` -```r -# Create a fake locked-out area -no_protection <- create_area(planning_rast, 20, 95) -# Plot fake locked-out area -terra::plot(no_protection) -``` +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The square, locked-in area is included in both scenarios. -![Map of potential locked-out areas.](figure/unnamed-chunk-5-1.png) - -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). +``` r +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(as.polygons(seamounts), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) +``` -```r -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_out = no_protection) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" +![plot of chunk prioritization-plot](figure/prioritization-plot-1.png) -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +# Implement locked-out areas using a `raster` input -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_out = no_protection) +What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude the "locked_area" from the prioritization results. -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) -# Solve the prioritization -solution_rast <- solve(problem_rast) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76296 rows, 56073 columns and 214385 nonzeros -#> Model fingerprint: 0xf7b92fc5 -#> Variable types: 0 continuous, 56073 integer (56073 binary) -#> Coefficient statistics: -#> Matrix range [1e+00, 5e+02] -#> Objective range [2e-02, 5e+02] -#> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3733.8850000 -#> Presolve removed 2633 rows and 1557 columns -#> Presolve time: 0.57s -#> Presolved: 73663 rows, 54516 columns, 207723 nonzeros -#> Variable types: 0 continuous, 54516 integer (54516 binary) -#> Root relaxation presolved: 73663 rows, 54516 columns, 207723 nonzeros -#> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... -#> -#> Warning: Markowitz tolerance tightened to 0.5 -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 80813 3.7266671e+03 0.000000e+00 1.275446e+03 5s -#> 89493 3.7240545e+03 0.000000e+00 1.496065e+03 10s -#> 94678 3.7234499e+03 0.000000e+00 1.518754e+03 15s -#> 98708 3.7226503e+03 0.000000e+00 1.557963e+03 20s -#> 103033 3.7218083e+03 0.000000e+00 2.200248e+03 25s -#> 108113 3.7216433e+03 0.000000e+00 7.935639e+03 30s -#> Concurrent spin time: 0.00s -#> -#> Solved with dual simplex -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 55179 3.7213346e+03 0.000000e+00 0.000000e+00 35s -#> -#> Root relaxation: objective 3.721335e+03, 55179 iterations, 34.54 seconds (101.60 work units) -#> -#> Nodes | Current Node | Objective Bounds | Work -#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time -#> -#> 0 0 3721.33461 0 49528 3733.88500 3721.33461 0.34% - 35s -#> -#> Explored 1 nodes (55285 simplex iterations) in 35.60 seconds (103.27 work units) -#> Thread count was 9 (of 10 available processors) -#> -#> Solution count 1: 3733.89 -#> -#> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.733885000000e+03, best bound 3.721335000000e+03, gap 0.3361% +We can now plot the solutions, and overlaying the outlines of the seamounts and locked-out area (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The square, locked-out area (red outline) is excluded in both scenarios. -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) -# Show the results -terra::plot(result_rast) +``` r +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(rbind(as.polygons(seamounts), as.polygons(locked_area)), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) ``` -![Map of proposed solution with locked-out areas.](figure/unnamed-chunk-6-1.png) - -Great! The area that we don't want to include in the solution is removed! - -# Implement locked-in and locked-out areas using a `raster` input +![plot of chunk prioritization-plot2](figure/prioritization-plot2-1.png) -We can also run a scenario in which we have both locked-in and locked-out areas. - - -```r -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_in = mpa_location, locked_out = no_protection) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" +# Implement locked-in areas using an `sf` input -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +We can repeat the examples using `sf` input data rather than raster. First we polygonize the raster inputs -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_in = mpa_location, locked_out = no_protection) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) +``` r +#create the planning grid, which is also the same as the cost grid since we are using planning units all with cost = 1 +pu_sf <- as.polygons(pu_raster, aggregate = FALSE) |> + st_as_sf() -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +features_sf <- as.polygons(fish_distributions, aggregate = FALSE) |> + st_as_sf() -# Solve the prioritization -solution_rast <- solve(problem_rast) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76297 rows, 56073 columns and 214786 nonzeros -#> Model fingerprint: 0x8acc3e12 -#> Variable types: 0 continuous, 56073 integer (56073 binary) -#> Coefficient statistics: -#> Matrix range [1e+00, 5e+02] -#> Objective range [2e-02, 5e+02] -#> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3741.0600000 -#> Presolve removed 4596 rows and 2855 columns -#> Presolve time: 0.63s -#> Presolved: 71701 rows, 53218 columns, 202309 nonzeros -#> Variable types: 0 continuous, 53218 integer (53218 binary) -#> Root relaxation presolved: 71701 rows, 53218 columns, 202309 nonzeros -#> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... -#> -#> Warning: Markowitz tolerance tightened to 0.5 -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 78109 3.7276375e+03 0.000000e+00 1.808055e+03 5s -#> 86897 3.7269661e+03 0.000000e+00 3.529364e+02 10s -#> 97327 3.7256267e+03 0.000000e+00 1.437909e+03 15s -#> 107276 3.7231045e+03 0.000000e+00 3.567818e+03 20s -#> 111807 3.7230532e+03 0.000000e+00 3.981599e+03 25s -#> 116500 3.7230355e+03 0.000000e+00 4.431877e+03 30s -#> 119866 3.7230222e+03 0.000000e+00 2.720234e+03 35s -#> 123492 3.7230029e+03 0.000000e+00 8.703959e+03 40s -#> 126136 3.7229923e+03 0.000000e+00 3.281546e+03 45s -#> Concurrent spin time: 0.00s -#> -#> Solved with dual simplex -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 64341 3.7230710e+03 0.000000e+00 0.000000e+00 49s -#> -#> Root relaxation: objective 3.723071e+03, 64341 iterations, 48.49 seconds (144.87 work units) -#> -#> Nodes | Current Node | Objective Bounds | Work -#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time -#> -#> 0 0 3723.07098 0 49364 3741.06000 3723.07098 0.48% - 49s -#> -#> Explored 1 nodes (64433 simplex iterations) in 49.59 seconds (146.66 work units) -#> Thread count was 9 (of 10 available processors) -#> -#> Solution count 1: 3741.06 -#> -#> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.741060000000e+03, best bound 3.723070980814e+03, gap 0.4809% +seamounts_sf <- as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + st_as_sf() -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +#replace NAs with zeroes +seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 -# Show the results -terra::plot(result_rast) +locked_area_sf <- as.polygons(locked_area, aggregate = FALSE) |> + st_as_sf() ``` -![Map of proposed solution with both locked-in and locked-out areas.](figure/unnamed-chunk-7-1.png) - -It looks like we got the results we wanted - the locked-in areas were included in the prioritization output as an area to protect, and the locked-out areas were not included! - -# Implement locked-in areas using an `sf` input - -To generate our data for this example, we will use data from Bermuda generated via `oceandatr`. - - +Let's check our features, seamounts and locked-in area look ok: -We can generate a locked-in layer that would represent an existing MPA. - -```r -# Create a fake MPA -mpa_location <- create_area(planning_sf, 20, 50) - -# Plot fake MPA -plot(mpa_location, main = NULL, border = F) +``` r +plot(cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(locked_area_sf))) ``` -![Map of potential locked-in areas.](figure/unnamed-chunk-9-1.png) +![plot of chunk unnamed-chunk-6](figure/unnamed-chunk-6-1.png) -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). +Now we can go through the same process of data prepration using `patchwise`, then runnning a prioritization using `prioritizr` -```r -# Create seamount patches - seamount areas that touch are considered the same patch -patches_sf <- patchwise::create_patches(seamounts_sf, spatial_grid = planning_sf) +``` r +# Create seamount patches +patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_in = mpa_location) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf, locked_in = locked_area_sf) -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = locked_area_sf) -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = mpa_location) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) - -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Create the prioritization problem +problem_patches_sf <- problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_sf) |> + add_binary_decisions() |> + add_default_solver() # Solve the prioritization -solution_sf <- solve(problem_sf) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76981 rows, 56170 columns and 218036 nonzeros -#> Model fingerprint: 0xfea2ae20 -#> Variable types: 0 continuous, 56170 integer (56170 binary) +solution_patches_sf_tbl <- solve(problem_patches_sf) +#> Set parameter Username +#> Set parameter LicenseID to value 2865081 +#> Set parameter TimeLimit to value 2147483647 +#> Set parameter MIPGap to value 0.1 +#> Set parameter Presolve to value 2 +#> Set parameter Threads to value 1 +#> Academic license - for non-commercial use only - expires 2027-09-17 +#> Gurobi Optimizer version 13.0.3 build v13.0.3rc0 (linux64 - "Linux Mint 22.3") +#> +#> CPU model: Intel(R) Core(TM) i7-8700T CPU @ 2.40GHz, instruction set [SSE2|AVX|AVX2] +#> Thread count: 6 physical cores, 12 logical processors, using up to 1 threads +#> +#> Non-default parameters: +#> TimeLimit 2147483647 +#> MIPGap 0.1 +#> Presolve 2 +#> Threads 1 +#> +#> Optimize a model with 162 rows, 1603 columns and 3551 nonzeros (Min) +#> Model fingerprint: 0x87bff6b1 +#> Model has 1603 linear objective coefficients +#> Variable types: 0 continuous, 1603 integer (1603 binary) #> Coefficient statistics: -#> Matrix range [1e+00, 6e+02] -#> Objective range [2e-02, 6e+02] +#> Matrix range [1e+00, 5e+01] +#> Objective range [1e+00, 5e+01] #> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3732.8950000 -#> Presolve removed 2637 rows and 1617 columns -#> Presolve time: 0.85s -#> Presolved: 74344 rows, 54553 columns, 211669 nonzeros -#> Variable types: 0 continuous, 54553 integer (54553 binary) -#> Found heuristic solution: objective 3732.8900000 -#> Root relaxation presolve removed 334 rows and 0 columns -#> Root relaxation presolved: 74010 rows, 54553 columns, 210667 nonzeros -#> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... -#> -#> Warning: Markowitz tolerance tightened to 0.25 -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 73906 3.7232193e+03 0.000000e+00 2.534610e+03 5s -#> 84130 3.7226720e+03 0.000000e+00 1.378961e+04 10s -#> 91314 3.7213965e+03 0.000000e+00 5.732673e+03 15s -#> 96380 3.7208580e+03 0.000000e+00 9.779234e+02 20s -#> 101294 3.7206395e+03 0.000000e+00 1.848140e+03 25s -#> 107134 3.7200285e+03 0.000000e+00 1.159798e+03 30s -#> 113612 3.7199862e+03 0.000000e+00 2.998862e+03 35s -#> 117246 3.7199741e+03 0.000000e+00 1.646831e+03 40s -#> 120256 3.7199589e+03 0.000000e+00 6.978084e+02 45s -#> Concurrent spin time: 0.00s +#> RHS range [1e+00, 2e+02] #> -#> Solved with dual simplex +#> Found heuristic solution: objective 297.0000000 +#> Presolve removed 10 rows and 1438 columns +#> Presolve time: 0.00s +#> Presolved: 152 rows, 165 columns, 646 nonzeros +#> Variable types: 0 continuous, 165 integer (150 binary) +#> Root relaxation presolved: 152 rows, 165 columns, 646 nonzeros #> -#> Root simplex log... #> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 58008 3.7200450e+03 0.000000e+00 0.000000e+00 50s -#> -#> Root relaxation: objective 3.720045e+03, 58008 iterations, 49.27 seconds (144.48 work units) +#> Root relaxation: objective 2.302000e+02, 83 iterations, 0.00 seconds (0.00 work units) #> #> Nodes | Current Node | Objective Bounds | Work #> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time #> -#> 0 0 3720.04500 0 50171 3732.89000 3720.04500 0.34% - 50s +#> 0 0 230.20000 0 14 297.00000 230.20000 22.5% - 0s +#> H 0 0 239.0000000 230.20000 3.68% - 0s #> -#> Explored 1 nodes (58105 simplex iterations) in 50.55 seconds (146.73 work units) -#> Thread count was 9 (of 10 available processors) +#> Explored 1 nodes (86 simplex iterations) in 0.01 seconds (0.01 work units) +#> Thread count was 1 (of 12 available processors) #> -#> Solution count 2: 3732.89 3732.9 +#> Solution count 2: 239 297 #> #> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.732890000000e+03, best bound 3.720045002821e+03, gap 0.3441% - -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) +#> Best objective 2.390000000000e+02, best bound 2.310000000000e+02, gap 3.3473% -# Show the results -plot(result_sf, main = NULL, border = F) +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) ``` -![Map of proposed solution with locked-in areas.](figure/unnamed-chunk-10-1.png) +For comparison, we will run a prioritization without using `patchwise` -The areas in yellow are those selected by the prioritization for protection. We can clearly see the locked-in area (the yellow square) was included! -# Implement locked-out areas using an `sf` input +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The locked-in area is included in both solutions. -What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude this area from the prioritization results. +``` r +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(as.polygons(seamounts, aggregate = TRUE), col = "red")) +``` -```r -# Create a fake locked-out area -no_protection <- create_area(planning_sf, 20, 95) +![plot of chunk prioritization-plot-sf](figure/prioritization-plot-sf-1.png) -# Plot fake locked-out area -plot(no_protection, main = NULL, border = F) +# Implement locked-in and locked-out areas using an `sf` input + +As a final example we will include both a locked-in and a locked-out area. We will use the same data as before, but create a new square area for the locked-out area. + + +``` r +locked_out_area <- as.polygons(ext(25, 35, 3, 13)) |> + rasterize(pu_raster, background = 0) |> + setNames("Locked-out") |> + as.polygons(aggregate = FALSE) |> + st_as_sf() ``` -![Map of potential locked-out areas.](figure/unnamed-chunk-11-1.png) +Plot our data + -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). +``` r +plot(cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(locked_area_sf), st_drop_geometry(locked_out_area))) +``` + +![plot of chunk unnamed-chunk-8](figure/unnamed-chunk-8-1.png) -```r -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_out = no_protection) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) +``` r +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf, locked_in = locked_area_sf, locked_out = locked_out_area) -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_out = no_protection) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = locked_area_sf, locked_out = locked_out_area) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Create the prioritization problem +problem_patches_sf <- problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_sf) |> + add_binary_decisions() |> + add_default_solver() # Solve the prioritization -solution_sf <- solve(problem_sf) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76981 rows, 56170 columns and 218035 nonzeros -#> Model fingerprint: 0x908efcd6 -#> Variable types: 0 continuous, 56170 integer (56170 binary) +solution_patches_sf_tbl <- solve(problem_patches_sf) +#> Set parameter Username +#> Set parameter LicenseID to value 2865081 +#> Set parameter TimeLimit to value 2147483647 +#> Set parameter MIPGap to value 0.1 +#> Set parameter Presolve to value 2 +#> Set parameter Threads to value 1 +#> Academic license - for non-commercial use only - expires 2027-09-17 +#> Gurobi Optimizer version 13.0.3 build v13.0.3rc0 (linux64 - "Linux Mint 22.3") +#> +#> CPU model: Intel(R) Core(TM) i7-8700T CPU @ 2.40GHz, instruction set [SSE2|AVX|AVX2] +#> Thread count: 6 physical cores, 12 logical processors, using up to 1 threads +#> +#> Non-default parameters: +#> TimeLimit 2147483647 +#> MIPGap 0.1 +#> Presolve 2 +#> Threads 1 +#> +#> Optimize a model with 163 rows, 1603 columns and 3651 nonzeros (Min) +#> Model fingerprint: 0x242e4b84 +#> Model has 1603 linear objective coefficients +#> Variable types: 0 continuous, 1603 integer (1603 binary) #> Coefficient statistics: -#> Matrix range [1e+00, 6e+02] -#> Objective range [2e-02, 6e+02] +#> Matrix range [1e+00, 5e+01] +#> Objective range [1e+00, 5e+01] #> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3732.8450000 -#> Presolve removed 2781 rows and 1608 columns -#> Presolve time: 0.64s -#> Presolved: 74200 rows, 54562 columns, 211019 nonzeros -#> Variable types: 0 continuous, 54562 integer (54562 binary) -#> Root relaxation presolved: 74200 rows, 54562 columns, 211019 nonzeros +#> RHS range [1e+00, 2e+02] #> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... +#> Found heuristic solution: objective 297.0000000 +#> Presolve removed 11 rows and 1438 columns +#> Presolve time: 0.00s +#> Presolved: 152 rows, 165 columns, 646 nonzeros +#> Variable types: 0 continuous, 165 integer (150 binary) +#> Root relaxation presolved: 152 rows, 165 columns, 646 nonzeros #> -#> Warning: Markowitz tolerance tightened to 0.5 #> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 79902 3.7235278e+03 0.000000e+00 1.105348e+03 5s -#> 88841 3.7209518e+03 0.000000e+00 3.043537e+03 10s -#> 93297 3.7199971e+03 0.000000e+00 3.033930e+03 15s -#> 97537 3.7198826e+03 0.000000e+00 6.369931e+02 20s -#> 101633 3.7198377e+03 0.000000e+00 1.484675e+04 25s -#> 105681 3.7198301e+03 0.000000e+00 1.369798e+03 30s -#> Concurrent spin time: 1.31s -#> -#> Solved with dual simplex -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 41579 3.7198669e+03 0.000000e+00 0.000000e+00 32s -#> -#> Root relaxation: objective 3.719867e+03, 41579 iterations, 30.79 seconds (90.95 work units) +#> Root relaxation: objective 2.313333e+02, 65 iterations, 0.00 seconds (0.00 work units) #> #> Nodes | Current Node | Objective Bounds | Work #> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time #> -#> 0 0 3719.86693 0 49565 3732.84500 3719.86693 0.35% - 31s +#> 0 0 231.33333 0 4 297.00000 231.33333 22.1% - 0s +#> H 0 0 232.0000000 231.33333 0.29% - 0s +#> 0 0 231.33333 0 4 232.00000 231.33333 0.29% - 0s #> -#> Explored 1 nodes (41687 simplex iterations) in 31.89 seconds (92.79 work units) -#> Thread count was 9 (of 10 available processors) +#> Explored 1 nodes (65 simplex iterations) in 0.01 seconds (0.01 work units) +#> Thread count was 1 (of 12 available processors) #> -#> Solution count 1: 3732.85 +#> Solution count 2: 232 297 #> #> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.732845000000e+03, best bound 3.719866927108e+03, gap 0.3477% +#> Best objective 2.320000000000e+02, best bound 2.320000000000e+02, gap 0.0000% -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) - -# Show the results -plot(result_sf, main = NULL, border = F) +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) ``` -![Map of proposed solution with locked-out areas.](figure/unnamed-chunk-12-1.png) - -Great! The area that we don't want to include in the solution is removed! - -# Implement locked-in and locked-out areas using an `sf` input - -We can also run a scenario in which we have both locked-in and locked-out areas. - +For comparison, we will run a prioritization without using `patchwise` -```r -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_in = mpa_location, locked_out = no_protection) -#> [1] "Processing patch 1 of 7" -#> [1] "Processing patch 2 of 7" -#> [1] "Processing patch 3 of 7" -#> [1] "Processing patch 4 of 7" -#> [1] "Processing patch 5 of 7" -#> [1] "Processing patch 6 of 7" -#> [1] "Processing patch 7 of 7" -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = mpa_location, locked_out = no_protection) +The selected planning units all exclude the locked-out area (red outlined square), but include the locked-in area (selected cells in square pattern). The solution using `patchwise` includes a whole seamount (also outlined in red) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +``` r +locked_out_area_agg <- vect(locked_out_area) |> + aggregate(by = "Locked-out", count = FALSE) |> + subset(`Locked-out` == 1, NSE = TRUE) -# Solve the prioritization -solution_sf <- solve(problem_sf) -#> Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[arm]) -#> -#> CPU model: Apple M1 Max -#> Thread count: 10 physical cores, 10 logical processors, using up to 9 threads -#> -#> Optimize a model with 76982 rows, 56170 columns and 218436 nonzeros -#> Model fingerprint: 0xf4048946 -#> Variable types: 0 continuous, 56170 integer (56170 binary) -#> Coefficient statistics: -#> Matrix range [1e+00, 6e+02] -#> Objective range [2e-02, 6e+02] -#> Bounds range [1e+00, 1e+00] -#> RHS range [1e+00, 4e+03] -#> Found heuristic solution: objective 3732.9700000 -#> Presolve removed 4318 rows and 2857 columns -#> Presolve time: 0.81s -#> Presolved: 72664 rows, 53313 columns, 206757 nonzeros -#> Variable types: 0 continuous, 53313 integer (53313 binary) -#> Found heuristic solution: objective 3732.9650000 -#> Root relaxation presolve removed 334 rows and 0 columns -#> Root relaxation presolved: 72330 rows, 53313 columns, 205755 nonzeros -#> -#> Deterministic concurrent LP optimizer: primal and dual simplex -#> Showing first log only... -#> -#> Warning: Markowitz tolerance tightened to 0.25 -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 73407 3.7242230e+03 0.000000e+00 2.207557e+03 5s -#> 82887 3.7230283e+03 0.000000e+00 1.081976e+03 10s -#> 94453 3.7219914e+03 0.000000e+00 8.368012e+04 15s -#> 104318 3.7210228e+03 0.000000e+00 2.816210e+04 20s -#> 108720 3.7204223e+03 0.000000e+00 8.221710e+03 25s -#> 113456 3.7203600e+03 0.000000e+00 1.704823e+04 30s -#> 119904 3.7201420e+03 0.000000e+00 1.317371e+04 35s -#> 125763 3.7201045e+03 0.000000e+00 1.655345e+03 40s -#> 128133 3.7200997e+03 0.000000e+00 2.980871e+03 45s -#> 130433 3.7200965e+03 0.000000e+00 8.521033e+02 50s -#> Concurrent spin time: 1.78s -#> -#> Solved with dual simplex -#> -#> Root simplex log... -#> -#> Iteration Objective Primal Inf. Dual Inf. Time -#> 60783 3.7201850e+03 0.000000e+00 0.000000e+00 51s -#> -#> Root relaxation: objective 3.720185e+03, 60783 iterations, 49.83 seconds (137.67 work units) -#> -#> Nodes | Current Node | Objective Bounds | Work -#> Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time -#> -#> 0 0 3720.18505 0 48931 3732.96500 3720.18505 0.34% - 51s -#> -#> Explored 1 nodes (60901 simplex iterations) in 51.06 seconds (139.89 work units) -#> Thread count was 9 (of 10 available processors) -#> -#> Solution count 2: 3732.97 3732.97 -#> -#> Optimal solution found (tolerance 1.00e-01) -#> Best objective 3.732965000000e+03, best bound 3.720185046242e+03, gap 0.3424% - -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) - -# Show the results -plot(result_sf, main = NULL, border = F) +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(rbind(as.polygons(seamounts, aggregate = TRUE), locked_out_area_agg), col = "red")) ``` -![Map of proposed solution with locked-in and locked-out areas.](figure/unnamed-chunk-13-1.png) - -It looks like we got the results we wanted - the locked-in areas were included in the prioritization output as an area to protect, and the locked-out areas were not included! +![plot of chunk prioritization-plot-sf2](figure/prioritization-plot-sf2-1.png) diff --git a/vignettes/locked_in_out.Rmd.orig b/vignettes/locked_in_out.Rmd.orig index 21dda92..323218c 100644 --- a/vignettes/locked_in_out.Rmd.orig +++ b/vignettes/locked_in_out.Rmd.orig @@ -34,367 +34,314 @@ A quick note about `prioritizr` - sometimes there are several cells that are equ To generate our data for this example, we will use data from Bermuda generated via `oceandatr`. -```{r} +```{r message=FALSE} # Load libraries library(patchwise) +library(terra) +library(sf) +library(prioritizr) +``` -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <- '+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' +We will import some basic planning data, same as in the `patchwise` readme. For this example we will use a a 40 x 40 raster planning grid which has a cost value of 1 for each cell, and the following features that will be targeted in the prioritization: -# Create a planning grid -planning_rast <- spatialgridr::get_grid(area, projection = projection) +- 4 random, binary (0 or 1) data layers, representing species distributions; for this example we are calling them fish +- 3 "patches" of contiguous features; for this example we are calling them seamounts -# Grab all relevant data -features_rast <- oceandatr::get_features(spatial_grid = planning_rast) +```{r data-setup} +#import planning units/ cost raster +pu_raster <- rast(system.file("extdata/pu_raster.tif", package = "patchwise")) -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_rast <- stats::setNames(planning_rast, "cost") +#import fish distributions +fish_distributions <- rast(system.file("extdata/spp_distributions.tif", package = "patchwise")) -# Separate seamount data - we want to protect entire patches -seamounts_rast <- features_rast[["seamounts"]] -features_rast <- features_rast[[names(features_rast)[names(features_rast) != "seamounts"]]] +#import seamounts +seamounts <- rast(system.file("extdata/seamounts.tif", package = "patchwise")) ``` -We can generate a locked-in layer that would represent an existing MPA. +We will generate a square within our planning area that will be used as our locked in area ```{r} -# Create a function to produce a square area -create_area <- function(spatial_grid, mpa_size, edge_offset) { - if(class(spatial_grid)[1] %in% c("RasterLayer", "SpatRaster")) { - terra::rast(spatial_grid, - vals = c(rep(NA, (terra::ncol(spatial_grid)*edge_offset)), - rep(c(rep(NA, edge_offset), - rep(1, mpa_size), - rep(NA, (terra::ncol(spatial_grid)-(edge_offset+mpa_size)))), - mpa_size), - rep(NA, (terra::ncol(spatial_grid)*(terra::nrow(spatial_grid)-edge_offset-mpa_size))))) - } else if(class(spatial_grid)[1] == "sf") { - sf_to_grid <- sf::st_make_grid(spatial_grid, cellsize = sqrt(sf::st_area(spatial_grid[1,]))) - sf_centroid <- sf::st_centroid(sf_to_grid) %>% - sf::st_coordinates(.) %>% - as.data.frame() %>% - dplyr::mutate(X = round(X, digits = 4), - Y = round(Y, digits = 4)) - - # sf::st_make_grid changes the orientation of how cells are read - change them back - sf_to_grid <- sf_to_grid %>% - as.data.frame() %>% - dplyr::bind_cols(sf_centroid) %>% - dplyr::arrange(dplyr::desc(Y), X) %>% - dplyr::select(-X, -Y) - - cols_sf_to_grid <- length(unique(sf_centroid$X)) - rows_sf_to_grid <- length(unique(sf_centroid$Y)) - - sf_to_grid <- sf_to_grid %>% - sf::st_sf() %>% - dplyr::mutate(layer = c(rep(NA, (cols_sf_to_grid*edge_offset)), - rep(c(rep(NA, edge_offset), - rep(1, mpa_size), - rep(NA, (cols_sf_to_grid-(edge_offset+mpa_size)))), - mpa_size), - rep(NA, (cols_sf_to_grid*(rows_sf_to_grid-edge_offset-mpa_size))))) - - # the planning sf only includes polygons for the planning area (we need to remove them) so the number of cells matches the planning grid - sf_to_grid <- sf_to_grid[sf_to_grid$geometry %in% spatial_grid[[1]], ] - rownames(sf_to_grid) <- NULL - sf_to_grid %>% - dplyr::relocate(layer, 1) - } -} +locked_area <- as.polygons(ext(5, 15, 25, 35)) |> + rasterize(pu_raster, background = 0) |> + setNames("Locked-in") ``` -```{r fig.cap = "Map of potential locked-in areas."} -# Create a fake MPA -mpa_location <- create_area(planning_rast, 20, 50) +Let's look at our fish distributions, seamounts and locked-in area: -# Plot fake MPA -terra::plot(mpa_location) +```{r visualize-data} +plot(c(fish_distributions, seamounts, locked_area)) ``` -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). - -```{r fig.cap = "Map of proposed solution with locked-in areas."} -# Create seamount patches - seamount areas that touch are considered the same patch -patches_rast <- patchwise::create_patches(seamounts_rast) - -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_in = mpa_location) - -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) - -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_in = mpa_location) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) +Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts) - this is the same process as shown in the main `patchwise` readme, but with the locked-in area. -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +```{r message=FALSE} +# Create seamount patches +patches_rast <- create_patches(seamounts) -# Solve the prioritization -solution_rast <- solve(problem_rast) +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_rast <- create_patch_df(spatial_grid = pu_raster, features = fish_distributions, patches = patches_rast, costs = pu_raster, locked_in = locked_area) -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_rast <- features_targets(targets = rep(0.2, (nlyr(fish_distributions) + 1)), features = fish_distributions, pre_patches = seamounts, locked_in = locked_area) -# Show the results -terra::plot(result_rast) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) ``` -The areas in green are those selected by the prioritization for protection. We can clearly see the locked-in area (the green square) was included! - -# Implement locked-out areas using a `raster` input +With all the data preparation now done, we can run a prioritization using `prioritizr`: -What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude this area from the prioritization results. +```{r prioritization, include=FALSE} +# Create the prioritization problem +problem_patches <- problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_rast) |> + add_binary_decisions() |> + add_default_solver() -```{r fig.cap = "Map of potential locked-out areas."} -# Create a fake locked-out area -no_protection <- create_area(planning_rast, 20, 95) +# Solve the prioritization +solution_patches_tbl <- solve(problem_patches) -# Plot fake locked-out area -terra::plot(no_protection) +# Convert the solution into a raster using the patchwise function `convert_solution()` +sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df_rast, spatial_grid = pu_raster) |> + setNames("With patchwise") ``` -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). - -```{r fig.cap = "Map of proposed solution with locked-out areas."} -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_out = no_protection) - -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +For comparison, we will run a prioritization without using `patchwise` -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_out = no_protection) +```{r prioritization-no-patches, include=FALSE} +features_rast <- c(fish_distributions, seamounts) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) - -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +problem_no_patches <- problem(x = pu_raster, features = features_rast) |> + add_locked_in_constraints(locked_area) |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() # Solve the prioritization -solution_rast <- solve(problem_rast) +solution_no_patches <- solve(problem_no_patches) |> + setNames("Without patchwise") +``` -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The square, locked-in area is included in both scenarios. -# Show the results -terra::plot(result_rast) +```{r prioritization-plot} +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(as.polygons(seamounts), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) ``` -Great! The area that we don't want to include in the solution is removed! - -# Implement locked-in and locked-out areas using a `raster` input +# Implement locked-out areas using a `raster` input -We can also run a scenario in which we have both locked-in and locked-out areas. +What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude the "locked_area" from the prioritization results. -```{r fig.cap = "Map of proposed solution with both locked-in and locked-out areas."} -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_rast <- patchwise::create_patch_df(spatial_grid = planning_rast, features = features_rast, patches = patches_rast, costs = cost_rast, locked_in = mpa_location, locked_out = no_protection) +```{r include=FALSE} +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_rast <- create_patch_df(spatial_grid = pu_raster, features = fish_distributions, patches = patches_rast, costs = pu_raster, locked_out = locked_area) -# Create boundary matrix for prioritizr -boundary_matrix_rast <- patchwise::create_boundary_matrix(spatial_grid = planning_rast, patches = patches_rast, patch_df = patches_df_rast) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_rast <- features_targets(targets = rep(0.2, (nlyr(fish_distributions) + 1)), features = fish_distributions, pre_patches = seamounts, locked_out = locked_area) -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_rast <- patchwise::features_targets(targets = rep(0.2, (terra::nlyr(features_rast) + 1)), features = features_rast, pre_patches = seamounts_rast, locked_in = mpa_location, locked_out = no_protection) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_rast <- constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_rast <- patchwise::constraints_targets(feature_targets = targets_rast, patch_df = patches_df_rast) +# Create the prioritization problem +problem_patches <- problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_rast) |> + add_binary_decisions() |> + add_default_solver() -# Run the prioritization -problem_rast <- prioritizr::problem(x = patches_df_rast, features = constraints_rast$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_rast) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_rast) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Solve the prioritization +solution_patches_tbl <- solve(problem_patches) + +# Convert the solution into a raster using the patchwise function `convert_solution()` +sol_rast_patches <- convert_solution(solution = solution_patches_tbl, patch_df = patches_df_rast, spatial_grid = pu_raster) |> + setNames("With patchwise") +problem_no_patches <- problem(x = pu_raster, features = features_rast) |> + add_locked_out_constraints(locked_area) |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() # Solve the prioritization -solution_rast <- solve(problem_rast) +solution_no_patches <- solve(problem_no_patches) |> + setNames("Without patchwise") +``` -# Convert the prioritization into a more digestible format -result_rast <- patchwise::convert_solution(solution = solution_rast, patch_df = patches_df_rast, spatial_grid = planning_rast) +We can now plot the solutions, and overlaying the outlines of the seamounts and locked-out area (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The square, locked-out area (red outline) is excluded in both scenarios. -# Show the results -terra::plot(result_rast) +```{r prioritization-plot2} +plot(c(solution_no_patches, sol_rast_patches), + fun = function()lines(rbind(as.polygons(seamounts), as.polygons(locked_area)), col = "red"), + plg = list(legend = c("Not selected", "Selected"))) ``` -It looks like we got the results we wanted - the locked-in areas were included in the prioritization output as an area to protect, and the locked-out areas were not included! - # Implement locked-in areas using an `sf` input -To generate our data for this example, we will use data from Bermuda generated via `oceandatr`. +We can repeat the examples using `sf` input data rather than raster. First we polygonize the raster inputs -```{r include=FALSE} -# Load libraries -library(patchwise) - -# Choose area of interest (Bermuda EEZ) -area <- oceandatr::get_area(area_name = "Bermuda", mregions_column = "territory1") -projection <- '+proj=laea +lon_0=-64.8108333 +lat_0=32.3571917 +datum=WGS84 +units=m +no_defs' - -# Create a planning grid -planning_sf <- spatialgridr::get_grid(area, projection = projection, option = "sf_square") +```{r} +#create the planning grid, which is also the same as the cost grid since we are using planning units all with cost = 1 +pu_sf <- as.polygons(pu_raster, aggregate = FALSE) |> + st_as_sf() -# Grab all relevant data -features_sf <- oceandatr::get_features(spatial_grid = planning_sf) +features_sf <- as.polygons(fish_distributions, aggregate = FALSE) |> + st_as_sf() -# Create a "cost" to protecting a cell - just a uniform cost for this example -cost_sf <- features_sf %>% - dplyr::mutate(cost = 1) %>% - dplyr::select(cost) +seamounts_sf <- as.polygons(seamounts, aggregate = FALSE, na.rm = FALSE) |> + st_as_sf() -# Separate seamount data - we want to protect entire patches -seamounts_sf <- features_sf %>% - dplyr::select(seamounts) +#replace NAs with zeroes +seamounts_sf[is.na(seamounts_sf$Seamounts), "Seamounts"] <- 0 -features_sf <- features_sf %>% - dplyr::select(-seamounts) +locked_area_sf <- as.polygons(locked_area, aggregate = FALSE) |> + st_as_sf() ``` -We can generate a locked-in layer that would represent an existing MPA. +Let's check our features, seamounts and locked-in area look ok: -```{r fig.cap = "Map of potential locked-in areas."} -# Create a fake MPA -mpa_location <- create_area(planning_sf, 20, 50) - -# Plot fake MPA -plot(mpa_location, main = NULL, border = F) +```{r} +plot(cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(locked_area_sf))) ``` -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). - -```{r fig.cap = "Map of proposed solution with locked-in areas."} -# Create seamount patches - seamount areas that touch are considered the same patch -patches_sf <- patchwise::create_patches(seamounts_sf, spatial_grid = planning_sf) +Now we can go through the same process of data prepration using `patchwise`, then runnning a prioritization using `prioritizr` -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_in = mpa_location) +```{r patchwise-prep-sf} +# Create seamount patches +patches_sf <- create_patches(seamounts_sf, spatial_grid = pu_sf) -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf, locked_in = locked_area_sf) -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = mpa_location) +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = locked_area_sf) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +# Create the prioritization problem +problem_patches_sf <- problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_sf) |> + add_binary_decisions() |> + add_default_solver() # Solve the prioritization -solution_sf <- solve(problem_sf) +solution_patches_sf_tbl <- solve(problem_patches_sf) -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) - -# Show the results -plot(result_sf, main = NULL, border = F) +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) ``` -The areas in yellow are those selected by the prioritization for protection. We can clearly see the locked-in area (the yellow square) was included! - - -# Implement locked-out areas using an `sf` input +For comparison, we will run a prioritization without using `patchwise` -What if we want to note of an area that we cannot consider for protection? We can include a "locked-out" area. We repeat the example above, but this time we will exclude this area from the prioritization results. +```{r prioritization-no-patches-sf, include=FALSE} +sf_planning_data <- cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(pu_sf), st_drop_geometry(locked_area_sf)>0) #prioritizr needs TRUE/ FALSE for the locked-in column -```{r fig.cap = "Map of potential locked-out areas."} -# Create a fake locked-out area -no_protection <- create_area(planning_sf, 20, 95) +problem_no_patches_sf <- problem(x = sf_planning_data, + features = names(subset(st_drop_geometry(sf_planning_data), select = -c(cost, Locked.in))), + cost_column = "cost") |> + add_locked_in_constraints("Locked.in") |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() -# Plot fake locked-out area -plot(no_protection, main = NULL, border = F) +solution_no_patches_sf <- solve(problem_no_patches_sf) ``` -Now we can continue using `patchwise` to group seamount areas so that entire seamount ranges are protected (and not portions of seamounts). +We can now plot the solutions, and overlaying the outlines of the seamounts (in red), we can see that one entire seamount is included in the solution that used `patchwise`, whereas the solution without `patchwise` selects a few planning units in each seamount. Note that in the solution with `patchwise`, planning units that overlap seamounts that are not entirely selected have also been selected to meet targets for other features (fish distributions). The locked-in area is included in both solutions. -```{r fig.cap = "Map of proposed solution with locked-out areas."} -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_out = no_protection) - -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) - -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_out = no_protection) +```{r prioritization-plot-sf} +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(as.polygons(seamounts, aggregate = TRUE), col = "red")) +``` -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +# Implement locked-in and locked-out areas using an `sf` input -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +As a final example we will include both a locked-in and a locked-out area. We will use the same data as before, but create a new square area for the locked-out area. -# Solve the prioritization -solution_sf <- solve(problem_sf) +```{r} +locked_out_area <- as.polygons(ext(25, 35, 3, 13)) |> + rasterize(pu_raster, background = 0) |> + setNames("Locked-out") |> + as.polygons(aggregate = FALSE) |> + st_as_sf() +``` -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) +Plot our data -# Show the results -plot(result_sf, main = NULL, border = F) +```{r} +plot(cbind(features_sf, st_drop_geometry(seamounts_sf), st_drop_geometry(locked_area_sf), st_drop_geometry(locked_out_area))) ``` -Great! The area that we don't want to include in the solution is removed! -# Implement locked-in and locked-out areas using an `sf` input +```{r patchwise-prep-sf2} +# Create patches dataframe - this creates constraints so that entire seamount patches are protected +patches_df_sf <- create_patch_df(spatial_grid = pu_sf, features = features_sf, patches = patches_sf, costs = pu_sf, locked_in = locked_area_sf, locked_out = locked_out_area) -We can also run a scenario in which we have both locked-in and locked-out areas. +# Create targets for protection - 20% for each feature (including 20% of whole seamounts) +targets_sf <- features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = locked_area_sf, locked_out = locked_out_area) -```{r fig.cap = "Map of proposed solution with locked-in and locked-out areas."} -# Create patches dataframe - this creates several constraints so that entire seamount units are protected together -patches_df_sf <- patchwise::create_patch_df(spatial_grid = planning_sf, features = features_sf, patches = patches_sf, costs = cost_sf, locked_in = mpa_location, locked_out = no_protection) +# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount patches +constraints_sf <- constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) -# Create boundary matrix for prioritizr -boundary_matrix_sf <- patchwise::create_boundary_matrix(spatial_grid = planning_sf, patches = patches_sf, patch_df = patches_df_sf) +# Create the prioritization problem +problem_patches_sf <- problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") |> + add_min_set_objective() |> + add_manual_targets(constraints_sf) |> + add_binary_decisions() |> + add_default_solver() -# Create targets for protection - let's just do 20% for each feature (including 20% of whole seamounts) -targets_sf <- patchwise::features_targets(targets = rep(0.2, ncol(features_sf)), features = features_sf, pre_patches = seamounts_sf, locked_in = mpa_location, locked_out = no_protection) +# Solve the prioritization +solution_patches_sf_tbl <- solve(problem_patches_sf) -# Add these targets to targets for protection for the "constraints" we introduced to protect entire seamount units -constraints_sf <- patchwise::constraints_targets(feature_targets = targets_sf, patch_df = patches_df_sf) +# Convert the solution into an sf object using the patchwise function `convert_solution()` +solution_sf_patches <- convert_solution(solution = solution_patches_sf_tbl, patch_df = patches_df_sf, spatial_grid = pu_sf) +``` -# Run the prioritization -problem_sf <- prioritizr::problem(x = patches_df_sf, features = constraints_sf$feature, cost_column = "cost") %>% - prioritizr::add_min_set_objective() %>% - prioritizr::add_manual_targets(constraints_sf) %>% - prioritizr::add_binary_decisions() %>% - prioritizr::add_boundary_penalties(penalty = 0.000002, data = boundary_matrix_sf) %>% - prioritizr::add_default_solver(gap = 0.1, threads = parallel::detectCores()-1) +For comparison, we will run a prioritization without using `patchwise` + +```{r prioritization-no-patches-sf2, include=FALSE} +sf_planning_data <- cbind(features_sf, + st_drop_geometry(seamounts_sf), + st_drop_geometry(pu_sf), + st_drop_geometry(locked_area_sf)>0, + st_drop_geometry(locked_out_area)>0) #prioritizr needs TRUE/ FALSE for the locked-in column + +problem_no_patches_sf <- problem(x = sf_planning_data, + features = names(subset(st_drop_geometry(sf_planning_data), select = -c(cost, Locked.in, Locked.out))), + cost_column = "cost") |> + add_locked_in_constraints("Locked.in") |> + add_locked_out_constraints("Locked.out") |> + add_min_set_objective() |> + add_relative_targets(0.2) |> + add_binary_decisions() |> + add_default_solver() + +solution_no_patches_sf <- solve(problem_no_patches_sf) +``` -# Solve the prioritization -solution_sf <- solve(problem_sf) +The selected planning units all exclude the locked-out area (red outlined square), but include the locked-in area (selected cells in square pattern). The solution using `patchwise` includes a whole seamount (also outlined in red) -# Convert the prioritization into a more digestible format -result_sf <- patchwise::convert_solution(solution = solution_sf, patch_df = patches_df_sf, spatial_grid = planning_sf) +```{r prioritization-plot-sf2} +locked_out_area_agg <- vect(locked_out_area) |> + aggregate(by = "Locked-out", count = FALSE) |> + subset(`Locked-out` == 1, NSE = TRUE) -# Show the results -plot(result_sf, main = NULL, border = F) +cbind(solution_no_patches_sf[,"solution_1"], st_drop_geometry(solution_sf_patches[, "protected"])) |> + setNames(c("Without patchwise", "With patchwise", "geometry")) |> + vect() |> + plot(1:2, + type = "interval", + plg = list(legend = c("Not selected", "Selected")), + fun = function()lines(rbind(as.polygons(seamounts, aggregate = TRUE), locked_out_area_agg), col = "red")) ``` - -It looks like we got the results we wanted - the locked-in areas were included in the prioritization output as an area to protect, and the locked-out areas were not included!