Merge pull request #5 from USFWS/mc-dev

Merge mc-dev with main

.gitignore

@@ -21,3 +21,14 @@ po/*~
 README.html
 DISCLAIMER.html
 NEWS.html
+# Rasters
+/docs/data/raster/nlcd/forest_distance.tif
+/docs/data/raster/nlcd/water_distance.tif
+/docs/data/raster/nlcd/nlcd.tif
+/docs/data/raster/nlcd/nlcd_ak_2016.tif
+/docs/data/raster/nlcd/nlcd.tif.aux.xml
+/docs/data/raster/nlcd/forest_distance.tif
+/docs/data/raster/ned/
+# revealjs files (many of them)
+**/docs/_mccrea_files/
+**/docs/index_files/

docs/R/create_map.R

@@ -0,0 +1,98 @@
+#' Create a leaflet map of occupancy within a refuge
+#'
+#' @param ras a `terra::ras` raster of psi estimates
+#' @param s  an `sf::st_point` of the sites surveyed
+#' @param r an `sf` multipolygon of the refuge boundary
+#' @para p whether to map the predicted value of psi ("Predicted") or SEs ("SE")
+#' @param h the height of the leaflet map returned
+#' @param w the width of the leaflet map returned
+#'
+#' @return a leaflet map
+#'
+#' @import RColorBrewer
+#' @import leaflet
+#' @import terra
+#' 
+#' @export
+#'
+#' @example 
+#' \dontrun{
+#' create_map(ras = psi, s = sites, r = tetlin, p = "Predicted", h = 650, w = 300)
+#' }
+
+create_map <- function(ras, 
+                       s, 
+                       r,
+                       p = "Predicted",
+                       h = NULL,
+                       w = NULL) {
+  if (p == "Predicted") {
+    x <- c(round(minmax(ras)[[1,1]],2), 
+           round(minmax(ras)[[2,1]],2))
+    grp <- "Psi"
+    ras <- ras$Predicted
+  } else if (p == "SE") {
+    x <- c(round(minmax(ras)[[1,2]],2), 
+           round(minmax(ras)[[2,2]],2))
+    grp <-"SE"
+    ras <- ras$SE
+  }
+
+
+  pal_rev <- colorNumeric(RColorBrewer::brewer.pal(5, "Spectral"), 
+                          x, 
+                          reverse = TRUE, 
+                          na.color = "#00000000")
+  pal <- colorNumeric(RColorBrewer::brewer.pal(5, "Spectral"), 
+                      x)
+
+  leaflet(height = h, width = w) |> 
+    addTiles() |>
+    addRasterImage(ras, 
+                   colors = pal_rev, 
+                   maxBytes = 10168580, 
+                   opacity = 0.75, 
+                   group = grp) |>
+    addCircleMarkers(data = s, lat = ~Y, lng = ~X, 
+                     radius = 0.5, 
+                     color = "black", 
+                     group = "sites") |>
+    addPolygons(data = r,
+                fill = FALSE,
+                color = "black",
+                group = "Tetlin",
+                weight = 0.5) |>
+    addLayersControl(overlayGroups = c(grp, 
+                                       "sites", 
+                                       "Tetlin"),
+                     options = layersControlOptions(collapsed = FALSE)) |>
+    addLegend(pal = pal,
+              values = x,
+              title = grp,
+              labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))) |>
+    addMiniMap(height = 100, 
+               width = 100) |>
+    addScaleBar()
+}
+
+
+
+base_map <- function(s, r, h = NULL, w = NULL) {
+  leaflet(height = h, width = w) |> 
+    addTiles() |>
+    addCircleMarkers(data = s, lat = ~Y, lng = ~X, 
+                     radius = 0.5, 
+                     color = "black", 
+                     group = "sites") |>
+    addPolygons(data = r,
+                fill = FALSE,
+                color = "black",
+                group = "Tetlin",
+                weight = 0.5) |>
+    addLayersControl(overlayGroups = c("sites", 
+                                       "Tetlin"),
+                     options = layersControlOptions(collapsed = FALSE)) |>
+    addMiniMap(height = 100, 
+               width = 100) |>
+    addScaleBar()
+}

docs/R/get_covariates.R

@@ -0,0 +1,39 @@
+## Get refuge boundary ----
+
+get_tetlin <- function(saveit = FALSE) {
+  source("./docs/R/get_sites.R")
+  tetlin <- get_refuge("Tetlin National Wildlife Refuge")
+  if (saveit) sf::st_write(tetlin, "./docs/data/shapefile/tetlin.shp")
+  tetlin
+}
+
+## Get rasters of covariates ----
+
+# Required packages
+library(FedData)
+library(terra)
+library(tidyverse)
+
+# Distance to forest cover
+nlcd <- terra::rast("./docs/data/raster/nlcd/nlcd.tif")
+forest <- terra::segregate(nlcd, classes = 42)  # Extract the forest layer
+forest <- terra::classify(forest, rcl = matrix(c(1,0,1,NA), nrow = 2, ncol = 2))  # Reclassify 0 as NA
+terra::writeRaster(forest, "./docs/data/raster/nlcd/forest.tif", overwrite = TRUE)  # Save it and load for efficiency
+forest <- terra::rast("./docs/data/raster/nlcd/forest.tif")
+forest_distance <- terra::distance(forest)
+forest_distance <- project(forest_distance, "EPSG: 4326")  # Reproject to WGS84
+forest_distance <- mask(crop(forest_distance, ext(tetlin)), tetlin)
+names(forest_distance) <- "forest"
+terra::writeRaster(forest_distance, "./docs/data/raster/nlcd/forest_distance.tif", overwrite = TRUE)  # Save it and load for efficiency
+
+# Distance to water
+nlcd <- terra::rast("./docs/data/raster/nlcd/nlcd.tif")
+water <- terra::segregate(nlcd, classes = 11)  # Extract the forest layer
+water <- terra::classify(water, rcl = matrix(c(1,0,1,NA), nrow = 2, ncol = 2))  # Reclassify 0 as NA
+terra::writeRaster(water, "./docs/data/raster/nlcd/water.tif", overwrite = TRUE)  # Save it and load for efficiency
+water <- terra::rast("./docs/data/raster/nlcd/water.tif")
+water_distance <- terra::distance(water)
+water_distance <- project(water_distance, "EPSG: 4326")  # Reproject to WGS84
+water_distance <- mask(crop(water_distance, ext(tetlin)), tetlin)
+names(water_distance) <- "water"
+terra::writeRaster(water_distance, "./docs/data/raster/nlcd/water_distance.tif", overwrite = TRUE)  # Save it and load for efficiency

docs/R/simulate_data.R

@@ -0,0 +1,101 @@
+
+
+sim_data <- function(n.occ = 8, saveit = FALSE) {  # Create a simulated unmarked data frame for a single season occupancy model
+  library(unmarked)
+  library(tidyverse)
+  library(terra)
+
+  ## Load covariate data
+  tetlin <- sf::st_read("./docs/data/shapefile/tetlin.shp", quiet = TRUE)
+  forest <- terra::rast("./docs/data/raster/nlcd/forest_distance.tif")
+  water <- terra::rast("./docs/data/raster/nlcd/water_distance.tif")
+  sites <- read.csv("./docs/data/csv/sites.csv")
+
+  set.seed(123)
+
+  # Sites
+  M <- nrow(sites)
+
+  # Sampling occasions
+  J <- n.occ
+
+  # Create a matrix of NAs
+  y <- matrix(NA, M, J)
+
+  # Create covariate data
+  site_covs <- sites %>% select(forest, water)
+
+  # Create unmarked data frame
+  umf <- unmarkedFrameOccu(y = y, siteCovs = site_covs)
+  sc <- scale(site_covs)
+  siteCovs(umf) <- sc
+
+  # Define our coefficients (detection increases with forest dist., occupancy (state) increases with water dist.)
+  cf <- list(det = c(0, 0.5), state = c(0, 0.4, -0.4))  
+
+  # Run simulation
+  unmarked_df <- simulate(umf, model = occu, formula = ~forest ~ water + forest, coefs = cf)
+
+  if (saveit) {
+    save(unmarked_df, file = "./docs/data/rdata/unmarked_df.Rdata")
+    dat_csv <- cbind(sites, unmarked_df[[1]]@y) |> 
+      rename_with(~str_c("y.", .), -(1:4))  # rename obs columns
+    write.csv(dat_csv, file = "./docs/data/csv/dat.csv")
+  }
+
+  unmarked_df
+}
+
+
+
+fit_model <- function(unmarked_df) {  # Fit a single season occupancy model
+  library(unmarked)
+  # load("./docs/data/rdata/unmarked_df.Rdata")
+  mod <- unmarked::occu(formula = ~forest ~ water + forest, 
+                        data = unmarked_df[[1]])
+}
+
+
+map_occ <- function(m = mod,  # Create a raster of species occurrence probability from an unmarked single season occupancy model
+                    pred = FALSE, 
+                    saveit = FALSE) {  
+  library(terra)
+  library(unmarked)
+  library(tidyverse)
+
+  # Load covariate rasters
+  forest <- terra::rast("./docs/data/raster/nlcd/forest_distance.tif")
+  water <- terra::rast("./docs/data/raster/nlcd/water_distance.tif")
+  sites <- read.csv("./docs/data/csv/sites.csv")
+
+  # Crop to refuge boundary
+  forest <- mask(crop(forest, ext(tetlin)), tetlin)
+  water <- mask(crop(water, ext(tetlin)), tetlin)
+
+  # Scale covariates
+  sc <- sites %>% 
+    select(forest, water) %>% 
+    scale()
+
+  # Attributes of scales dataset
+  s <- rbind(attr(sc, "scaled:center"), attr(sc, "scaled:scale"))
+
+  # Scale the raster data and combine into a single raster
+  water.s <- (water - s[[1,2]]) / s[[2,2]]
+  forest.s <- (forest - s[[1,1]]) / s[[1,2]]
+
+  if (pred) {
+    # Generate a raster of predicted values from the model
+    wf <- c(water.s, forest.s)  # Combine rasters
+    psi <- unmarked::predict(mod, type = "state", newdata = wf)
+  } else {
+    # Predict from the beta values for the state model (detection)
+    beta <- coef(mod, type = "state")
+    logit.psi <- beta[1] + beta[2]*water.s + beta[3]*forest.s
+    psi <- exp(logit.psi) / (1 + exp(logit.psi))
+  }
+
+  if (saveit) terra::writeRaster(psi, "./docs/data/raster/psi/psi.tif", overwrite = TRUE)
+
+  psi
+}

docs/R/spatial_helpers.R

@@ -0,0 +1,101 @@
+#' Get a FWS refuge boundary from the FWS AGOL feature server
+#'
+#' @param orgname the name of the refuge
+#'
+#' @return an sf multipolygon object 
+#' 
+#' @import httr
+#' @import sf
+#' 
+#' @export
+#'
+#' @example get_refuge("Tetlin National Wildlife Refuge")
+get_refuge <- function(orgname) {
+
+  orgname <- toupper(orgname)
+
+  message(paste("Downloading boundary layer for", orgname))
+
+  url <- httr::parse_url("https://services.arcgis.com/QVENGdaPbd4LUkLV/arcgis/rest/services")
+  url$path <- paste(url$path, "National_Wildlife_Refuge_System_Boundaries/FeatureServer/0/query", sep = "/")
+  url$query <- list(where = paste("ORGNAME =", paste0("'",orgname,"'")),
+                    returnGeometry = "true",
+                    f = "pgeojson"
+  )
+  request <- httr::build_url(url)
+  prop <- sf::st_read(request)
+
+  message("Done.")
+
+  return(prop)
+}
+
+
+#' Generate a random sample of points within a multipolygon and extract raster data
+#'
+#' @param refuge a sf multipolygon boundary of a refuge
+#' @param n the sample size
+#'
+#' @return
+#' 
+#' @import sf
+#' @import terra
+#' 
+#' @export
+#'
+#' @example get_sites(tetlin)
+get_sites <- function(refuge,
+                      n = 100) {
+  # Generate sample
+  sites <- sf::st_sample(refuge, n, type = "random", exact = TRUE) |>
+    sf::st_as_sf(quite = TRUE)
+  }
+
+
+extract_cov <- function(sites, forest, water) {
+  # # Import rasters
+  # forest <- terra::rast("./docs/data/raster/nlcd/forest_distance.tif")
+  # water <- terra::rast("./docs/data/raster/nlcd/water_distance.tif")
+
+  # Extract raster data at sites
+  sites <- data.frame(sf::st_coordinates(sites),
+                      forest = terra::extract(forest, sites)$forest,
+                      water = terra::extract(water, sites)$water)
+  }
+
+
+#'  Crop a spatraster (NLCD) to an sf multipoloygon (refuge boundary)
+#'
+#' @param nlcd a `spatraster` object (NLCD)
+#' @param refuge an `sf` multipolygon (refuge boundary)
+#'
+#' @return a `spatraster` cropped to the refuge boundary
+#' 
+#' @import terra
+#' @import magrittr
+#' @import sf
+#' @import dplyr
+#' @import FedData
+#' 
+#' @export
+#'
+#' @example crop_nlcd(nlcd, tetlin)
+crop_nlcd <- function(nlcd, refuge) {
+
+  # Crop raster
+  nlcd <- nlcd %>%
+    terra::crop(., sf::st_transform(refuge, sf::st_crs(terra::crs(.))),
+                snap = "out") %>%
+    terra::as.factor()
+
+  # Assign NLCD colors to the classes
+  levels(nlcd) <- FedData::nlcd_colors() %>% as.data.frame()
+  terra::coltab(nlcd) %<>% 
+    magrittr::extract2(1) %>% 
+    dplyr::filter(value %in% FedData::nlcd_colors()$ID)
+
+  return(nlcd)
+}
+
+
+

docs/_extensions/quarto-ext/fontawesome/_extension.yml

@@ -0,0 +1,7 @@
+title: Font Awesome support
+author: Carlos Scheidegger
+version: 1.1.0
+quarto-required: ">=1.2.269"
+contributes:
+  shortcodes:
+    - fontawesome.lua