Merge pull request #18 from EcosystemEcologyLab/read-tiles-better

Deal with tiles more efficiently

R/plot_slopes.R

@@ -23,8 +23,9 @@
 #' @examples
 #' plot_slopes(tar_read(slope_liu_agb), limits = c(-1, 1), width = 5, height = 5)
 #' 
-plot_slopes <- function(slope_rast, region, target_name = NULL, save = TRUE, ext = "png", outdir = "output/figs", limits = NULL, ...) {
+plot_slopes <- function(slope_rast, region, target_name = NULL, save = TRUE, ext = "png", outdir = "output/slopes/figs", limits = NULL, ...) {
 
+  fs::dir_create(outdir)
   #varnames are lost in current version of geotargets, so pull from target name instead
   target_name <- as.character(rlang::ensym(slope_rast))
   title <- dplyr::case_when(

R/plot_summary_stats.R

@@ -14,24 +14,26 @@ plot_summary_stats <- function(data) {
       subset == "forest" ~ "Forest Service",
       subset == "wilderness" ~ "Wilderness",
       subset == "grazing" ~ "Grazing Alotments",
+      subset == "pima" ~ "Pima County",
       .default = subset
     )) |> 
     ggplot(aes(x = product, y = median, color = median > 0)) +
     facet_wrap(vars(subset), scales = "free_y") +
     geom_hline(yintercept = 0, linetype = 3) +
     geom_point() +
-    geom_linerange(aes(ymin = q.25, ymax = q.75), linewidth = 1.3) +
-    geom_linerange(aes(ymin = q.1, ymax = q.9), linewidth = 0.7) +
+    geom_linerange(aes(ymin = q25, ymax = q75), linewidth = 1.3) +
+    geom_linerange(aes(ymin = q10, ymax = q90), linewidth = 0.7) +
     labs(
       y = "Median Slope (Mg ha<sup>-1</sup> yr<sup>-1</sup>)",
       x = " ",
-      color = "Median AGB Trend"
+      color = "Median AGB Trend",
+      caption = "Error bars represent interquartile range (thick) and 10th-90th percentile (thin)."
     ) +
     scale_color_manual(values = c("TRUE"="darkgreen", "FALSE"="darkred"),
                        labels = c("TRUE" = "Increasing", "FALSE" = "Decreasing")) +
     theme_minimal() +
     theme(
       axis.text.x = element_text(angle = 45, hjust = 1),
-      axis.title = element_markdown()
+      axis.title = ggtext::element_markdown()
     ) 
 }

R/read_agb.R

@@ -8,17 +8,41 @@ read_agb <- function(path, vect) {
 
     out <- 
       terra::rast(path, win = terra::ext(vect)) |> 
-      terra::crop(vect, mask = TRUE)
+      terra::crop(vect, mask = TRUE, overwrite = TRUE)
   }
   if (fs::is_dir(path)) { #if tiles, merge them
-    #TODO apply window to these as well so unecessary tiles don't get read in
     tifs <- fs::dir_ls(path, glob = "*.tif") 
     vrt_crs <- terra::crs(terra::rast(tifs[1])) 
     vect <- vect |> terra::project(vrt_crs)
 
+    # for each tile, either drop it (if not in vect), keep it whole (if
+    # completely inside vect), or crop it to borders of vect.  In benchmarking
+    # this ended up being faster than using a vrt() and then cropping and
+    # masking.
+    # TODO: could potentially speed this up with furrr, but might be tricky with `targets`?
+    rast_list <- 
+      purrr::map(tifs, terra::rast) |> 
+      purrr::map(\(x) process_tile(x, vect)) |> 
+      purrr::compact()
+
     out <-
-      vrt(tifs, set_names = TRUE) |> 
-      crop(vect, mask = TRUE)
+      rast_list |>
+      sprc() |>
+      merge()
   }
   out
+}
+
+# for a single tile in a list of tiles
+process_tile <- function(rast, vect) {
+  #if tile is not at all within AZ, return NULL
+  if (terra::relate(ext(rast), ext(vect), relation =  "disjoint")) {
+    return(NULL)
+    #if tile is completely within AZ, return as-is
+  } else if (terra::relate(ext(rast), vect, relation = "within")) {
+    return(rast)
+    #otherwise, crop to AZ
+  } else {
+    return(crop(rast, vect, mask = TRUE, overwrite = TRUE))
+  }
 }

R/read_az_landuse.R

@@ -1,5 +1,15 @@
-#re-project and crop shapefiles to AZ
+#generic re-project and crop shapefiles to AZ
 read_az_landuse <- function(landuse_file, az) {
   terra::vect(landuse_file) |> 
-    terra::project(az) |> terra::crop(az)
-}
+    terra::project(az) |> 
+    terra::crop(az)
+}
+
+#for wilderness, subset to just national wilderness
+read_az_wilderness <- function(file, az) {
+  wilderness <- terra::vect(file)
+  wilderness[wilderness$AREATYPE == "NATIONAL WILDERNESS AREA"] |> 
+    terra::project(az) |> 
+    terra::crop(az)
+
+}

R/summarize_slopes.R

@@ -1,24 +1,26 @@
 summarize_slopes <- function(raster, sub_vect) {
   raster_name <- deparse(substitute(raster))
   sub_vect_name <- deparse(substitute(sub_vect))
-  sub_vect <- terra::project(sub_vect, raster)
-  raster <- terra::mask(raster, sub_vect)
+  raster <- raster[["slope"]] #only use coef layer, not p-value layer
+  sub_vect <- 
+    terra::project(sub_vect, raster) |> 
+    sf::st_as_sf() |> 
+    sf::st_combine() #converts to single multipolygon so only one row of summary stats per subset x product
 
-  as.data.frame(raster, na.rm = TRUE) |>
-    dplyr::summarize(
-      mean = mean(slope),
-      sd = sd(slope),
-      min = min(slope),
-      max = max(slope),
-      median = median(slope),
-      q.025 = quantile(slope, 0.025),
-      q.1 = quantile(slope, 0.1),
-      q.25 = quantile(slope, 0.25),
-      q.75 = quantile(slope, 0.75),
-      q.9 = quantile(slope, 0.9),
-      q.975 = quantile(slope, 0.975),
-      n = n()
-    ) |> 
-    mutate(subset = sub_vect_name, product = raster_name, .before = mean) 
+  exactextractr::exact_extract(
+    raster, sub_vect, fun = c(
+      "count",
+      "min",
+      "max",
+      "mean",
+      "stdev",
+      "median",
+      "quantile"
+    ),
+    quantiles = c(0.025, 0.1, 0.25, 0.75, 0.9, 0.975)
+  ) |> 
+    #fix mangled quantile column names
+    dplyr::rename(q02.5 = q02, q97.5 = q97) |> 
+    dplyr::mutate(subset = sub_vect_name, product = raster_name, .before = 1) 
 }
 

R/tar_write_csv.R

@@ -0,0 +1,6 @@
+#wrapper around readr::write_csv that returns the file path so it works with tar_file()
+tar_write_csv <- function(x, file, ...) {
+  fs::dir_create(fs::path_dir(file))
+  readr::write_csv(x, file)
+  file
+}