Fix bug in preprocess_species()

Distributions are now clipped to the clip polygon, not just its extent.

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: BirdFlowR
 Title: Predict and Visualize Bird Movement
-Version: 0.1.0.9056
+Version: 0.1.0.9060
 Authors@R: 
     c(person("Ethan", "Plunkett",  email = "plunkett@umass.edu", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0003-4405-2251")),

NEWS.md

@@ -1,8 +1,18 @@
+
+# BirdFlowR 0.1.0.9060
+2024-05-15 
+
+* Update `preprocess_species()` `crs` argument so that it accepts objects 
+of class `crs` as produced by [sf::st_crs()].
+* Fix bug in `preprocess_species()` that prevented clipping to irregular 
+  polygons.
+
+
 # BirdFlowR 0.1.0.9059
 2024-05-15
 
 * Added `callaghan_abundance` dataset on species populations.
-* Added `get_population()`  function.
+* Added `get_population()`  function.                 
 
 # BirdFlowR 0.1.0.9058
 2024-05-07

R/preprocess_species.R

@@ -328,6 +328,11 @@ preprocess_species <- function(species = NULL,
   if (is.null(crs)) {
     crs <- terra::crs(mp$custom_projection)
   } else {
+
+    # Handle objects of class crs (defined in sf)
+    if(inherits(crs, "crs"))
+      crs <- crs$wkt
+
     crs <- terra::crs(crs)
   }
 

R/process_rasters.R

@@ -165,6 +165,13 @@ process_rasters <- function(res,
   abunds_lci <- terra::project(abunds_lci, mask, method =  project_method)
   bf_msg(" done.\n")
 
+  # Clip data
+  if(!is.null(clip)){
+    abunds <- terra::mask(abunds, clip)
+    abunds_uci <- terra::mask(abunds_uci, clip)
+    abunds_lci <- terra::mask(abunds_lci, clip)
+  }
+
   # aggregate to target resolution
   if (factor != 1) {
     bf_msg("Resampling to target resolution (", res, " km)\n")

data-raw/callaghan_abundance.R

@@ -27,9 +27,17 @@ names(a) <- n
 
 # Get eBird codes from auk
 t <- auk::ebird_taxonomy
+
+# Based on scientific name
 mv <- match(a$scientific_name, t$scientific_name)
 a$species_code <- t$species_code[mv]
 
+# Based on common name
+unmatched <- is.na(a$species_code)
+
+a$species_code[unmatched] <- t$species_code[match(a$common_name[unmatched], t$common_name)]
+
+
 # Determine which ones are in the current eBird version
 ebird_ver <- ebirdst::ebirdst_version()$version_year
 r <- ebirdst::ebirdst_runs

tests/testthat/helper-get_americas.R

@@ -0,0 +1,19 @@
+
+# Function to return americas - used for setting up "big run" clipping boundary
+# here so we can reproduce the clipping issue in that run.
+get_americas <- function(clip_to_mainland_us = FALSE, include_hawaii = FALSE){
+  earth <- rnaturalearth::ne_countries(scale = 50)
+  americas <- earth[grep("America", earth$continent), , drop = FALSE]
+  if(clip_to_mainland_us){
+    extent <- c(ymax = 50, ymin = 25, xmin = -130, xmax = -55 )
+    americas <- sf::st_crop(americas, extent)
+    americas <- americas[americas$name == "United States", , drop = FALSE]
+  }
+  # Drop Hawaii
+  if(!include_hawaii){
+    clip <- sf::st_bbox(c(ymax = 25, ymin = 15, xmin = -165, xmax = -150 )) |> sf::st_as_sfc()
+    sf::st_crs(clip) <- "EPSG:4326"
+    americas <- sf::st_difference(americas, clip)
+  }
+  americas
+}

tests/testthat/test-preprocess_species.R

@@ -164,6 +164,7 @@ test_that("preprocess_species() works with clip", {
   }
 
   terra::crs(clip) <- terra::crs(proj)
+
   if (interactive()) {
     # Plot "Full" abundance for "example_data" and our clipping polygon
     a <- preprocess_species(ebirdst_example_species(), hdf5 = FALSE, res = 30)
@@ -211,3 +212,86 @@ test_that("preprocess_species() works with crs arg", {
                              crs = birdflow_crs)
   )
 })
+
+
+test_that("preprocess_species() works with clip and crs", {
+
+  # This was a problem that manifested while running lots of species with a
+  # clip polygon that buffered the Americas and while also using a custom
+  # crs.
+
+  # I'm not certain but it's possible that prior to 0.1.0.9060 clipping was
+  # only done to the rectangular extent of the clip. and certain that
+  # it failed in this case with both an irregular clip polygon and a custom
+  # crs.
+
+  # Note this test is slow and requires an active ebirdst access key
+  # See ebirdst::set_ebirdst_access_key()
+
+  skip("Slow test to documents an old bug. Always skipped.")
+
+  species <- "banswa"
+
+  # Define custom CRS
+  lat0 <- 30
+  lon0 <- -95
+  false_easting <- 5.5e6
+  false_northing <- 9.5e6
+  crs <-
+    paste0(
+      'PROJCRS["Custom_Lambert_Azimuthal",
+    BASEGEOGCRS["WGS 84",
+        DATUM["World Geodetic System 1984",
+            ELLIPSOID["WGS 84",6378137,298.257223563,
+                LENGTHUNIT["metre",1]],
+            ID["EPSG",6326]],
+        PRIMEM["Greenwich",0,
+            ANGLEUNIT["Degree",0.0174532925199433]]],
+    CONVERSION["unnamed",
+        METHOD["Lambert Azimuthal Equal Area",
+            ID["EPSG",9820]],
+        PARAMETER["Latitude of natural origin",', lat0, ',
+            ANGLEUNIT["Degree",0.0174532925199433],
+            ID["EPSG",8801]],
+        PARAMETER["Longitude of natural origin",', lon0,',
+            ANGLEUNIT["Degree",0.0174532925199433],
+            ID["EPSG",8802]],
+        PARAMETER["False easting",', false_easting,',
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8806]],
+        PARAMETER["False northing",', false_northing, ',
+            LENGTHUNIT["metre",1],
+            ID["EPSG",8807]]],
+    CS[Cartesian,2],
+        AXIS["(E)",east,
+            ORDER[1],
+            LENGTHUNIT["metre",1,
+                ID["EPSG",9001]]],
+        AXIS["(N)",north,
+            ORDER[2],
+            LENGTHUNIT["metre",1,
+                ID["EPSG",9001]]]]') |> sf::st_crs()
+
+
+  # Make Clipping polygon
+  suppressWarnings({
+    a <- get_americas(include_hawaii = FALSE) |> sf::st_transform(crs)
+    clip <- sf::st_union(a) |>  sf::st_buffer(300000)
+  })
+
+
+  bf <- preprocess_species(species, res = 200, hdf5 = FALSE, clip = clip,
+                           crs = crs, skip_quality_checks = TRUE)
+
+  r <- rasterize_distr(get_distr(bf, 20), bf)
+
+  if(FALSE)
+    plot(r)
+
+  # expect uppper right corner to be NA
+  expect_true(is.na(r[1, ncol(r)]))
+
+})
+
+
+