Merge pull request #26 from Spatiotemporal-Exposures-and-Toxicology/d…

…ev03

0.0.6 patch

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: scomps
 Title: Scalable R geospatial computation
-Version: 0.0.6.11212023
+Version: 0.0.6.11222023
 Authors@R: 
     person("Insang", "Song", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-8732-3256"))

NAMESPACE

@@ -23,12 +23,11 @@ export(initate_log)
 export(rast_short)
 export(set_clip_extent)
 export(sp_index_grid)
-export(sp_indexing)
 export(switch_packbound)
 export(validate_and_repair_vectors)
+import(exactextractr)
 import(future)
 import(future.apply)
-import(progressr)
 importFrom(dplyr,across)
 importFrom(methods,is)
 importFrom(rlang,sym)

R/processing.R

@@ -5,14 +5,12 @@
 #' @author Insang Song
 #' @param pnts sf or SpatVector object
 #' @param buffer_r numeric(1). buffer radius. this value will be automatically multiplied by 1.25
-#' @param nqsegs integer(1). the number of points per a quarter circle; SOON TO BE DEPRECATED
 #' @param target_input sf or SpatVector object to be clipped
 #' @return A clipped sf or SpatVector object.
 #' @export
 clip_as_extent <- function(
   pnts,
   buffer_r,
-  nqsegs = NULL,
   target_input) {
   if (any(sapply(list(pnts, buffer_r, target_input), is.null))) {
     stop("One or more required arguments are NULL. Please check.\n")
@@ -22,9 +20,7 @@ clip_as_extent <- function(
 
   if (detected_pnts != detected_target) {
     warning("Inputs are not the same class.\n")
-    target_input <- switch(detected_target,
-      sf = terra::vect(target_input),
-      terra = sf::st_as_sf(target_input))
+    target_input <- switch_packbound(target_input)
   }
 
   ext_input <- set_clip_extent(pnts, buffer_r)
@@ -200,25 +196,30 @@ extract_with_buffer_kernel <- function(
 #' 
 #' @description For simplicity, it is assumed that the coordinate systems of the points and the raster are the same. Kernel function is not yet implemented. 
 #' @param polys sf/SpatVector object. Polygons.
-#' @param surf stars/SpatRaster object. A raster of whatnot a summary will be calculated
+#' @param surf SpatRaster object. A raster from which a summary will be calculated
 #' @param id character(1). Unique identifier of each point.
-#' @param func a function taking one argument. For example, function(x) mean(x, na.rm = TRUE) or \(x) mode(x, na.rm = TRUE)
+#' @param func a generic function name in string or a function taking two arguments that are
+#'  compatible with \code{\link[exactextractr]{exact_extract}}.
+#'  For example, "mean" or or \code{\(x, w) weighted.mean(x, w, na.rm = TRUE)}
 #' @param na.rm logical(1). NA values are omitted when summary is calculated.
+#' @param grid_ref A character or sf/SpatVector object. To subset \code{polys} in \code{distribute_*} functions.
 #' @return a data.frame object with function value
 #' @author Insang Song \email{geoissong@@gmail.com}
 #' @importFrom rlang sym
 #' @importFrom dplyr across
+#' @import exactextractr
 #' @export
 extract_with_polygons <- function(
     polys,
     surf,
     id,
-    func = mean,
-    na.rm = TRUE
+    func = "mean",
+    na.rm = TRUE,
+    grid_ref = NULL
     ) {
   # type check
   stopifnot("Check class of the input points.\n" = any(methods::is(polys, "sf"), methods::is(polys, "SpatVector")))
-  stopifnot("Check class of the input raster.\n" = any(methods::is(surf, "stars"), methods::is(surf, "SpatRaster")))
+  stopifnot("Check class of the input raster.\n" = methods::is(surf, "SpatRaster"))
   stopifnot(is.character(id))
 
   cls_polys <- check_packbound(polys)
@@ -228,46 +229,25 @@ extract_with_polygons <- function(
     polys <- switch_packbound(polys)
   }
 
-  extract_with_polygons_sf <- function(polys, surf, id, func) {
-    extracted <- stars::st_extract(x = surf, at = polys, FUN = func)
-    # extracted = extracted |>
-    #   group_by(!!sym(id)) |>
-    #   summarize(across(-!!sym(id), ~func)) |>
-    #   ungroup()
-    return(extracted)
-  }
-
-  extract_with_polygons_terra <- function(polys, surf, id, func) {
-    extracted <- terra::extract(surf, polys, fun = func, ID = TRUE)
-    extracted$ID <- unlist(polys[[id]])
-    colnames(extracted)[1] <- id
-    extracted <- extracted |>
-      dplyr::group_by(!!rlang::sym(id)) |>
-      dplyr::summarize(dplyr::across(-1, ~func)) |>
-      dplyr::ungroup()
-    return(extracted)
+  if (!is.null(grid_ref)) {
+    polys <- polys[grid_ref, ]
   }
 
-  extracted_poly <- switch(
-    cls_surf,
-    sf = extract_with_polygons_sf(
-      polys = polys,
-      surf = surf,
-      id = id,
-      func = func),
-    terra = extract_with_polygons_terra(
-      polys = polys,
-      surf = surf,
-      id = id,
-      func = func)
-  )
+  extracted_poly <-
+    exactextractr::exact_extract(
+      x = surf,
+      y = sf::st_as_sf(polys),
+      fun = func,
+      force_df = TRUE,
+      append_cols = id
+    )
   return(extracted_poly)
 }
 
 
 #' Extract raster values with point buffers or polygons
 #'
-#' @param vector SpatVector object.
+#' @param vector sf/SpatVector object.
 #' @param raster SpatRaster object.
 #' @param id character(1). Unique identifier of each point.
 #' @param func function taking one numeric vector argument.
@@ -308,8 +288,17 @@ extract_with <- function(
 #'  the nearest points in threshold will be selected.
 #'  Default is \code{2 * sedc_bandwidth}.
 #' @param target_fields character(varying). Field names in characters.
-#' @note sf implementation is pending. Only available for terra.
-#' Currently the function internally converts sf objects to terra.
+#' @return data.frame (tibble) object with input field names with
+#'  a suffix \code{"_sedc"} where the sums of EDC are stored.
+#'  Additional attributes are attached for the EDC information.
+#'    - attr(result, "sedc_bandwidth"): the bandwidth where
+#'  concentration reduces to approximately five percent
+#'    - attr(result, "sedc_threshold"): the threshold distance
+#'  at which emission source points are excluded beyond that
+#' @note Distance calculation is done with terra functions internally.
+#'  Thus, the function internally converts sf objects in \code{point_*} arguments
+#'  to terra.
+#'  The optimal EDC should be carefully chosen by users.
 #' @author Insang Song
 #' @export
 calculate_sedc <-
@@ -335,8 +324,8 @@ calculate_sedc <-
     # select egrid_v only if closer than 3e5 meters from each aqs
     point_from_buf <-
       terra::buffer(
-                    point_from_buf,
-                    threshold = threshold,
+                    point_from,
+                    width = threshold,
                     quadsegs = 90)
     point_to <- point_to[point_from_buf, ]
     point_to$to_id <- len_point_to
@@ -370,8 +359,11 @@ calculate_sedc <-
                       list(sedc = ~sum(w_sedc * ., na.rm = TRUE)))
       ) |>
       dplyr::ungroup()
+
+    attr(near_from_to, "sedc_bandwidth") <- sedc_bandwidth
+    attr(near_from_to, "sedc_threshold") <- threshold
 
-    invisible(near_from_to)
+    return(near_from_to)
 }
 
 
@@ -441,14 +433,4 @@ aw_covariates <- function(
 
 }
 
-# ncbuf = terra::intersect(vect(ppb), vect(nc))
-# ncbuf_a = ncbuf
-# ncbuf_a$segarea = expanse(ncbuf_a)
-# ncbuf_k = data.frame(ncbuf_a) |>
-#   dplyr::group_by(id) |>
-#   dplyr::summarize(across(is.numeric,
-#                ~weighted.mean(., w = segarea))) |>
-#   dplyr::ungroup()
-
-#ncbufagg = terra::aggregate(ncbuf, by = 'id', fun = weighted.mean, w = ncbuf_a$segarea)
 

R/scale_process.R

@@ -91,16 +91,6 @@ distribute_process_grid <-
       future.packages = c("terra", "sf", "dplyr", "scomps", "exactextractr"))
     results_distributed <- do.call(dplyr::bind_rows, results_distributed)
 
-
-    # print(results_distributed)
-    # print(names(results_distributed))
-    # results_distributed <- results_distributed[!is.na(unlist(results_distributed[[detected_id]])),]
-
-    # post-processing
-    # names(results_distributed)[1] <- par_fun[[detected_id]]
-    # results_distributed[[par_fun[[detected_id]]]] <-
-    #   unlist(par_fun[[detected_point]][[par_fun[[detected_id]]]])
-
     return(results_distributed)
 }
 
@@ -142,13 +132,14 @@ distribute_process_grid <-
 #' # distribute_process_hierarchy()
 #' @import future
 #' @import future.apply
-#' @import progressr
 #' @export
-distribute_process_hierarchy <- function(
-  regions,
-  split_level = NULL,
-  fun_dist,
-  ...) {
+distribute_process_hierarchy <-
+  function(
+    regions,
+    split_level = NULL,
+    fun_dist,
+    ...
+  ) {
   par_fun <- list(...)
 
   if (!any(length(split_level) == 1, length(split_level) == nrow(regions))) {
@@ -228,7 +219,6 @@ distribute_process_hierarchy <- function(
 #' # distribute_process_multirasters()
 #' @import future
 #' @import future.apply
-#' @import progressr
 #' @export
 distribute_process_multirasters <- function(
   filenames,
@@ -269,12 +259,6 @@ distribute_process_multirasters <- function(
                   future.packages =
                   c("terra", "sf", "dplyr", "scomps", "future"))
   results_distributed <- do.call(dplyr::bind_rows, results_distributed)
-  # results_distributed <-
-  #   results_distributed[!is.na(results_distributed[["ID"]]), ]
-  # post-processing
-  # names(results_distributed)[1] <- par_fun[[detected_id]]
-  # results_distributed[[par_fun[[detected_id]]]] <-
-  #   unlist(par_fun[[detected_point]][[par_fun[[detected_id]]]])
 
   return(results_distributed)
 }