From c241a915379956999b730fc826433c474bf1a09b Mon Sep 17 00:00:00 2001 From: Tobias Stephan <73840881+tobiste@users.noreply.github.com> Date: Sat, 25 May 2024 13:32:26 -0400 Subject: [PATCH] update --- NAMESPACE | 6 ++ R/spatial_interpolation.R | 195 ++++++++++++++++++++++------------- man/compact_grid.Rd | 21 ++++ man/spatial_interpolation.Rd | 15 ++- 4 files changed, 159 insertions(+), 78 deletions(-) create mode 100644 man/compact_grid.Rd diff --git a/NAMESPACE b/NAMESPACE index 85dc29c..a89d365 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/spatial_interpolation.R b/R/spatial_interpolation.R index 49d1d74..eb22e90 100644 --- a/R/spatial_interpolation.R +++ b/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() +} diff --git a/man/compact_grid.Rd b/man/compact_grid.Rd new file mode 100644 index 0000000..ed40dfe --- /dev/null +++ b/man/compact_grid.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_interpolation.R +\name{compact_grid} +\alias{compact_grid} +\title{Compact grid} +\usage{ +compact_grid(grid) +} +\arguments{ +\item{x}{output of \code{\link[=stress2grid]{stress2grid()}}, \code{\link[=PoR_stress2grid]{PoR_stress2grid()}}, or \code{\link[=kernel_dispersion]{kernel_dispersion()}}} +} +\value{ +\code{sf} object +} +\description{ +Filter spatial interpolation containing a range of search radii or kernel +half widths to find smallest wavelength (R) with the least spherical sd. +} +\seealso{ +\code{\link[=spatial_interpolation]{spatial_interpolation()}}, \code{\link[=v_mean]{v_mean()}}, \code{\link[=v_delta]{v_delta()}} +} diff --git a/man/spatial_interpolation.Rd b/man/spatial_interpolation.Rd index 7ee8f98..dff577e 100644 --- a/man/spatial_interpolation.Rd +++ b/man/spatial_interpolation.Rd @@ -17,7 +17,8 @@ spatial_interpolation( dist_weight = c("inverse", "linear"), idp = 1, dist_threshold = 0.01, - R = 10, + R_range = seq(1, 10, 1), + compact = TRUE, ... ) } @@ -44,7 +45,9 @@ spatial_interpolation( \item{dist_threshold}{Numeric. Distance weight to prevent overweight of data nearby (0 to 1). Default is \code{0.1}} -\item{R}{Numeric value or vector specifying the kernel half-width, i.e. the search radius (in km). Default is \code{1}} +\item{R_range}{Numeric value or vector specifying the kernel half-width, i.e. the search radius (in km). Default is \code{1}} + +\item{compact}{logical.} } \value{ list @@ -58,11 +61,13 @@ Based on \code{\link[tectonicr:stress2grid]{tectonicr::stress2grid()}} \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) } } \seealso{ -\code{\link[tectonicr:stress2grid]{tectonicr::stress2grid()}}, \code{\link[=v_mean]{v_mean()}}, \code{\link[=v_sd]{v_sd()}} +\code{\link[tectonicr:stress2grid]{tectonicr::stress2grid()}}, \code{\link[=v_mean]{v_mean()}}, \code{\link[=v_delta]{v_delta()}} }