update

NAMESPACE

@@ -19,6 +19,7 @@ export(as.plane)
 export(best_fit_plane)
 export(cartesian_to_acosvec)
 export(center)
+export(compact_grid)
 export(correct_pair)
 export(defgrad_from_axisangle)
 export(defgrad_from_comp)
@@ -93,15 +94,18 @@ importFrom(Directional,rfb)
 importFrom(Directional,rkent)
 importFrom(Directional,vmf.mle)
 importFrom(dplyr,arrange)
+importFrom(dplyr,as_tibble)
 importFrom(dplyr,bind_rows)
 importFrom(dplyr,filter)
 importFrom(dplyr,full_join)
+importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,near)
 importFrom(dplyr,pull)
 importFrom(dplyr,rename)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
+importFrom(dplyr,ungroup)
 importFrom(expm,expm)
 importFrom(expm,logm)
 importFrom(ggforce,geom_circle)
@@ -124,11 +128,13 @@ importFrom(sf,st_crs)
 importFrom(sf,st_is)
 importFrom(sf,st_make_grid)
 importFrom(sf,st_transform)
+importFrom(stats,aggregate)
 importFrom(tectonicr,deg2rad)
 importFrom(tectonicr,dist_greatcircle)
 importFrom(tectonicr,rad2deg)
 importFrom(tibble,as_tibble)
 importFrom(tibble,tibble)
+importFrom(tidyr,drop_na)
 importFrom(tidyr,pivot_longer)
 importFrom(tidyr,pivot_wider)
 importFrom(tidyselect,starts_with)

R/spatial_interpolation.R

@@ -1,5 +1,5 @@
 #' Spatial interpolation
-#' 
+#'
 #' Inverse distance weighted spatial interpolation of plane or line objects.
 #'
 #' @param x numeric vector, array, or object of class `"line"` or `"plane"`
@@ -13,21 +13,24 @@
 #' @param dist_weight Distance weighting method which should be used: `"linear"`, or `"inverse"` (the default).
 #' @param idp Numeric. The weighting power of inverse distance. When set to `0`, no weighting is applied.
 #' @param dist_threshold Numeric. Distance weight to prevent overweight of data nearby (0 to 1). Default is `0.1`
-#' @param R Numeric value or vector specifying the kernel half-width, i.e. the search radius (in km). Default is `1`
-#' 
+#' @param R_range Numeric value or vector specifying the kernel half-width, i.e. the search radius (in km). Default is `1`
+#' @param compact logical.
+#'
 #' @returns list
-#' @importFrom tectonicr dist_greatcircle  
+#' @importFrom tectonicr dist_greatcircle
 #' @importFrom sf st_transform st_coordinates st_is st_as_sf st_make_grid st_crs st_bbox
 #' @details Based on [tectonicr::stress2grid()]
-#' 
-#' @seealso [tectonicr::stress2grid()], [v_mean()], [v_sd()]
+#'
+#' @seealso [tectonicr::stress2grid()], [v_mean()], [v_delta()]
 #' @export
 #' @examples
 #' \dontrun{
 #' data <- read_strabo_JSON("E:/Lakehead/Field work/StraboSpot_07_02_2023.json", dataset = "TS")
-#' ps <- data$data |> dplyr::mutate(dipdir = (strike + 90) %% 360) |> dplyr::filter(type == "planar_orientation" & !(feature_type %in% c("other", "vector", "option_13")))
+#' ps <- data$data |>
+#'   dplyr::mutate(dipdir = (strike + 90) %% 360) |>
+#'   dplyr::filter(type == "planar_orientation" & !(feature_type %in% c("other", "vector", "option_13")))
 #' ps_vec <- structr::as.plane(cbind(ps$dipdir, ps$dip))
-#' structr::spatial_interpolation(x=ps_vec, coords = ps, gridsize = .05, R = 15, dist_threshold = 0.01, threshold = Inf)
+#' structr::spatial_interpolation(x = ps_vec, coords = ps, gridsize = .05, R_range = seq(1, 10, 1), dist_threshold = 0.01, threshold = Inf)
 #' }
 spatial_interpolation <- function(x,
                                   coords,
@@ -41,10 +44,11 @@ spatial_interpolation <- function(x,
                                   dist_weight = c("inverse", "linear"),
                                   idp = 1.0,
                                   dist_threshold = 0.01,
-                                  R = 10,
+                                  R_range = seq(1, 10, 1),
+                                  compact = TRUE,
                                   ...) {
   # stopifnot(inherits(coords, "sf"), is.numeric(threshold), is.numeric(arte_thres),
-  #   arte_thres > 0, is.numeric(dist_threshold), is.numeric(R), is.numeric(idp), 
+  #   arte_thres > 0, is.numeric(dist_threshold), is.numeric(R), is.numeric(idp),
   # )
 
   transform <- FALSE
@@ -59,7 +63,7 @@ spatial_interpolation <- function(x,
   dist_weight <- match.arg(dist_weight)
 
   # pre-allocating
-  #N <- nrow(x)
+  # N <- nrow(x)
   lat <- lon <- numeric()
 
   coords <- coords |>
@@ -100,80 +104,125 @@ spatial_interpolation <- function(x,
   stopifnot(inherits(grid, "sf"), any(sf::st_is(grid, "POINT")))
   G <- sf::st_coordinates(grid)
 
-  N <- numeric(nrow(G))
+  R <- N <- numeric(nrow(G))
 
   SH <- c()
   for (i in seq_along(G[, 1])) {
     distij <- tectonicr::dist_greatcircle(G[i, 2], G[i, 1], datas[, 2], datas[, 1])
 
     if (min(distij) <= arte_thres) {
-      # for (k in seq_along(R)) {
-      # R_search <- R
-      ids_R <-
-        which(distij <= R) # select those that are in search radius
-
-      N_in_R <- length(ids_R)
-
-      if (N_in_R < min_data) {
-        # not enough data within search radius
-        sd_vec <- NA
-        mean_vec <- cbind(NA, NA, NA)
-        mdr <- NA
-      } else if (N_in_R == 1) {
-        sd_vec <- 0
-        mean_vec <- datas[ids_R, 3:5]
-        mdr <- distij[ids_R] / R
-      } else {
-        mdr <- mean(distij[ids_R], na.rm = TRUE) / R
-        dist_threshold_scal <- R * dist_threshold
-
-        if (dist_weight == "linear") {
-          w <- R + 1 - max(dist_threshold_scal, distij[ids_R])
+      for (k in seq_along(R_range)) {
+        R <- R_range[k]
+        ids_R <-
+          which(distij <= R) # select those that are in search radius
+
+        N_in_R <- length(ids_R)
+
+        if (N_in_R < min_data) {
+          # not enough data within search radius
+          sd_vec <- NA
+          mean_vec <- cbind(NA, NA, NA)
+          mdr <- NA
+        } else if (N_in_R == 1) {
+          sd_vec <- 0
+          mean_vec <- datas[ids_R, 3:5]
+          mdr <- distij[ids_R] / R
         } else {
-          w <- 1 / (max(dist_threshold_scal, distij[ids_R]))^idp
+          mdr <- mean(distij[ids_R], na.rm = TRUE) / R
+          dist_threshold_scal <- R * dist_threshold
+
+          if (dist_weight == "linear") {
+            w <- R + 1 - max(dist_threshold_scal, distij[ids_R])
+          } else {
+            w <- 1 / (max(dist_threshold_scal, distij[ids_R]))^idp
+          }
+
+          # mean value
+          mean_vec <- v_mean(datas[ids_R, 3:5], w)
+          sd_vec <- v_delta(datas[ids_R, 3:5], w)
+        }
+        SH.ik <- c(
+          lon = G[i, 1],
+          lat = G[i, 2],
+          x = mean_vec[1],
+          y = mean_vec[2],
+          z = mean_vec[3],
+          delta = sd_vec / DEG2RAD(),
+          R = R,
+          mdr = mdr,
+          N = N_in_R
+        )
+
+        if (SH.ik["delta"] <= threshold & !is.na(SH.ik["delta"])) {
+          SH <- rbind(SH, SH.ik)
         }
-
-        # mean value
-        mean_vec <- v_mean(datas[ids_R, 3:5], w)
-        sd_vec <- v_delta(datas[ids_R, 3:5], w)
-      }
-      SH.ik <- c(
-        lon = G[i, 1],
-        lat = G[i, 2],
-        x = mean_vec[1],
-        y = mean_vec[2],
-        z = mean_vec[3],
-        sd = sd_vec / DEG2RAD(),
-        mdr = mdr,
-        N = N_in_R
-      )
-
-      if (SH.ik['sd'] <= threshold & !is.na(SH.ik['sd'])) {
-        SH <- rbind(SH, SH.ik)
       }
-      # }
     }
   }
 
-  res_coords <- data.frame(X = SH[, 1], Y = SH[, 2])
-
-  if (transform) {
-    vec <- to_spherical(cbind(x = SH[, 3], y = SH[, 4], z = SH[, 5]), class(x))
-  } else {
-    vec <- cbind(x = SH[, 3], y = SH[, 4], z = SH[, 5])
-  }
-
-  stats <- data.frame(
-    sd = SH[, "sd"],
-    mdr = SH[, "mdr"],
-    N = SH[, "N"]
-  )
-
-  res <- list(
-    coords = res_coords,
-    mean = vec,
-    stats = stats
-  )
+  vec <- to_spherical(cbind(x = SH[, 3], y = SH[, 4], z = SH[, 5]))
+
+  res <- dplyr::as_tibble(SH) |>
+    dplyr::rename(lon = lon.X, lat = lat.Y) |>
+    dplyr::mutate(N = as.integer(N)) |>
+    sf::st_as_sf(coords = c("lon", "lat"), crs = sf::st_crs(x), remove = FALSE) |>
+    dplyr::group_by(R)
+
+  res$dipdir <- vec[, 1]
+  res$dip <- vec[, 2]
+
+  if(compact) res <- compact_grid(res)
+
+  # res_coords <- data.frame(X = SH[, 1], Y = SH[, 2])
+  #
+  # if (transform) {
+  #   vec <- to_spherical(cbind(x = SH[, 3], y = SH[, 4], z = SH[, 5]), class(x))
+  # } else {
+  #   vec <- cbind(x = SH[, 3], y = SH[, 4], z = SH[, 5])
+  # }
+  #
+  # stats <- data.frame(
+  #   sd = SH[, "sd"],
+  #   mdr = SH[, "mdr"],
+  #   N = SH[, "N"]
+  # )
+  #
+  # res <- list(
+  #   coords = res_coords,
+  #   mean = vec,
+  #   stats = stats
+  # )
 
   return(res)
 }
+
+#' Compact grid
+#'
+#' Filter spatial interpolation containing a range of search radii or kernel
+#' half widths to find smallest wavelength (R) with the least spherical sd.
+#'
+#' @param x output of [stress2grid()], [PoR_stress2grid()], or [kernel_dispersion()]
+#' @returns \code{sf} object
+#'
+#' @importFrom dplyr ungroup mutate select left_join as_tibble
+#' @importFrom tidyr drop_na
+#' @importFrom sf st_as_sf
+#' @importFrom stats aggregate
+#' @seealso [spatial_interpolation()], [v_mean()], [v_delta()]
+#'
+#' @export
+compact_grid <- function(grid) {
+  lon <- lat <- x <- y <- z <- dipdir <- dip <- R <- numeric()
+  group <- character()
+
+  data <- grid |>
+    dplyr::ungroup() |>
+    dplyr::as_tibble() |>
+    tidyr::drop_na(x) |>
+    dplyr::mutate(group = paste(lon, lat))
+
+  aggregate(R ~ group, data, min, na.rm = TRUE) |>
+    dplyr::left_join(data, by = c("group", "R")) |>
+    dplyr::select(-group) |>
+    sf::st_as_sf()
+}