Merge branch 'main' of https://github.com/tobiste/tectonicr into main

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: tectonicr
 Title: Analyzing the Orientation of Maximum Horizontal Stress
-Version: 0.2.98
+Version: 0.3.0
 Authors@R: 
     person("Tobias", "Stephan", , "tobias.stephan1@yahoo.com", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-9290-014X"))

NAMESPACE

@@ -65,6 +65,7 @@ export(mean_resultant_length)
 export(model_shmax)
 export(norm_chisq)
 export(orthodrome)
+export(parse_wsm_quality)
 export(prd_err)
 export(projected_pb_strike)
 export(pvm)

NEWS.md

@@ -1,3 +1,5 @@
+# tectonicr 0.3.0
+
 # tectonicr 0.2.98 _2024-04-07_
 
 * minor fixes

R/interpolation.R

@@ -146,11 +146,11 @@ stress2grid <- function(x,
 
   # WSM method weighting (from 0 to 5)
   if (method_weighting && "type" %in% colnames_x) {
-    method_weights <- data.frame(
-      type = c("FMS", "FMF", "BO", "DIF", "HF", "GF", "GV", "OC", NA),
-      w_method = c(4, 5, 5, 5, 4, 5, 4, 2, 1) / 5
+    parse_method <- setNames(
+      c(4, 5, 5, 5, 4, 5, 4, 2, 1) / 5,
+      c("FMS", "FMF", "BO", "DIF", "HF", "GF", "GV", "OC", NA)
     )
-    x <- dplyr::left_join(x, method_weights)
+    w_method <- parse_method[x$type]
   } else {
     w_method <- rep(1, length(azi))
   }
@@ -161,8 +161,7 @@ stress2grid <- function(x,
     rep(1, length(azi))
   }
 
-  x_coords <- # sf::st_transform(x, crs = "WGS84") |>
-    sf::st_coordinates(x) |>
+  x_coords <- sf::st_coordinates(x) |>
     as.data.frame()
 
   datas <- data.frame(
@@ -176,12 +175,13 @@ stress2grid <- function(x,
   if (quality_weighting) {
     datas <- tidyr::drop_na(datas, w_quality)
   }
+  datas <- as.matrix(datas)
 
   if (is.null(grid)) {
     # Regular grid
     if (is.null(lon_range) || is.null(lat_range)) {
-      lon_range <- range(datas$lon, na.rm = TRUE)
-      lat_range <- range(datas$lat, na.rm = TRUE)
+      lon_range <- range(datas[,1], na.rm = TRUE)
+      lat_range <- range(datas[,2], na.rm = TRUE)
     }
 
     grid <- sf::st_bbox(
@@ -203,14 +203,11 @@ stress2grid <- function(x,
   stopifnot(inherits(grid, "sf"), any(sf::st_is(grid, "POINT")))
   G <- sf::st_coordinates(grid)
 
-
-  # lat.X <- lat.Y <- lon.Y <- lon.X <-
   R <- N <- numeric(nrow(G))
 
-
   SH <- c()
   for (i in seq_along(G[, 1])) {
-    distij <- dist_greatcircle(G[i, 2], G[i, 1], datas$lat, datas$lon, ...)
+    distij <- dist_greatcircle(G[i, 2], G[i, 1], datas[,2], datas[,1], ...)
 
     if (min(distij) <= arte_thres) {
       for (k in seq_along(R_range)) {
@@ -226,7 +223,7 @@ stress2grid <- function(x,
           meanSH <- mdr <- NA
         } else if (N_in_R == 1) {
           sd <- 0
-          meanSH <- datas$azi[ids_R]
+          meanSH <- datas[ids_R, 3]
           mdr <- distij[ids_R] / R_search
         } else {
           mdr <- mean(distij[ids_R], na.rm = TRUE) / R_search
@@ -240,13 +237,13 @@ stress2grid <- function(x,
             w_distance <- rep(1, length(ids_R))
           }
 
-          w <- w_distance * datas$w_quality[ids_R] * datas$w_method[ids_R]
+          w <- w_distance * datas[ids_R, 5] * datas[ids_R, 4]
 
           # mean value
           if (stat == "median") {
-            stats <- wcmedian(datas$azi[ids_R], w)
+            stats <- wcmedian(datas[ids_R, 3], w)
           } else {
-            stats <- wcmean(datas$azi[ids_R], w)
+            stats <- wcmean(datas[ids_R, 3], w)
           }
           meanSH <- as.numeric(stats[1])
           sd <- as.numeric(stats[2])

R/various.R

@@ -3,7 +3,7 @@
 #' Assigns numeric values of the precision of each measurement to the
 #' categorical quality ranking of the World Stress Map (A, B, C, D).
 #'
-#' @param regime Either a string or a character vector of WSM quality ranking
+#' @param x Either a string or a character vector of WSM quality ranking
 #'
 #' @returns \code{"integer"} or vector of type \code{"integer"}
 #'
@@ -13,17 +13,16 @@
 #' indicator. *World Stress Map Technical Report* **16-01**, GFZ German Research
 #' Centre for Geosciences. \doi{10.2312/wsm.2016.001}
 #'
-#' @export
-#'
+#' @name parse_wsm
 #' @examples
-#' quantise_wsm_quality(c("A", "B", "C", "D", NA))
+#' parse_wsm_quality(c("A", "B", "C", "D", NA))
 #' data("san_andreas")
-#' quantise_wsm_quality(san_andreas$quality)
-quantise_wsm_quality <- function(regime) {
-  as.numeric(sapply(X = regime, FUN = regime2unc))
-}
+#' parse_wsm_quality(san_andreas$quality)
+NULL
 
-regime2unc <- function(x) {
+#' @rdname parse_wsm
+#' @export
+parse_wsm_quality <- function(x) {
   c(
     "A" = 15,
     "B" = 20,
@@ -32,6 +31,12 @@ regime2unc <- function(x) {
   )[x]
 }
 
+#' @rdname parse_wsm
+#' @export
+quantise_wsm_quality <- function(x) {
+  .Deprecated(parse_wsm_quality)
+  as.numeric(sapply(X = x, FUN = parse_wsm_quality))
+}
 
 
 #' Quick analysis of a stress data set

data-raw/DATASET.R

@@ -7,7 +7,7 @@ wsm2016 <-
     regime = ifelse(regime == "TF", "T", regime),
     regime = ifelse(regime == "SS", "S", regime),
     regime = ifelse(regime == "U", NA, regime),
-    quality.quant = tectonicr::quantise_wsm_quality(quality),
+    quality.quant = tectonicr::parse_wsm_quality(quality),
     unc = ifelse(is.na(sd), quality.quant, sd),
     unc = ifelse(unc > quality.quant, quality.quant, unc),
     unc = ifelse(unc == 0, 15, unc),

tests/testthat/test_that.R

@@ -92,7 +92,7 @@ test.weights <- 1 / c(5, 1, 2, 4)
 test_that("Output of functions is as expected", {
   expect_equal(longitude_modulo(-361), -1)
   expect_equal(abs_vel(0.21, 0, r = 1), 0)
-  expect_equal(quantise_wsm_quality(c("A", "E", "F", "G", 5)), c(15, NA, NA, NA, NA))
+  expect_equal(as.numeric(parse_wsm_quality(c("A", "E", "F", "G", 5))), c(15, NA, NA, NA, NA))
   expect_equal(circular_median(c(15, 16)), 15.5)
   expect_equal(circular_median(c(15, 15, 16)), 15)
   expect_equal(circular_IQR(c(15, 16, 15, 15)), 1)

vignettes/E_interpolation.Rmd

@@ -67,7 +67,6 @@ ggplot(mean_SH) +
   ) +
   facet_wrap(~R)
 ```
-<!--  ![](interpolation.png) -->
 
 ### PoR coordinate system
 
@@ -96,7 +95,6 @@ ggplot(mean_SH_PoR) +
   ) +
   facet_wrap(~R)
 ```
-<!--  ![](interpolation_PoR.png) -->
 
 ## Rasterize the interpolation
 
@@ -134,7 +132,6 @@ ggplot(mean_SH_PoR_reduced) +
   scale_alpha("Standard deviation", range = c(1, .25)) +
   coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat))
 ```
-<!-- ![](interpolation_voronoi.png) -->
 
 The map highlights **stress anomalies** which show misfits to the direction of 
 tested plate boundary force.