Calc NLCD ratio

DESCRIPTION

@@ -8,7 +8,7 @@ Authors@R:
              email = "kyle.messier@nih.gov"))
 Description: Near Real Time air pollution model results and code produced by the SET group. It is fully tested, versioned, and open source and open access.
 Depends: R (>= 4.1.0)
-Imports: dplyr, sf, stats, terra, methods, BART, data.table, spData
+Imports: dplyr, sf, stats, terra, methods, BART, data.table, spData, exactextractr, ggplot2, scatterpie
 Suggests: covr, withr, knitr, rmarkdown, testthat (>= 3.0.0), stars, sftime
 Encoding: UTF-8
 VignetteBuilder: knitr, rmarkdown

R/calc_nlcd_ratio.R

@@ -0,0 +1,86 @@
+# needs terra, exactextractr, spData, sf
+
+globalVariables(c("us_states", "%>%"))
+
+
+#' Compute land cover classes ratio in circle buffers around points
+#'
+#' @param data_vect terra::SpatVector of points geometry
+#' @param buf_radius numeric (non-negative) giving the
+#' radius of buffer around points
+#' @param year numeric giving the year of NLCD data used
+#' @param nlcd_path character giving nlcd data path
+#' @export
+calc_nlcd_ratio <- function(data_vect,
+                            buf_radius = 1000,
+                            year = 2021,
+                            nlcd_path) {
+  # check inputs
+  if (!is.numeric(buf_radius)) {
+    stop("buf_radius is not a numeric.")
+  }
+  if (buf_radius <= 0) {
+    stop("buf_radius has not a likely value.")
+  }
+  if (!is.numeric(year)) {
+    stop("year is not a numeric.")
+  }
+  if (class(data_vect)[1] != "SpatVector") {
+    stop("data_vect is not a terra::SpatVector.")
+  }
+  if (!is.character(nlcd_path)) {
+    stop("nlcd_path is not a character.")
+  }
+  if (!file.exists(nlcd_path)) {
+    stop("nlcd_path does not exist.")
+  }
+  # open nlcd file corresponding to the year
+  nlcd_file <- list.files(nlcd_path,
+                          pattern = paste0("nlcd_", year, "_.*.tif$"),
+                          full.names = TRUE)
+  if (length(nlcd_file) == 0) {
+    stop("NLCD data not available for this year.")
+  }
+  nlcd <- rast(nlcd_file)
+  # select points within mainland US and reproject on nlcd crs if necessary
+  # need spData library
+  data(us_states, package = "spData")
+  us_main <- sf::st_union(us_states) %>%
+    terra::vect() %>%
+    terra::project(y = crs(data_vect))
+  data_vect_b <- data_vect %>%
+    terra::intersect(x = us_main)
+  if (!terra::same.crs(data_vect_b, nlcd)) {
+    data_vect_b <- terra::project(data_vect_b, crs(nlcd))
+  }
+  # create circle buffers with buf_radius
+  bufs_pol <- terra::buffer(data_vect_b, width = buf_radius) %>%
+    sf::st_as_sf()
+  # ratio of each nlcd class per buffer
+  nlcd_at_bufs <- exactextractr::exact_extract(nlcd,
+                                               st_geometry(bufs_pol),
+                                               fun = "frac",
+                                               stack_apply = TRUE,
+                                               progress = FALSE)
+  # select only the columns of interest
+  nlcd_at_bufs <- nlcd_at_bufs[names(nlcd_at_bufs)[grepl("frac_",
+                                                         names(nlcd_at_bufs))]]
+  # change column names
+  fpath <- system.file("extdata", "nlcd_classes.csv", package = "NRTAPmodel")
+  nlcd_classes <- read.csv(fpath)
+  nlcd_names <- names(nlcd_at_bufs) %>%
+    sub(pattern = "frac_", replacement = "") %>%
+    as.numeric()
+  nlcd_names <- nlcd_classes[nlcd_classes$value %in% nlcd_names, c("class")]
+  new_names <- sapply(
+    nlcd_names,
+    function(x) {
+      paste0("frac_", x, "_", year, "_", buf_radius, "m")
+    }
+  )
+  names(nlcd_at_bufs) <- new_names
+  # merge data_vect with nlcd class fractions (and reproject)
+  new_data_vect <- cbind(data_vect_b, nlcd_at_bufs) %>%
+    terra::project(crs(data_vect))
+  return(new_data_vect)
+}

inst/extdata/nlcd_classes.csv

@@ -0,0 +1,17 @@
+"","value","class","names","col"
+"1",0,"Unc","Unclassified","white"
+"2",11,"WTR","Open Water","#476ba1"
+"3",21,"OSD","Developed, Open Space","#decaca"
+"4",22,"LID","Developed, Low Intensity","#d99482"
+"5",23,"MID","Developed, Medium Intensity","#ee0000"
+"6",24,"HID","Developed, High Intensity","#ab0000"
+"7",31,"BRN","Barren Land","#b3aea3"
+"8",41,"DFO","Deciduous Forest","#68ab63"
+"9",42,"EFO","Evergreen Forest","#1c6330"
+"10",43,"MFO","Mixed Forest","#b5ca8f"
+"11",52,"SHB","Shrub/Scrub","#ccba7d"
+"12",71,"GRS","Herbaceous","#e3e3c2"
+"13",81,"PAS","Hay/Pasture","#dcd93d"
+"14",82,"CRP","Cultivated Crops","#ab7028"
+"15",90,"WDW","Woody Wetlands","#bad9eb"
+"16",95,"EHW","Emergent Herbaceous Wetlands","#70a3ba"