Updating the strata to include the finalized BCRs including new BCR numbers

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bbsBayes2
 Type: Package
 Title: Hierarchical Bayesian Analysis of North American BBS Data
-Version: 1.2026.0
+Version: 1.2026.1
 Authors@R: c(person("Brandon P.M.", "Edwards", role = c("aut", "cre"),
                      email = "brandonedwards3@cmail.carleton.ca"),
               person("Adam C.", "Smith", role = "aut",
@@ -30,7 +30,7 @@ Imports:
     units (>= 0.8.0),
     withr (>= 2.5.0),
     cmdstanr
-Depends: R (>= 3.5)
+Depends: R (>= 4.1)
 SystemRequirements: CmdStan (https://mc-stan.org/users/interfaces/cmdstan)
 URL: https://github.com/bbsBayes/bbsBayes2,
     https://bbsbayes.github.io/bbsBayes2/
@@ -52,10 +52,8 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 LazyDataCompression: bzip2
-RoxygenNote: 7.3.3
 Roxygen: list(markdown = TRUE)
 Suggests: 
-    cmdstanr,
     covr,
     knitr,
     rmarkdown,
@@ -66,3 +64,4 @@ Suggests:
 Additional_repositories: https://mc-stan.org/r-packages/, https://ropensci.r-universe.dev
 Config/testthat/edition: 3
 VignetteBuilder: knitr
+Config/roxygen2/version: 8.0.0

R/plot-indices.R

@@ -110,7 +110,7 @@ plot_indices <- function(indices = NULL,
     to_plot <- indices[which(indices$region == i), ]
 
     samples_name <- to_plot %>%
-      dplyr::select(.data$region_type,.data$region) %>%
+      dplyr::select("region_type","region") %>%
       dplyr::distinct() %>%
       as.character() %>%
       paste(collapse = "_")

R/prepare-spatial.R

@@ -197,13 +197,17 @@ prepare_spatial <- function(prepared_data,
 
     v <- centre_union %>%
       sf::st_voronoi(envelope = bbox) %>%
-      sf::st_cast()
+      sf::st_cast() |>
+      sf::st_make_valid()
 
     vint <- sf::st_intersection(v, cov_hull_buf) %>%
       sf::st_cast("POLYGON") %>%
       sf::st_sf() %>%
       sf::st_join(centres, join = sf::st_contains) %>%
-      dplyr::arrange(.data[["strata_name"]])
+      dplyr::arrange(.data[["strata_name"]]) #%>%
+      # dplyr::group_by(.data[["strata_name"]]) %>%
+      # dplyr::summarise() %>%
+      # dplyr::arrange(.data[["strata_name"]])
 
     nb_db <- suppressWarnings(spdep::poly2nb(vint,
                             row.names = vint[["strata_name"]],

R/stratify.R

@@ -45,7 +45,7 @@
 #'   errors. For example, the map associated with the standard `bbs`
 #'   stratification excludes 3,877 surveys on 72 routes when this argument is
 #'   NULL. All of these 72 routes have starting locations on the coasts.
-#'   Setting this argument to 3000 (any route within 3 km of a polygon) ensures
+#'   Setting this argument to 4000 (any route within 4 km of a polygon) ensures
 #'   all of these coastal routes are included. Default `NULL`.
 #'   Caution: Users should be cautious of using this argument if the strata map
 #'   does not represent the full landmass of Canada and the United States. For
@@ -276,16 +276,29 @@ stratify <- function(by,
 
   if(stratify_type != "custom" & !use_map) {
     # Stratify by established groups
-    routes <- dplyr::mutate(
-      routes,
-      strata_name = dplyr::case_when(
-        stratify_by == "prov_state" ~ .data$st_abrev,
-        stratify_by == "bcr_old" ~ paste0("BCR", .data$bcr),
-        stratify_by == "latlong" ~ paste0(trunc(.data$latitude),
-                                    "_",
-                                    trunc(.data$longitude)),
-        stratify_by %in% c("bbs_usgs", "bbs_cws") ~
-          paste0(.data$country, "-", .data$st_abrev, "-", .data$bcr)))
+
+    if(stratify_by == "prov_state"){
+      routes <- dplyr::mutate(
+        routes,
+        strata_name = .data$st_abrev)
+    } else if(stratify_by == "bcr_old") {
+      routes <- dplyr::mutate(
+        routes,
+        strata_name = paste0("BCR", .data$bcr))
+    } else if(stratify_by == "latlong") {
+      routes <- dplyr::mutate(
+        routes,
+        strata_name = paste0(trunc(.data$latitude),
+                             "_",
+                             trunc(.data$longitude)))
+    } else if( stratify_by %in% c("bbs_usgs", "bbs_cws")) {
+      routes <- dplyr::mutate(
+        routes,
+        strata_name = paste0(.data$country, "-",
+                             .data$st_abrev, "-",
+                             .data$bcr))
+    }
+
 
     if(stratify_type == "subset") {
       meta_strata <- strata_custom %>%

README.md

@@ -7,12 +7,11 @@
 
 <!-- badges: end -->
 
-
-
-# bbsBayes2 1.2026.0 Released June, 2026 - Access to 2026 data release.
+# bbsBayes2 1.2026.1 Released July, 2026 - Updated BCRs and BCR names.
 
 Details on new releases below. 
 
+
 # bbsBayes2 <img src="man/figures/logo.png" align="right"/>
 
 bbsBayes2 is a package for performing hierarchical Bayesian analysis of North

data-raw/data_strata.R

@@ -25,17 +25,26 @@ library(rmapshaper) # allows simplification of the polygon boudnaries in a way
 # Updated BCR stratification ----------------------------------------------
 # updated BCR definitions that includes changes to Northern BCRS
 # circa 2025
-
-bcr_new <- "data-raw/maps_orig/bcr_2025_lakes12_statprov3.gpkg" %>%
-  read_sf() %>%
+# bcr_new <- read_sf("https://gis.birdscanada.org/server/services/Hosted/Bird_Conservation_Regions_by_Administrative_Divisions/MapServer/WFSServer?service=wfs&request=getcapabilities") |>
+# "https://gis.birdscanada.org/server/services/Hosted/Bird_Conservation_Regions_by_Administrative_Divisions/MapServer/WFSServer?service=wfs&request=getcapabilities"
+# tst <- ggplot()+
+#   geom_sf(data = tmp,
+#           aes(fill = BCR_Code))
+# tst
+#
+bcr_blank <- "data-raw/maps_orig/bcr2026_statprov.gpkg" %>%
+read_sf() %>%
   filter(!is_lake,
-         country_code %in% c("CA","US","MX"), #  removing France, and Hawaii
-         statprov_code != "HI") %>%
+         !is.na(COUNTRY),
+         COUNTRY %in% c("CAN","USA","MEX"), #  removing France, and Hawaii
+         STATEABB != "US-HI") %>%
   sf::st_transform(crs = bbsBayes2::equal_area_crs) %>%
   rmapshaper::ms_simplify(keep = 0.05,
                           keep_shapes = TRUE) %>%
-  sf::st_make_valid() %>%
-  mutate(strata_name = paste0("BCR",bcr_label))   %>%
+  sf::st_make_valid()
+
+bcr_new <- bcr_blank %>%
+  mutate(strata_name = paste0("bcr",bcr_label))   %>%
   group_by(strata_name) %>%
   summarise(.groups = "drop") %>%
   #st_cast("POLYGON") %>%
@@ -54,23 +63,30 @@ st_write(bcr_new,
 # New BCR bbs stratification ----------------------------------------------
 
 
-bbs_new <- "data-raw/maps_orig/bcr_2025_lakes12_statprov3.gpkg" %>%
-  read_sf() %>%
-  sf::st_transform(crs = bbsBayes2::equal_area_crs)
-
-bbs_new <- rmapshaper::ms_simplify(bbs_new, keep = 0.05,
-                                keep_shapes = TRUE) %>%
-mutate(strata_name = paste0(country_code,"-",statprov_code,"-",bcr_label)) %>%
-  st_make_valid() %>% # As needed when st_is_valid() fails
-  filter(!is_lake,
-         country_code %in% c("CA","US"), #  removing Mexico, France, and Hawaii
-         statprov_code != "HI") %>%
-  group_by(strata_name,bcr_label) %>%
+# bbs_new <- "data-raw/maps_orig/bcr2026_statprov.gpkg" %>%
+#   read_sf() %>%
+#   sf::st_transform(crs = bbsBayes2::equal_area_crs)
+
+bbs_new <- bcr_blank %>%
+mutate(country_code = ifelse(COUNTRY == "USA",
+                             "US",
+                             COUNTRY),
+       country_code = ifelse(COUNTRY == "CAN",
+                             "CA",
+                             country_code),
+       country_code = ifelse(COUNTRY == "MEX",
+                             "MX",
+                             country_code),
+       prov_state = str_extract(STATEABB,
+                                   "(?<=-)\\w+"),
+       strata_name = paste0(country_code,"-",prov_state,"-",bcr_label)) %>%
+  group_by(strata_name,
+           bcr_label,
+           prov_state,
+           country_code) %>%
   summarise(.groups = "drop") %>%
   mutate(area_sq_km = as.numeric(units::set_units(st_area(.), "km^2"))) %>%
-  mutate(country_code = str_extract(strata_name, "^(CA)|(US)|(MX)"),
-         bcr = bcr_label,
-         prov_state = str_extract(strata_name, "(?<=-)[A-Z]{2}"),
+  mutate(bcr = bcr_label,
          country = case_when(country_code == "CA" ~ "Canada",
                              country_code == "US" ~ "United States of America",
                              country_code == "MX" ~ "Mexico"),
@@ -191,7 +207,7 @@ st_write(strata_bcr_old,
 # updated to include the land area included in "bbs" stratification
 
 
-bbs_boundary <- bcr_new %>%
+bbs_boundary <- bcr_blank %>%
   summarise()%>%
   st_transform(crs = 4326)%>%
   st_make_valid()

inst/models/first_diff_hier_bbs_CV.stan

@@ -176,11 +176,11 @@ transformed parameters {
       }else{ // for years 2020 and 2021 (so differences from 2019-2020 and 2020-2021)
     if(y_2020[t-1]){// will only be positive in 2020 negative in 2021 allows splitting based on p_2020
     beta[,t-1] = beta_2020*p_2020  + BETA_2020*p_2020 ; // in 2020 strata-parameters forced to 0
-    YearEffect[t] = YearEffect[t-1] - BETA_2020*p_2020 ; // hyperparameter trajectory interesting to monitor but no direct inference
+    YearEffect[t] = YearEffect[t-1] + BETA_2020*p_2020 ; // hyperparameter trajectory interesting to monitor but no direct inference
 
       }else{
       beta[,t-1] = beta_2020*(1-p_2020)  + BETA_2020*(1-p_2020) ; // in 2020 strata-parameters forced to 0
-    YearEffect[t] = YearEffect[t-1] - BETA_2020*(1-p_2020) ; // hyperparameter trajectory interesting to monitor but no direct inference
+    YearEffect[t] = YearEffect[t-1] + BETA_2020*(1-p_2020) ; // hyperparameter trajectory interesting to monitor but no direct inference
 
        }
  }

inst/models/first_diff_spatial_bbs_CV.stan

@@ -177,11 +177,11 @@ transformed parameters {
       }else{ // for years 2020 and 2021 (so differences from 2019-2020 and 2020-2021)
     if(y_2020[t-1]){// will only be positive in 2020 negative in 2021 allows splitting based on p_2020
     beta[,t-1] = beta_2020*p_2020  + BETA_2020*p_2020 ; // in 2020 strata-parameters forced to 0
-    YearEffect[t] = YearEffect[t-1] - BETA_2020*p_2020 ; // hyperparameter trajectory interesting to monitor but no direct inference
+    YearEffect[t] = YearEffect[t-1] + BETA_2020*p_2020 ; // hyperparameter trajectory interesting to monitor but no direct inference
 
       }else{
       beta[,t-1] = beta_2020*(1-p_2020)  + BETA_2020*(1-p_2020) ; // in 2020 strata-parameters forced to 0
-    YearEffect[t] = YearEffect[t-1] - BETA_2020*(1-p_2020) ; // hyperparameter trajectory interesting to monitor but no direct inference
+    YearEffect[t] = YearEffect[t-1] + BETA_2020*(1-p_2020) ; // hyperparameter trajectory interesting to monitor but no direct inference
 
        }
  }