From f402d9a9dde8625317aa4b9ff71f08069faaf6a5 Mon Sep 17 00:00:00 2001 From: Tobias Stephan <73840881+tobiste@users.noreply.github.com> Date: Mon, 9 Dec 2024 18:14:01 -0500 Subject: [PATCH] update --- R/coordinates.R | 8 +- R/gg_stereonet.R | 154 ++++++++++++++------------- R/math.R | 177 +++++++++++++++++--------------- R/mohR.R | 2 +- README.Rmd | 4 +- README.md | 4 +- man/figures/README-stereo-1.png | Bin 12978 -> 12897 bytes man/fisher_statistics.Rd | 2 + man/ggstereo.Rd | 18 ++-- man/ggstereocontour.Rd | 30 +++--- man/or_eigen.Rd | 2 +- man/prepare-ggplot.Rd | 41 +++++--- man/stats.Rd | 6 +- vignettes/Intro.Rmd | 18 ++-- 14 files changed, 248 insertions(+), 218 deletions(-) diff --git a/R/coordinates.R b/R/coordinates.R index b8496bb..7dcc128 100644 --- a/R/coordinates.R +++ b/R/coordinates.R @@ -23,13 +23,13 @@ vec2lin0 <- function(x, y, z) { # nz <- sapply(n[, 3], function(x) ifelse(x < 0, -x, x)) nz <- n[, 3] # cbind( - azimuth = atan2d(n[, 2], n[, 1]) - plunge = asind(nz) + azimuth <- atan2d(n[, 2], n[, 1]) + plunge <- asind(nz) # ) res <- mapply(correct_inc, azi = azimuth, inc = plunge) |> t() rownames(res) <- rownames(x) - colnames(res) <- c('azimuth', 'plunge') + colnames(res) <- c("azimuth", "plunge") res } @@ -43,7 +43,7 @@ vec2fol0 <- function(x, y, z) { res <- mapply(correct_inc, azi = dip_direction, inc = dip) |> t() rownames(res) <- rownames(x) - colnames(res) <- c('dip_direction', 'dip') + colnames(res) <- c("dip_direction", "dip") res } diff --git a/R/gg_stereonet.R b/R/gg_stereonet.R index c88f5cf..91446b6 100644 --- a/R/gg_stereonet.R +++ b/R/gg_stereonet.R @@ -20,23 +20,32 @@ #' #' @examples #' if (require("mapproj")) { -#' x <- Plane(120, 85) -#' ggstereo() + -#' ggplot2::geom_point(data = gg(x), ggplot2::aes(x, y), color = "red") + -#' ggplot2::geom_path(data = ggl(x), ggplot2::aes(x, y), color = "red") +#' x <- Plane(120, 85) +#' ggstereo() + +#' ggplot2::geom_point(data = gg(x), ggplot2::aes(x, y), color = "red") + +#' ggplot2::geom_path(data = ggl(x), ggplot2::aes(x, y), color = "red") #' -#' x2 <- Line(120, 5) -#' ggstereo() + -#' ggplot2::geom_point(data = gg(x2), ggplot2::aes(x, y), color = "darkgreen") + -#' ggplot2::geom_path(data = ggl(x2, d = 8), ggplot2::aes(x, y, group = group), color = "darkgreen") -#' -#' x3 <- Plane(137, 71) -#' ggstereo() + -#' ggplot2::geom_point(data = gg(x3), ggplot2::aes(x, y), color = "darkgreen") + -#' ggplot2::geom_path(data = ggl(x3, d = 90), ggplot2::aes(x, y, group = group), color = "darkgreen", lwd = 1) + -#' ggplot2::geom_path(data = ggl(x3, d = 90 + 11), ggplot2::aes(x, y, group = group, color = 'sde <90')) + -#' ggplot2::geom_path(data = ggl(x3, d = 90 - 11), ggplot2::aes(x, y, group = group, color = 'sde >90')) -#' } +#' x2 <- Line(120, 5) +#' ggstereo() + +#' ggplot2::geom_point(data = gg(x2), ggplot2::aes(x, y), color = "darkgreen") + +#' ggplot2::geom_path(data = ggl(x2, d = 8), ggplot2::aes(x, y, group = group), color = "darkgreen") +#' +#' x3 <- Plane(137, 71) +#' ggstereo() + +#' ggplot2::geom_point(data = gg(x3), ggplot2::aes(x, y), color = "darkgreen") + +#' ggplot2::geom_path( +#' data = ggl(x3, d = 90), +#' ggplot2::aes(x, y, group = group), color = "darkgreen", lwd = 1 +#' ) + +#' ggplot2::geom_path( +#' data = ggl(x3, d = 90 + 11), +#' ggplot2::aes(x, y, group = group, color = "sde <90") +#' ) + +#' ggplot2::geom_path( +#' data = ggl(x3, d = 90 - 11), +#' ggplot2::aes(x, y, group = group, color = "sde >90") +#' ) +#' } NULL #' @rdname prepare-ggplot @@ -67,7 +76,7 @@ ggl <- function(x, ..., d = 90, n = 1e3) { stopifnot(is.spherical(x)) if (n %% 2 > 0) n <- n + 1 if (is.plane(x) | is.fault(x)) { - #x[, 1] <- 180 + x[, 1] + # x[, 1] <- 180 + x[, 1] x[, 2] <- 90 - x[, 2] } @@ -82,7 +91,7 @@ ggl <- function(x, ..., d = 90, n = 1e3) { res <- matrix(ncol = 3, nrow = n * nx) |> as.data.frame() - colnames(res) <- c("x", "y", "id") + colnames(res) <- c("x", "y", "id") for (i in seq_along(x[, 1])) { D <- cbind(azimuth = seq(0, 360, l = n), plunge = rep(90 - d[i], n)) |> @@ -95,13 +104,13 @@ ggl <- function(x, ..., d = 90, n = 1e3) { D1 <- vrotate(D, zaxis, deg2rad(strike)) rotangle <- vangle(zaxis, as.line(x[i, ])) - if(d[i] < 90 & is.plane(x)){ + if (d[i] < 90 & is.plane(x)) { d[i] <- 180 - d[i] D1 <- -D1 x[i, 1] <- x[i, 1] + 180 - x[i, 2] <- 90- x[i, 2] + x[i, 2] <- 90 - x[i, 2] } - + if (d[i] < 90) { k <- -1 } else { @@ -116,32 +125,32 @@ ggl <- function(x, ..., d = 90, n = 1e3) { if (d[i] != 90 & d[i] > as.line(x[i, ])[2]) { # upper hemisphere - #flag <- TRUE - #D_rotrot <- vrotate(D_rot, zaxis, deg2rad(180)) + # flag <- TRUE + # D_rotrot <- vrotate(D_rot, zaxis, deg2rad(180)) # D_rotrot_fixed <- list(az = D_rotrot[, 1], inc = D_rotrot[, 2]) |> # fix_inc() D_rotrot_fixed <- D_fixed - #prec <- sqrt(.Machine$double.eps) + # prec <- sqrt(.Machine$double.eps) dangle <- vangle(Line(D_fixed$az, D_fixed$inc), as.line(x[i, ])) - #cond <- d[i] >= (dangle + prec) | d[i] <= (dangle - prec) + # cond <- d[i] >= (dangle + prec) | d[i] <= (dangle - prec) cond <- dplyr::near(d[i], dangle) D_fixed$az[cond] <- D_rotrot_fixed$az[cond] D_fixed$inc[cond] <- D_rotrot_fixed$inc[cond] - } #else { - #cond <- rep(FALSE, n) - #} + } # else { + # cond <- rep(FALSE, n) + # } D_fixed2 <- data.frame(x = (180 - D_fixed$az), y = D_fixed$inc) if (d[i] < 90) { D_fixed3 <- D_fixed2 - #D_fixed3$cond <- cond + # D_fixed3$cond <- cond } else { D_fixed3 <- utils::tail(D_fixed2, n = n / 2) - #D_fixed3$cond <- FALSE + # D_fixed3$cond <- FALSE } D_fixed3$id <- i @@ -159,9 +168,9 @@ ggl <- function(x, ..., d = 90, n = 1e3) { } #' Stereoplot Perimeter -#' +#' #' Adds a frame to the stereographic projection -#' +#' #' @param n resolution of frame #' #' @param color,fill,lwd Graphical parameters @@ -173,7 +182,7 @@ ggframe <- function(n = 1e4, color = "black", fill = NA, lwd = 1, ...) { prim.l1 <- seq(0, 180, length = n / 2) prim.l2 <- seq(-180, 0, length = n / 2) prim.long <- c(prim.l1, prim.l2) - + prim_df <- data.frame(prim.long, prim.lat) geom_polygon(aes(x = prim.long, y = prim.lat), data = prim_df, color = color, fill = fill, lwd = lwd, ..., inherit.aes = FALSE) @@ -211,10 +220,10 @@ ggstereo_grid <- function(d = 10, rot = 0, ...) { #' Stereonet using ggplot #' -#' @param data Default dataset to use for plot. If not already a data.frame, -#' will be converted to one by [ggplot2::fortify()]. If not specified, must be +#' @param data Default dataset to use for plot. If not already a data.frame, +#' will be converted to one by [ggplot2::fortify()]. If not specified, must be #' supplied in each layer added to the plot. -#' @param mapping Default list of aesthetic mappings to use for plot. If not +#' @param mapping Default list of aesthetic mappings to use for plot. If not #' specified, must be supplied in each layer added to the plot. #' @param earea logical. Whether the projection is equal-area ("Schmidt net") #' (`TRUE`, the default), or equal-angle ("Wulff net") (`FALSE`). @@ -232,17 +241,17 @@ ggstereo_grid <- function(d = 10, rot = 0, ...) { #' #' @examples #' if (require("mapproj")) { -#' test_data <- rbind( -#' rvmf(100, mu = Line(90, 45), k = 10), -#' rvmf(50, mu = Line(0, 0), k = 20) -#' ) |> as.line() +#' test_data <- rbind( +#' rvmf(100, mu = Line(90, 45), k = 10), +#' rvmf(50, mu = Line(0, 0), k = 20) +#' ) |> as.line() #' -#' ggstereo(grid = TRUE) + -#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) +#' ggstereo(grid = TRUE) + +#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) #' -#' ggstereo(earea = FALSE, centercross = TRUE) + -#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) -#' } +#' ggstereo(earea = FALSE, centercross = TRUE) + +#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) +#' } ggstereo <- function(data = NULL, mapping = aes(), earea = TRUE, centercross = TRUE, grid = FALSE, grid.spacing = 10, grid.rot = 0, ...) { # if(earea){ # crs = "+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0" @@ -250,19 +259,20 @@ ggstereo <- function(data = NULL, mapping = aes(), earea = TRUE, centercross = T # crs = "+proj=stere +lat_0=90 +lon_0=0 +x_0=0 +y_0=0" # } rlang::check_installed("mapproj", reason = "to use `coord_map()`") - + ggplot(data = data, mapping = mapping) + - #theme_void() + + # theme_void() + theme( - title = element_text(face = "bold"), - panel.background = element_blank(), + plot.title = element_text(face = "bold"), + panel.background = element_blank(), panel.border = element_blank(), panel.grid = element_blank(), - axis.ticks = element_blank(), - axis.title = element_blank(), + axis.ticks = element_blank(), + axis.title = element_blank(), axis.text = element_blank(), - #legend.title = element_blank() - ) + { + # legend.title = element_blank() + ) + + { if (grid) { ggstereo_grid(d = grid.spacing, rot = grid.rot, color = "grey90", lwd = .2) } @@ -271,8 +281,8 @@ ggstereo <- function(data = NULL, mapping = aes(), earea = TRUE, centercross = T annotate("point", x = 0, y = 90, pch = as.numeric(centercross) * 3) + scale_y_continuous(limits = c(0, 90)) + scale_x_continuous(limits = c(-180, 180)) + - coord_map(ifelse(earea, "azequalarea", "stereographic"), orientation = c(90, 0, 0)) - # coord_sf(crs = crs, default_crs = crs) + coord_map(ifelse(earea, "azequalarea", "stereographic"), orientation = c(90, 0, 0)) + # coord_sf(crs = crs, default_crs = crs) } ignore_unused_imports <- function() { @@ -357,23 +367,23 @@ vmf_kerncontour <- function(u, hw = NULL, kernel_method = c("cross", "rot"), ngr #' @name ggstereocontour #' @examples #' if (require("mapproj")) { -#' test_data <- rbind( -#' rvmf(100, mu = Line(90, 45), k = 10), -#' rvmf(50, mu = Line(0, 0), k = 20) -#' ) |> as.line() -#' -#' ggstereo() + -#' geom_contourf_stereo(gg(test_data)) + -#' ggplot2::scale_fill_viridis_d(option = "A") + -#' # guides(fill = guide_colorsteps(barheight = unit(8, "cm"), show.limits = TRUE)) + -#' geom_contour_stereo(gg(test_data), color = "grey") + -#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y), color = "lightgrey") + -#' ggframe() +#' test_data <- rbind( +#' rvmf(100, mu = Line(90, 45), k = 10), +#' rvmf(50, mu = Line(0, 0), k = 20) +#' ) |> as.line() +#' +#' ggstereo() + +#' geom_contourf_stereo(gg(test_data)) + +#' ggplot2::scale_fill_viridis_d(option = "A") + +#' # guides(fill = guide_colorsteps(barheight = unit(8, "cm"), show.limits = TRUE)) + +#' geom_contour_stereo(gg(test_data), color = "grey") + +#' ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y), color = "lightgrey") + +#' ggframe() #' -#' ggstereo() + -#' geom_contourf_stereo(gg(test_data), norm = TRUE, bins = 50, threshold = .1) + -#' ggplot2::scale_fill_viridis_d(option = "A") -#' } +#' ggstereo() + +#' geom_contourf_stereo(gg(test_data), norm = TRUE, bins = 50, threshold = .1) + +#' ggplot2::scale_fill_viridis_d(option = "A") +#' } NULL #' @rdname ggstereocontour @@ -389,7 +399,7 @@ geom_contour_stereo <- function(data, ngrid = 200, hw = NULL, optimal_bw = c("cr res$Density <- normalize(res$Density) } res$Density[res$Density <= threshold] <- NA - + geom_contour(data = res, aes(x = -Long, y = Lat, z = Density), ...) } diff --git a/R/math.R b/R/math.R index 328ecac..d2a367a 100644 --- a/R/math.R +++ b/R/math.R @@ -47,15 +47,14 @@ vsum <- function(x, w = NULL) { #' @rdname operations #' @export vnorm <- function(x) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } xn <- x / vlength(x) if (transform) { - to_spherical(xn) + to_spherical(xn, class = class) } else { xn } @@ -64,11 +63,10 @@ vnorm <- function(x) { #' @rdname operations #' @export vcross <- function(x, y) { - transform <- FALSE + transform <- is.spherical(x) if (is.spherical(x)) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -98,7 +96,6 @@ vdot <- function(x, y) { if (is.spherical(x)) { class <- class(x) x <- to_vec(x) - # transform <- TRUE } else { x <- vec2mat(x) } @@ -114,11 +111,10 @@ vdot <- function(x, y) { #' @rdname operations #' @export vrotate <- function(x, rotaxis, rotangle) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -142,11 +138,10 @@ vrotate <- function(x, rotaxis, rotangle) { } vrotaxis <- function(x, y) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -177,16 +172,15 @@ vrotaxis <- function(x, y) { #' @rdname operations #' @export vangle <- function(x, y) { - transform <- FALSE - if (is.spherical(x)) { - classx <- class(x) + transform <- is.spherical(x) + if (transform) { + # classx <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } if (is.spherical(y)) { - classy <- class(y) + # classy <- class(y) y <- to_vec(y) } else { y <- vec2mat(y) @@ -204,11 +198,10 @@ vangle <- function(x, y) { #' @rdname operations #' @export vproject <- function(x, y) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -238,11 +231,10 @@ vproject_length <- function(x, y) { #' @rdname operations #' @export vreject <- function(x, y) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -275,11 +267,10 @@ v_orthogonalize <- function(x, y) { if (abs(vdot(x, y)) < 1.0e-4) { return(list(x, y)) } else { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { classx <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -353,11 +344,10 @@ v_orthogonalize <- function(x, y) { #' vtransform(vec, mat) vtransform <- function(x, A, norm = FALSE) { stopifnot(is.matrix(A)) - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -396,12 +386,24 @@ vresultant <- function(x, w = NULL, mean = FALSE) { as.numeric(w) } + transform <- is.spherical(x) + if (transform) { + v <- to_vec(x) + } else { + v <- vec2mat(x) + } + + R <- vsum(x, w) if (mean) { N <- sum(w) R <- R / N } - R + if (transform) { + to_spherical(R, class(x)) + } else { + R + } } #' Statistical estimators of the distribution of a set of vectors @@ -423,7 +425,7 @@ vresultant <- function(x, w = NULL, mean = FALSE) { #' if otherwise). For enough large sample it approaches the angular standard #' deviation (`"csd"`) of the Fisher statistics. #' -#' `v_rdegree` returns the degree of preferred orientation of vectors (in %). +#' `v_rdegree` returns the degree of preferred orientation of vectors, range: (0, 1). #' #' `v_sde` returns the spherical standard error (numeric). If the number of #' data is less than 25, if will print a additional message, that the output @@ -436,7 +438,9 @@ vresultant <- function(x, w = NULL, mean = FALSE) { #' @seealso [vresultant()], [fisher_statistics()] #' @name stats #' @examples +#' set.seed(1234) #' x <- rvmf(100, mu = Line(120, 50), k = 5) +#' v_mean(x) #' v_var(x) #' v_delta(x) #' v_rdegree(x) @@ -446,7 +450,7 @@ vresultant <- function(x, w = NULL, mean = FALSE) { #' fisher_statistics(x) #' #' #' weights: -#' x2 <- Line(c(0, 0), c(0, 90)) +#' x2 <- Line(c(0, 0), c(0, 90)) #' v_mean(x2) #' v_mean(x2, w = c(1, 2)) #' v_var(x2) @@ -456,10 +460,9 @@ NULL #' @rdname stats #' @export v_mean <- function(x, w = NULL) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { v <- to_vec(x) - transform <- TRUE } else { v <- vec2mat(x) } @@ -516,8 +519,8 @@ v_sd <- function(x, w = NULL) { #' @rdname stats #' @export v_delta <- function(x, w = NULL) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { v <- to_vec(x) } else { v <- vec2mat(x) @@ -526,7 +529,7 @@ v_delta <- function(x, w = NULL) { vlength() d <- acos(Rbar) - if (is.spherical(x)) { + if (transform) { rad2deg(d) } else { d @@ -541,18 +544,18 @@ v_rdegree <- function(x, w = NULL) { } else { v <- vec2mat(x) } - #N <- nrow(v) + w <- if (is.null(w)) { rep(1, times = nrow(v)) } else { as.numeric(w) } - + N <- sum(w) Rbar <- vresultant(vnorm(v), w, mean = FALSE) |> vlength() - 100 * (2 * Rbar - N) / N + (2 * Rbar - N) / N } #' @rdname stats @@ -568,9 +571,9 @@ v_sde <- function(x, w = NULL) { } else { as.numeric(w) } - + N <- sum(w) - + if (N < 25) warning("The standard error might not be a good estimator for N < 25") xbar <- vsum(v, w) / N Rbar <- vlength(xbar) @@ -608,16 +611,15 @@ v_confidence_angle <- function(x, w = NULL, alpha = 0.05) { } v_antipode <- function(x) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { class <- class(x) v <- to_vec(x) - transform <- TRUE } else { v <- vec2mat(x) } - #xa <- vrotate(v, c(1, 0, 0), pi) + # xa <- vrotate(v, c(1, 0, 0), pi) xa <- -v if (transform) { @@ -652,6 +654,7 @@ estimate_k <- function(x, w = NULL) { #' #' @param x numeric. Can be three element vector, three column array, or an #' object of class `"line"` or `"plane"` +#' @param w numeric. weights #' @returns list, with #' \describe{ #' \item{`"k"`}{estimated concentration parameter \eqn{\kappa} for the von Mises-Fisher @@ -664,10 +667,9 @@ estimate_k <- function(x, w = NULL) { #' x <- rvmf(100, mu = Line(120, 50), k = 5) #' fisher_statistics(x) fisher_statistics <- function(x, w = NULL) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { x <- to_vec(x) - transform <- TRUE } else { x <- vec2mat(x) } @@ -676,9 +678,9 @@ fisher_statistics <- function(x, w = NULL) { } else { as.numeric(w) } - + N <- sum(w) - + R <- vresultant(x, w) |> vlength() @@ -703,15 +705,15 @@ fisher_statistics <- function(x, w = NULL) { #' @note For non-unit vectors the interpolation is not uniform #' @details #' A Slerp path is the spherical geometry equivalent of a path along a line -#' segment in the plane; a great circle is a spherical geodesic. +#' segment in the plane; a great circle is a spherical geodesic. #' @export vslerp <- function(x, y, t) { - transform <- FALSE - if (is.spherical(x)) { + transform <- is.spherical(x) + if (transform) { + class <- class(x) x <- to_vec(x) } else { x <- vec2mat(x) - transform <- TRUE } if (is.spherical(y)) { y <- to_vec(y) @@ -722,7 +724,7 @@ vslerp <- function(x, y, t) { slerp <- (x * sin((1 - t) * theta) + y * sin(t * theta)) / sin(theta) if (transform) { - to_spherical(slerp, class(x)) + to_spherical(slerp, class) } else { slerp } @@ -785,7 +787,7 @@ ortensor <- function(x, norm = TRUE) { #' @export #' @examples #' set.seed(1) -#' mu <- rvmf(n = 1) |> vec2line() +#' mu <- rvmf(n = 1) |> vec2line() #' x <- rfb(100, mu = mu, k = 1, A = diag(c(10, 0, 0))) #' x_eigen <- or_eigen(x) #' x_eigen @@ -821,7 +823,7 @@ or_eigen <- function(x, scaled = FALSE) { #' #' @importFrom dplyr near #' @name strain_shape -#' +#' #' @details \describe{ #' \item{`stretch_ratios`}{Sqrt of eigenvalue ratios} #' \item{`strain_ratios`}{Log of stretch ratios} @@ -837,35 +839,35 @@ or_eigen <- function(x, scaled = FALSE) { #' \item{`MAD`}{maximum angular deviation (Kirschvink, 1980)} #' \item{`US`}{Uniformity statistic of Mardia (1972)} #' } -#' +#' #' @seealso [ortensor()], [or_eigen()], [fabric_indexes()] -#' +#' #' @returns list -#' -#' @references +#' +#' @references #' Flinn, Derek.(1963): "On the statistical analysis of fabric diagrams." Geological Journal 3.2: 247-253. -#' +#' #' Kirschvink, J. (1980): The least-squares line and plane and the analysis of palaeomagnetic data. Geophysical Journal International, 62(3), 699-718. -#' +#' #' Lisle, Richard J. (1985): "The use of the orientation tensor for the description and statistical testing of fabrics." Journal of Structural Geology 7.1: 115-117. -#' -#' Lode, Walter (1926): "Versuche über den Einfluß der mittleren Hauptspannung auf das Fließen der Metalle Eisen, Kupfer und Nickel“ +#' +#' Lode, Walter (1926): "Versuche über den Einfluß der mittleren Hauptspannung auf das Fließen der Metalle Eisen, Kupfer und Nickel“ #' (*"Experiments on the influence of the mean principal stress on the flow of the metals iron, copper and nickel"*], Zeitschrift für Physik, vol. 36 (November), pp. 913–939, DOI: 10.1007/BF01400222 -#' +#' #' Mardia, Kantilal Varichand. (1975): "Statistics of directional data." Journal of the Royal Statistical Society Series B: Statistical Methodology 37.3: 349-371. -#' +#' #' Nadai, A., and Hodge, P. G., Jr. (1963): "Theory of Flow and Fracture of Solids, vol. II." ASME. J. Appl. Mech. December 1963; 30(4): 640. https://doi.org/10.1115/1.3636654 -#' -#' Ramsay, John G. (1967): "Folding and fracturing of rocks." Mc Graw Hill Book Company 568. -#' +#' +#' Ramsay, John G. (1967): "Folding and fracturing of rocks." Mc Graw Hill Book Company 568. +#' #' Vollmer, Frederick W. (1990): "An application of eigenvalue methods to structural domain analysis." Geological Society of America Bulletin 102.6: 786-791. -#' +#' #' Vollmer, Frederick W. (2020): "Representing Progressive Fabric Paths on a Triangular Plot Using a Fabric Density Index and Crystal Axes Eigenvector Barycenters." Geological Society of America Abstracts. Vol. 52. -#' +#' #' Watterson, Juan. (1968): "Homogeneous deformation of the gneisses of Vesterland, south-west Greenland". No. 78. CA Reitzel. -#' +#' #' Woodcock, N. H. (1977): "Specification of fabric shapes using an eigenvalue method." Geological Society of America Bulletin 88.9: 1231-1236. -#' +#' #' @examples #' set.seed(1) #' mu <- Line(120, 50) @@ -932,10 +934,10 @@ or_shape_params <- function(x) { } else { kind <- "LS" } - + # Vollmer N <- nrow(x) - + P <- eig$values[1] - eig$values[2] # Point index (Vollmer, 1990) G <- 2 * (eig$values[2] - eig$values[3]) # Girdle index (Vollmer, 1990) R <- 3 * eig$values[3] # Random index (Vollmer, 1990) @@ -945,8 +947,8 @@ or_shape_params <- function(x) { us <- (15 * N / 2) * sum((eig$values[1] - 1 / 3)^2, (eig$values[2] - 1 / 3)^2, (eig$values[3] - 1 / 3)^2) # Uniformity statistic of Mardia D <- (us / (5 * N))^0.5 # D of Vollmer 2020 - - Vollmer <- c(P = P, G = G, R = R, B = B, C=C, I=I, D = D) + + Vollmer <- c(P = P, G = G, R = R, B = B, C = C, I = I, D = D) Lisle_intensity <- 7.5 * sum((eig$values - 1 / 3)^2) @@ -967,13 +969,13 @@ or_shape_params <- function(x) { Woodcock <- c(strength = e13, shape = K) Watterson_intensity <- Rxy + Ryz - 1 - # JPF + # JPF # hom.dens <- projection(hom.cpo, upper.proj(hom.cpo), stereonet) # # hom.kde <- if(bw != "NA"){kde2d(unlist(hom.dens[[3]][1]), unlist(hom.dens[[3]][2]), h = bw/100*2.4, lims = kde.lims)$z} else{kde2d(unlist(hom.dens[[3]][1]), unlist(hom.dens[[3]][2]), lims = kde.lims)$z} # hom.kde <- kde2d(unlist(hom.dens[[1]]), unlist(hom.dens[[2]]), h = bw / 100 * 2, lims = kde.lims)$z # hom.norm <- norm(hom.kde, type = "2") # dens.norm <- norm(kde, type = "2") - + list( stretch_ratios = stretch_ratios, strain_ratios = strain_ratios, @@ -1037,7 +1039,12 @@ or_shape_params <- function(x) { #' stereo_point(x_centered, col = "black") #' stereo_point(Line(c(0, 90, 180), c(0, 0, 90)), col = 2:4, lab = c("E3", "E2", "E1")) center <- function(x, max_vertical = FALSE) { - x_cart <- to_vec(x) + transform <- is.spherical(x) + if (transform) { + x_cart <- to_vec(x) + } else { + x_cart <- x + } x_or <- ortensor(x_cart) x_svd <- svd(x_or) # Singular Value Decomposition of a Matrix x_eigen <- list( @@ -1057,7 +1064,7 @@ center <- function(x, max_vertical = FALSE) { } else { x_cent <- x_trans } - if (is.spherical(x)) { + if (transform) { x_cent <- to_spherical(x_cent, class(x)) } x_cent diff --git a/R/mohR.R b/R/mohR.R index 9df6025..5afdb37 100644 --- a/R/mohR.R +++ b/R/mohR.R @@ -313,7 +313,7 @@ ggMohr <- function(s1, s2, s3, coulomb = c(70, 0.6), sliding = 0.81, units = "MP # theta.slope <- -atan(2*theta.f) ggplot() + - geom_circle(aes(x0 = circle13.m, y0 = 0, r = circle13.r), fill = fill, alpha = alpha) + + geom_circle(aes(x0 = circle13.m, y0 = 0, r = circle13.r), fill = fill, alpha = alpha) + geom_circle(aes(x0 = circle23.m, y0 = 0, r = circle23.r), fill = "white") + geom_circle(aes(x0 = circle12.m, y0 = 0, r = circle12.r), fill = "white") + { diff --git a/README.Rmd b/README.Rmd index 8b8d317..a984750 100644 --- a/README.Rmd +++ b/README.Rmd @@ -23,7 +23,7 @@ knitr::opts_chunk$set( `{structur}` is a free and open-source package for R that provides tools for structural geology. You can -- analyze and visualize orientation data of structural geology (including, stereographic projecitons, contouring, fabric plots, and statistics.) +- analyze and visualize orientation data of structural geology (including, stereographic projections, contouring, fabric plots, and statistics.) - analyze stress (including visualization of the magnitudes of stress in the Mohr circle and extracting the maximum horizontal stress of a 3D stress tensor) @@ -51,7 +51,7 @@ library(ggplot2) data(example_planes) planes <- Plane(example_planes$dipdir, example_planes$dip) -fabric <- or_shape_params(planes)$Vollmer['D'] +fabric <- or_shape_params(planes)$Vollmer["D"] ggstereo() + geom_contourf_stereo(gg(planes)) + diff --git a/README.md b/README.md index e29078d..821ab65 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ for structural geology. You can - analyze and visualize orientation data of structural geology - (including, stereographic projecitons, contouring, fabric plots, and + (including, stereographic projections, contouring, fabric plots, and statistics.) - analyze stress (including visualization of the magnitudes of stress in @@ -47,7 +47,7 @@ library(ggplot2) data(example_planes) planes <- Plane(example_planes$dipdir, example_planes$dip) -fabric <- or_shape_params(planes)$Vollmer['D'] +fabric <- or_shape_params(planes)$Vollmer["D"] ggstereo() + geom_contourf_stereo(gg(planes)) + diff --git a/man/figures/README-stereo-1.png b/man/figures/README-stereo-1.png index f818c712a389e53074f2862406c6ecc8862c584e..06e2a680b7bfddd36aa2817e30a163c8628af7a3 100644 GIT binary patch literal 12897 zcma)j2|QF^`1hT$lYPrDRI(<1FUdI@)( z#~)amdyyRgoC0(-u9`l2wK&2WJk=T~^58GtRi`XhvO+_?foNbQ#j7*mJe23GZ8ME7 zo}gJ?82S@hUl8zX>qtDxHz)}1KoD3s+Ul9pUXXpu5^#n6r+{m1C@5`MyUs!Y17J{m z2@QDKjDetD!UlwyA_3-NG{9r*1xV)2=`;8Zj|XhJ{R(AY_K_1&^-pQc;Ozc*scQ;=Ez5uyXwSd7GFA@)MS5l80cpOzJc6L-bci0IQ9E2xv0fjCzms{2h)P!;gGQ z*|OVlAH|Q%P`cbNT5s*9$?PBSe%{|@4ccDvZMHtz9%wdcQdeuJM+zU2m#35Je>-Mx z!agidUXiYfn2y+Lj!Be2foDZ7yXGeh^6(Pr!y;RwW}Jwyw7L4cqem7q{?X7;^ZZre zra#@SN0Tn0(v1yMbh_sjl;AW@??ge{;BzW3jg(&844Wzr-d{NC&EYJ^STQ+vE}4V~ zSEO7EF6efMZC?h64R#E=^5Cg&6_K??Q3w0AN^}(O?WVJOsOz(sZ4oCiA0pEyYM)gT z#&yR-S6?@3k14m%`%34qzr5zO8@fYs(84J$X91qOdNnrtdxp69WKZtk-76OUZM#}W zo4V_k+AixmW5KV@R;00xymB;Q(m|T-{=)d*tGoX1)t(r*8r4ouzhUF4Ww>0*kI!Vx zD6PXP=}2Q|RqNlRETg$4@h{91pMh#)E1sH9xT-WU3Le!r{jRT-Bz$_RV7mS#zK&~G zX1cL#awd}psv+y=4YCxbTInryPk6W_WZZMCu?f2HdvEHS|EA^sC1e8z-{`Zhp@J^Q zVBU!;R=*`G&}eO&gC7*FeVaNJQ)aO z=b~XP2F*Xm#hm8ppmQBB!1k~4`^Q!0XGIbJD{iUa>;IpD-DL;4u4AXbnz)%7hM^So zQ;P5d@-RG$ppG!0hjhYhPg1^_MaEH1$B2uP&I|krs@R(>+?Hn@UfQ@cGVil|iGiPD zboBOW?&3;w&o(Kve3sd!ahC0n?YtPKIX9hsc@}+YFVodbSy(pM@6x_rcx#maKI<#Z zcn|xd39U9qf71BoqlfalTNWoN?iZOl>AyxLa7*q!PxDh%dpGJ`-Y61pqwM~7e}a7H zbXs(z!o@|DlgJO7G3R4kc!~x{;RPyn#1m zVoZM{u(=5u;%8S3x%?4zTYr^`If20$J6b<2zf~f7n^5iQx#wZOB(I7e#3bkr82<)k%F2g;E-JLOm&?`y; zAg1=XXRxxQ;k(!9PlDn470QDD3JXu(nw^UOjUP@lX`B-6Y|pyllNM}-+8a5gNQWqX zvRZWh*g;xVc!nOcd|p)+bDuGGm+-h}!YQb`(r?p&nZzv=ckL2M!Dc|NLE=YhtIs*I zjz?Vxy=(9Ky7%WV%-OLyW*6$#M{nFq8!c9s;G1&^4fn7jAYi#wg(E5@-a|2`l4vaB zFzFx11%$I;l^_wc8p(oXp2T6?%^@eYd)Z3_+#llmOV`j|P)ne0$X%3+)vqK@WgXUK zTB}$Ij?Pz`-gYN2Ls}DV*sBTat>54(S3=TMZ(0qFx97qz1bjE(X=qyHDMllX5_~y_4t5^=W-&Cpro~OEvcZ8}*n}Q5@XHs@y#l*S|*7 zIY}-iYF#Rt(JI&+36H|`?Fh4@ek&9jelR-0la?cDvSpa7%nH3Ajk<=_h49I$BK2pB z&vT%~&6X8f6J9kDeG*FTi(ZFkbr;6+*fVI5O;&{G&t^Fe24+4`1uj6xCUvh;gk!0?t053 z$;LBXh=9zOo|21Z_4(Wwkx(V`q#b3di1Ug|C9&J6Cqslv@uG>u-7mGf1a(^}G2SJq zRx-yNti(IXnZfTwe*k`N*yshufZ|6`KaI#3@-`{lK0%Z3^??-!vDOt~I1uR_n=k$- zz>D6L^LQmY-GQ*3ZTV#xc`JQS3{^T!`LP^`q>|0&vsFT*eBT^UQN+-||LcV32*}x! z@8h>v?8q?qs7zRrku^W<;NI4IKj;)H-pV_!d1g1{F*aCwP)RmR7NX7Z`lo@GCxh_~ zm2;yYHiPWd7Hj!j6`C9L)j78od4Gc@tgBwmYyq)AW?q*;6z%KX1G70c^1@ZF1uQ9G zIp*}~l)`uQ8A^y(&x3+ap%aPRyDRncvEjkXbxXgktaj%zc0*U9Rn!`1>!RUue>19_;n^pW6EC zGnEqjlNwb>fVM#uvuusRnLBsL*>h5s);1pT=In~E=il>O{uQ=n^$QBDMe_=0dGJ}xiGkd%i(oV(l3+9DeY6uG zd2R&Mk(*l=^^R0f_NJ5vuTHzK=+g0ghOmrJt)S8ez2_JqNbp;Q#8@^?&!J2N_pc9 z#sL{<)xvJI#egFfmWMlXoL6k{eZHGa_`|PjopXe|HEC!lL8w#GiP*Yn)jwHHGVB%3 za4-&}XTiRSvrGRbReQDnJ{L78Qs;-C)ZI+g{W|@Tc~&s!Q4krqk9^&K^lrQ`HR%1U zO;c`R>L}avA|~RnCj#?fNm0css#_({w@D3jQ|({=hbvGpH{60n&;Rkkb{Ifk_3>8S z)3vo|`_oVw<^-|u61{2P_p$_iRcJF><`T9|?w<#=1jNn$04qYKCNS+r`e^@g7>;hj z0wwN!f5rbOYaD)D2hlBNFvUP|xT608JG)_y181ENxLjwl`CLqZ9I1S;Cg{yP zDkR$x9%9uNR$`(_QEBczQaAsasZhG%J{_mu=6E0SP*nl__jL+O*RL8tjBa$ei`1~7 zd2c*;U*UT$^2m;R#O(zCbTjy7QhnAl&U^J*gjawrRYa&ra{i4Bks^;P^XcFz_W}gn-{^Hqn4?(<;%8$##h+u!wv%-=h{pxhq z()IL&vqieuo333HgZjANTZqW&MuZ?H)eWEpVWgzMP|9R=_ldq%4J=ttdCtj(=z}`D!${ z#)#a%qhD*hW(4bqh>kMmsJ(f9Qk@sNX_H;WTRB5y?(;rT0Bme0J~4(0Nble0OJKx; z_55U9^$TH%?EGRgv6O$X^Y5^WDDL-KrSU3-I?#&hy$~zYobE(?jv?ePN*w>bb-1W1 z8j`6AUz|1Ijv9slVv7uG-<+oqsqOT^WbTRetFvEg8S%?uU7MRG%vg8i@>(Q8b<@~B z_#XEV4=IozJ@BD+>$fK?a7y3?TbN||LVrTG0lX7mR^0!6eaTw25UAdny7R_Jkd$1v zA}J92{GIZ(IcsvQ%K0jo*FpdTwj1|cGAl1<8Tp3vBZ@T3(2jj7%R2W?Px*E*0nLIH zD3YJQ3@hrFC|j@|$!jGzNUD1|59MDqi2PH1jlV+%98nd64kFR>KVf~P<4DjAxpO_l zkb8GE2LFa${_$FvjV~wb&_4epJi8-bU@=#F=9@m~Q)cZB*iYk7*a(a)o)?r55}$Uy zgr}8QiwHWj@K_X#xVE zMeZ{8y#tHhDPKPACV1QfHaEKFCumwi$tBFG&UUgBtMBjUo#BbAyr**d;^TzfYa%ZN z2W(A4hP0k{r3esqnX$S2BGd^?P*t%csFp%0D9s{xV$5$LI#a}5dD{+}ja?$+bl}cO zNV*$FA?t?3mQQ)`cFQr<<@WErp-;3um1v`uXsgtx8ijBC4(t0c=>?rJhn!imfUP?2OorTa2>s^XW#~F9Cg4RnF9lkwbVpEvg2Zii zts|I8H?qwkvzNjC zpf?Lybsjsg49Go}*i$QHLAcJ!yx}$BA=N=pg_F$k{~^7ErHe&vW^=y)hTAH?g`eiu+wo8ObD~G_`Do=7`>PBI(h_|hN$IMb83p8ZLOV$ zrTv1dumOtE*f^447G|1R3BZF}D`O4U#HRj6Z|fXPZ~hiJa&-gAJ^?}`#_XUhi4Iouw>dpr`! zOFs;`(4b4g%Cnj^QWqrQE#Y&6EYIX(|9Xm)n>o=S9D5{r69@Z!wCDjK7hQ> zlhaV+^H+JkyKmF+wJJFZJCE3IU_@RSk312T<*hWGVg^H662b}0*#W%4Q08{HbLC~yvGui?`ZIWPE`#=-I!)yeJC)enr#cysx?5(kv&mrZ2qtxcu zoTu04UB?^)knCX6#x%^?KseQ9MZo$A9P8_q;JgiB&ynrYTvk$9isTm`NrX#&Bggin zoV6{z`}Qjlj$+2ePRUogRsv=P$mOYW1-|H>FL5(jx61Fl8_|H)7(1u#t!VIN2SFXg zdl#EdYRQZWKW8X|Q@&)dsdB%76``pzEmk`2?wjm^5*M!;5gY^L0kcBnvK`cM7j@N% zYR{L~OL+!L*Dlw@Zo!W!x5OxIo~R{ZLe3Q7y<94rX@+-wLWpj3)}@yYYoJ2k-!_vc z4nH5={4YvWOWzEF;4rudDO@CnIFn))AGDZ34uhw%O~RZnYH}is}d8mh+H@ zDkP<`=2rFI4yePBV?B43W6!2}MR$rZkUC5*U|-dv_Zqv8xKd64238g@NOG{FD|7Tm z!pk!}gjFwIc!&IZ!!ZoxCNU}a-II@K_k#g{P!(4N&$#%)DZI3d6+zR7{;V{1NA6WA zn4)Use-cQ0D-G0vpBXp^m1Zy!9UROZpaoMHffyT;gTOpRdcyT`CW0&{fM^ge8c)&Q zHVp)*197rj&SC3m%W#>0|9&l!hGZ|#No{@EDc8yw-Ls8dHz7AKn*Y>to+KV>{Fmxe znH~o_3sHvBpqPMv(5yul7Jsf-XShgSg;CC=()Ydeem28D`IV|hh5>+Bxj;wN(uGvp zXH}Ifm%jJYW>nnysV#`YGlq&4U7-%7Tj<5$+4NMxMdQ-WR2{J2sjU*@bm8>^40jpm z`+63dBi=(d{oJmZyyH9XMta{R55tO*2}|FIlY_h3Zm1J8-%-#W6;Xj1+Ya)axaoEV z{OlyL76uv3!yr96!u2-{1lc4nj0!!}$?%_79N0If(>745pUDb1$WPqlPsHNc?$!XL}Lsqns6-ZG7zhB84rgKE|= zDjJq8#`Iw&7$xzwXpg)WpF%l+IA^yV_52}IK;jggbHH%L;%8Q31~d5DAk+PmRFV`J zj_0&uQ&x~Chj&0w+dvE9h@aBpv z)X$hfi^mio%TE@OXV8SDcH*QH=#n|Vv}qYbastm`Si#qgR!`Y+V(Kue7Je}I&2tly z*@2BVb|PE$>&zdPG$d*`2+IP9k4gY5cUQ>>wZa0z=+753Z@W2wx#)h01{lc_4hAy8 z->Bn$pM;Y}tU|2PGXVUL1#Z$*<4fM6UV1txze)>C=yrH-2`3H0!4VIO1pr?IU>%ju zcPk#NlI-Di_U;8l#r~Ydp3)7wg1qXh7%Lb8jagHQa@qVjSY-in;|KrwG3$xIx(!{(uze=V$L!#n)IQ^Bj@H-4!nUkZQQh>v>`v_&2b*5iLy8k`u5N2T%_+rwI-N zIGx~LNzwyu9zw4oFSj$lVooRJThT_6UmFt)>T!{#>sy+&(}+w_2^;m$1#s(W%&=-Q z1yca;G@?!MZPPhsa%Q3oS_WUg@8~!=U>r*N3p;Ua&8gGHRWjkm{=S)+N|I3W_|FHYl z)q2{xGGBq;ZgqTM#N<__1jZ_+9W07FNMpT*E~Hb3xjTQ)cTIH8#|=(-oxL(yXM$`~elf$kJ$6mQ+e^LIh*b5rbcOz9Cgv82DB{SjJvb=NXMm_NbG$ScaZI7(i03?~d^9$LK_%S`! z;QRR6zwLJRE%V!SyIOa0#76PAati!u3FetKgiVP_)sXHvvdlHT`uQS&^d+gEsc-!I zzGORcqXNJ;q5xB_H149qJM{yFhdaN_Dci|(6m|7Ug6hzVjkhWKET9Ba^~4te-Fv$xPq0-D%%G#rAKjd~GkxtQ z7*m}&S?awPZ%?oZ14$MZ0wi3PXgr*BrSd+cqij$FtJY(QSrCO$UA=uH$;BqC!hU`N zoQ!1x%zjvyT$19)q)=y4iT^XHjRn@=AdJM>v9N@n%qXPaIFR)6PJ#o-xTOg>@<+uI>$Y@pU!WVE`;VHir^ z0*(*L{&B%jG?4zAr3A@0%mH3vsN_I3-%}SU=QH`?U}(jjazzXml7ZNCqQ-PyQhk#B z)>~F*s;zm>M}8esN`o6OUyq_g*8py`xX)P9ZMo5xwyuUshp)|G88<24K{*RLd{zN+ z&5KX?^eZ+_V!9 z?j^3STf1QBw;y12QM16h^9<^CA46ZW0*Ta~xFpF|^iA^f5`a!25#j6Y$!8)O zG<0F-<=2{%P*9KR$i{vzn7~F4RKM}=I>*c-)0ev+?)YmOf z(8G`vXsyu(?_5#gMgpA`qb*@*|COgBr+#$wR4Ti0rZM1#vH|GFhAdm36wj5sOU;FD zzFZ^%7D!9GDT69!^%% zSfx1*54~%g+HjXCBUS~-Xhn-#Ptdwpi%%Vq%#dg_X|~1nJ~bSg95LKA)Tnsb$wFPg zOh;pbMf=Q~KZF?PY_@yu+Fsw@h`ApXyUjDoa_l9Bbu^3D#lF6lxtb!U-mR>@=#dd< zJh$@?R?ujA>Cn#AQS}nlf>qoibVNJ3;O)T~IrYi;ZKL->dVTMzM3zM8@ZbLtm%|8>&RG1a|=mZRcjQoajK-6PJqvdK2)1`n74 zFS?8s?JR`*7lq<#Yb`G4ukx!|+5$x;I0Pxd1G%#gZaB5&=?)Yk-B(ghm9tZNL(`PO zHkJktMguYQ#A1w~_A*0V962_>T26Vc3&xWrV`W8oT5AEPs2VCITd0*|*|9n}%>BtQ zhtp|?OV_(Jn5nJ|z{~ekd1xez-I}|M=wla49 zR+vU)_hLyAj1c?kyPgV-gaA0{@vt^STBzeYpc>mk_J&(*Tkf=}MpP9VgN9RR5M z=XTlAFh!w>EQ%sai^`3%44nEEbToze=io0pffRUl#j| zj>8qBbNX)ixQCi*5#9gV)pTTp=y}zq8=kqd?tNcYn3LZqG2|=lqm>p!8W7?C z2K6vnt(*0&QX;F~#!te(O#c;-R13odRKbx(&thTd4FyBFV`LgXPO}|}J_3*CbEQ)i4rq4MHMixq}*3~8rOH-zaU9abZq zw_r0SmO*$L|{-9Ok-H80}UXU~?XF)KVTzvHq>+|#I)F!G+e z?kge<_N}Ek{uPkD+Ieimc|CqSniaffQlW{H!F(`2huJC7g=HGAc|O3k=$xt5qjAi5 zRqjKdHGc8m6^2Z!0XbLEFpXTpD}~}LdicF_4-ZlK_Y`RHlGiocq)65XDO8)}^nY04 zbnp6By|I(S2pS$+sftS@Eq{#Z5SAa;WXu&$=qej{BO@+vLK|TJV ztXHQqbHi*Z>5Du|DuxTWrvDmee{nzC;{hrxc&U+8vb;;(&l}=YA)S_)(4?7IS1~G-l{lB8h`ULEF9=qlv*pVDgN&6NB``n87wXZ^&pSZQqS_C{v@I z6w>Mwn>aT@-)H=6ch0x!L<~)uQ<@N9Fl6J*={%dFMKdig!(ZT)oXL!-{+`zqsTq{B zB#$a#93eZsvT_(>GE%PslDbz<_WLJd-uFvYU-kF3d`5Z>&Zwfy91a^eMAy%rLY1)2 z(eFZSb7U*ULoga^Ry3N$&J+VnTl~8JDf!kpxlFo4eF^^TWX~}v)j(yDK;^+t`7M%i z##vtk4)EEASpbh?q~F;JMz5Z=+wIWDx)hMaH#f#FSKL^(1r5enS#FR(_0GkJHiQVcV7$7eU4ghuX= zuW8B=4o?M&OBY1`j3Q}6mP&sC{1@{Cf>nf14n}r`v5qIJj6F7>L3Hi+=2Yu2Wb7lX z(;t4>a{iGTsb1nGKzZ3I4y$Hb++pk{U>soxi+~{l8=a*3F2hC$Cay(JBH@!7#iXh$w@fx+2O z>Ce~?U;n>x_$m3ok3ld4E;OPD&x`)UJ#rT6c<1|mQsLL(ZGP$Bhfv55)y>2q$obMJ z6Y*nMW5e+>TN-QAvfM(vH(W!buA)h?v_p>FWTu*MbDjx7hB-+cR;ZU@EjcKFhs4F7 z^Y0rIY6P`!1IJg7^$QGngGGg+YK)6=(>5>8Iq*lq?cK1mp1R(xPub}2fYZ_=v?b5P z+5yzj{Y`j($kLuGyH4*TXn2;>E_|rVP?$Opa#s^K9}ru%f->PPp5+a!jsZydW*&D%RInu z&w!2NcRJK3l;M%||9A*vWl~vTX05%t2%eV!#h}t{p2E}xH*9``Ga%N2N08J-6C1!z zr?zzop>>)=uBo+#m)@=;hff@E?=F?sfKu3Drrk`_4dc_&UV-|_&I8tmi#qyROWYyJ zh^xyC%;eZ>a<>CG5qu?H&$R{nc41sjkF7esQALWkYc7`Xbq!0Uil_lh>XF6ORg-`^9sPXZ4C z{6k%)>Ro+$gmsRp4!B{spU;o$k&B{($dko?Fs&@r{7^BNJDeBZjf4DmHCv5}?YiHE z(TAYJMt_D3p;uJ1pMzQTsw0Z~r+5PVjK7shokHCr?3QyR1EvipGHHrlP8Bwu{ynRw zsMIE9%7CzL*8=*V4Bbo*xkw?C6XZ{I%z8?K{w&zPvj!^uZDM#CtQPAsO!t-MkI6sW zG`#gRlH2QrIU|!C&HyW;eHI~>F6y;kbszp65gvGghAxY#*!k#2m^ge<&TRmmf4Tkk z?mwiDElsRYJ?k!R!?*eH;WIc-K)CS$TO|{cVZD_;k(JJi*}kinkY>@iVt7S0;n9!k zJr{j@|F6cKS~Ox_AE?0LF6wzs1*04ZQrk7!5&)Q$MU975<=e(u;M3^N=pMB^xMGuP z$TMfY`}UuwX0*tK2bGr`e=D5cSUhV`G)vwvuuka(9F(wCtaCih=Z(*UgAW3n!tA~$ zR_!Z?&Z$uI5X7oEI+^<@Epaeio)_uUdJcP9==)-;#Ub4R8~r5 zcfQYky<#6n&YTOXtk9--k$u;zhJ<@Ex)V~eZFc9uPmWg7qc^p2rZ6j>qTH3)&+Spp zTCD8KCMFw)b9F+XZxt5a@slpZ%J1sdA44povJ~U4sup8gUIU#aO ze9bfj-tFOZyFB}vARTp|x+)=0x@usdhFSFur7uty%-$AC7xzv#e_BA(e{V*z_R31L zu^uJ`^ZSab+rs7QOn&*&g6KLo%O?o@Z5+a$i*2RN& zUcreHmjuNxVPoYvDKWyoQ%^50D4(W23#HWrGS5yPRW_F@@nXR`=6;%r4DtT$dDm*l z|Cwa@PzETc2zh6-Cz}^MHf|h%|HC1v>Y4ZEKr8a$OV^w{77FVTsmJod{M}At&nE`V zJdX@9_x$4-r@@TExsBKWB9@h+oi2&>bI|eJ<3MgPLIv1!Mr#_qL}%-j_v>X=nuuS~ z+PLyvp>us?j7ZxHMJa+~8S%ZyhprtD#}`;75Rx%=A2jyQ zhExKN{^)<$>pg{FFcmAQy~>TZ3d#A?)~n=T?ZhuRE0?N1y}b2sxTE*aSx{3YPOwdH zIVmaM9#R^S-1IFDxOOHfubS7BF1Y_<60-<`S7+@5G2xjN!6xd({eb$G3obs7yUtiP z+n7C_1D}^Nfz5A*g!>mhR+*Yt+J+hBv+61^7x|ih)a9vFvZ&;D`R?a^{*v>Zq(1uU zZ$#!r$XPDKY|+{oRk^^;8HAL6GayEj;c4fl^s5K1b*GcR&Y(pE8dhwBodPwx)!w(hR;wQyS*hK-Lq-Ou z$OMub8vh)q_Ax)ya}RcfR%5c;0c*XW_!CafQcv_a}?yo+UQ5zhT|!O##ltq5C#RYM0B`BVVLmVEi+!49W z7O^tzd>OIrv3=#yVnkCDBJ+Ckea;t4luXxg-?g1sxq7UKlsmUOlN!C8eYyrFJ5HtU zh+rGF9z~h+hEuY9_?byJnG=tys42nLr!N{<9EJSC*r*PbRz_^^fqFVjQ7PJByUR&Q zcZk-Dq4w%Zmjpv(A|!V-Q?3+rvR&>S*R|j&oRp=;kGCxmLdewdrtEQV)c!S*Hfhy$ zi+@U!WidMI`Qvb}{2l1MUj zZ}D$ZAmkvMeZqbcf1b1v4cehCdQ{!J_s$-~w1`yyiHMzSZVnsiBUCBHOtiC&N$>D24lnlC}pM0`^ zvZH(=;3LxDtcdDZSSio#y?EI(E*#UmXtY*oYQ$1yBBW%@p$q!RH3eO(DNTkJb|&7v z*CQd_<->JJ*|0N)R=>uLzYWt*fB0t~NrD1y7OkMCw0Vl>^=d9sc*|0; z=}teb)iSHN5)U8ys2yChSZd=%++P_tTy690`7=Q1&6mcm5BslokE@(F_(J_POE|$q zxuhXzy<#gd{*fjh)ZCDB|Krgk|CRkWQ3m9Tl7pV>)1J_J|I@A0Lr~#FT;~1gh7*!k zM?Sfc?x+QRY}n53@;?wUvaMIz7Bn{`*)rCgzb`GdrH7o2n<@GN>a|azAbhd;8sP)X%e4lRr6|H^N+?Yl0_37Ns z!`JK1gW=164>h(9t4^5uZ4jS9g?0lkU9-;STU>Bo>(|NVj)`s4XmWTgdvxzWz-HxJ z_v}trL+3*<>&T?hT8)ij&F{aS0z>y{zdj&HqW@Yd=gs}!(s=(hLPsJ*Ss|^rbzS<$ Qe-+WuG}I_lzy0{X0b|6RvH$=8 literal 12978 zcmZ{L2Uru&w`d9w2-2jep(#ig>4E{m(xrzc3IbLLAfgz0OQ-@;6p&695C|yn6Dd-} zf`CX1(xgl8(tCUHcklmx@B7}n$?nPQ%$c*hXHK1&L}HC~nUE+X1Oj2wyQXCZfxtir z1p1pEPOULAlcJ+uFk%fYw5gX62pR%$fS}P3G!@1}92_7HR9Fe8MTv7=p|(*Q&=7~e z?SDViLOj~R0qsDAmFW0*bUYOfS@M}U{7Z$}+2LPoDux3({x7k=B?pK2c!&5QhssKa zN-8`F<+F}QL*oCX74HCvr=n6RQBmS6so(g@q4=T7_#rAhDd)OT=>Vy8Kvz0ARK}zK z5}=|CQIDa@p_9s!%F2_W%9E4&`ud@G^iaIRP$l*Fx9h*{LnlK+CnqD9v{!FyZY4Nx zf6m%p9670UI2nrn+wX59wKo->XumjDOl1bc<7+pcLLkfye@|$$8~Hv2avq|mrEcM! zu{;uePq5h*vAU(88!)Sk!+mi-eNC+VmF75oqo+?Omw_vBuwILK`}ib?VPiOR$Hmdh zQRkDClB289>aqvQrVrU4586G0~ zdTt6z`9}`n4&JAKvMDh>fcc25b*a({JugaJUBjR`^3yeXWLLarhpptOKekE#d<9{MTk3 zzmE9dgIN2x!;^350#lRaQRwPP_CGWzo5pP93&Gc_lwIiY;( zjKg^W0cx0(hfs?`O(kySLtjv#0HbKh&v~?jAV9O^29ZBBev0r`PW!=W1 zqTe#E>a4*l?P#Lb)joX_%!1BkAM8noD_&VcQnGvNr#b(jNqod)j>7zpGN3Ru%sra42v+47Q1Ie9dx!-(g|z zsrrJU9VnufKh~(0eLO=QM}@U1;K=09yiuq9Z^3JIz4?1b=dQ+^@L%GOa~#@G-aqzo z9Em$2Z>h8W;Ae^BzT(hHTPXLY!vXF6shWj_`=UCkH4XT5GT~8Q#6+k}-6G_e8g+DB zfVj#nLbV(6N%x0&sb%S0zWVBz( zPJYJ5fCbw<5--R^Rf5bVV2#GJ?s~%87~PY)1Jsdw>`9GQkpJhWuOfkAvERZZshR|q z`mXY8$bM54R#_+y(JDDr&P|ukQPJD5}hi zY_0C9NU@l?%CfxdcmhA#I;cH9kygDJa|z5n84m)3#%%fLY7}d>PoKm}Gsmd|dn-P6 ze6I)FD2S16pLIl;@Uyd%Cl^ixlIF&aAZeDbO;%Q3H~%xXQKxdPA{9C(iMiRnL2H+{ z*THO&V%SNHiW~nl>*EveoWMUf-s@)C8t{3J*_lZ{ZFlEGf)A^>?$2+zC(r%|iC6TG z)!m05U@F}NWXQ$3D=DyV3FB+$Bh|u}f{(Ihr`3y$X4HSSC}2{VbiAXFcR8Kf8-%R1IDFq)(zL6!z5Nu%^d#Ysg6m{xN#u%Vzi^@9%;5o3%PGu29P9D$ZF$E*34>EK zi^DP*Iv`u`eQ9hg$MYG(d$cW&Z)en|w*ZnzRj;kv^7pL8wMP-Jm23f3>9Mk*lDLX$CEkya`!Bm5)jYFkD7O zbrAE(B%I@!@aw_T8%)E+fFhC-<`3&Cx3IR$j!{;f}RWwpFvRA3aU`UerVt<|a%_0!dm4KA-1#rG@B%*9^x*&Z6( zaupfF6!tga?G4>boJ{TCgQqJD%7ZN|bY6Rm8PY5)2hVIfY3U26X;{f$>`i1nl3#mg zpp1RW&2LHUm@wDXk9GUvjX7Vz!MdXtYx^UnvIn6+9|J{`OnB#Q&ZFYdQI1oQ?<=4R zm22CV=O!?gWVcZQ*1kxTckqBGfw{+)czN3V8jIFQA^31=p&QyWj9@AStXB8v#7#|p zNTEgwY{_C@fc`ZGtr4;M+lhSvD6AFZr!*5oj?h6f@9OvkbNe2b#ZYvbx9qMIr)Z&8 z|K&wlbI*6uB);s0^Lm~?{9lp;YG_vL}C1&JH z!(QnY>^+A#^ECPqN{kX{(|i;ja7BD^6sNda0_x;60fUCma>|s`Hr_m};sxKnO7E+@ zL~SzRwWc|4aLjFV`Jf5mUh0bwq(&o$)#=Q}SvF{xgWWYJ+B;kG@CDnSku0ItVH}~c zVsClcModJoGaEW8Xek`8{{+dJ5KOu7fD^v$WMlfL5YPrE_&s1kyj>M)F#={5Vpx&` zSkG0j0G?0Gf30yFlR(f9!0NsEg{a`Cw@uHzaGs+8l>LqACg`(O42G{Y!owPjl%^9u z|CWPeC#CDr;Q4YXdR&o-6gHiLq~6M%5Q~Wc6t9DJBn3Eab{>S&f1lW73Epu%OhejU z_04`z=hx-PgJm&Ru`WW*b@7);O|%9O_@PDM6(vfCy`o{sd1!))DoL@Q<6x6R4ufg8<#x}@BXi+0>tt%r0N{v%C>qPJNbc=qbl%}Ub}BQ`Dpy`Hl1&a+tuHL z>>gK0qhm~QLRHx(JiOTVCfjll&a~-Cnlx=uS#PCfH>q}DP4dQc?V@URn!eZQSnTp1 zbjN0kZxMErm;6T#MjjnG#aeUKESr9ABl+ZDFZ0Zt(lg)s^O;?fv{Amr}g$T4X2(kzkG<97QMG?0gj5}c|F-u#W~=O zpC3Psi92(2C}$Qh>RE6o-R-2t;~=2m(9bl+-PK^&Olz@o?(uUmUa~1p+x4){d?`8k zOv)AkZYhTOV8!~VCjFQ5rF;@aDH zsBZiehJ;!GXJ3hdFhOsw9pI|U z*D-lZ$50^ir3WF+zX?9F-pakggL#nKA&#eYTT9g@-zg$XN?FbBH*nJw9h{;mGIBsd z#2FqWZhI1wGq%J8AQ}Rg(&H z`~}E9cSYgi53j|pGqQO5p2RIFAI{%|<3ddovHkR;!WOhOg-kei!8j4(dAo+HS07Zz z%yri-OfhUd@-h~%A1l9WQGAbIxn18t0_~Urn8!XxXdjEbrD!^rtl{!~IiSV%=YL}- zQezyyJDT-x#-WaT*Dn{tdZwPzflJClmZHlcsKMR5Y#|XN(UJq~mT-FpH!;dwdq4}@ z{S8|k>FBj;yM~xrO1R8ok$$)NCGTA95m{pyk;diuDoZ{)qCcn(y&?(Z73Sat$Y>dk z1Lb5wnKG3&)Dj-Eaqg5w<=(UR4Qd(>7S7a2zu$Zf4LDwM=@B+CSJnl8T6ATV_vA-r ze}pA=y5}3I_yO&JBrl~7;MWD?G#(5xQ?#&r2NK^{}fAbrXf-hgjHf@iX^) zOHkN|fvR$4N7w7=(uzZRuk-p5YZgsxaqK|c>lgOxMm7EmK(ntr54nBzyhV|Q1j3ic za;^+0W~yVhAnUpU6bDX$2^43>SBK_}e%rSlo5)+06he*^=Woa8DRE?CG?(}WWi6Go}8Cw znA7FZn#r z6H$jjsQ4fYHqCyO0qf~NSv!?2fbITEm;-A9k4s<@ZHN@ds*_&MhAVLY_a|}!0 zA~WFMcf5Y~Tx9N%R(5N-gy=CR8q*f_l9bo)LG@>Q`&s05$n+K|-x+|G?p+Cr1;~;0 zYm5eLA{VSVoYjuc7ao-kz_(nY>#-q*dusIN8urY>`QIgktDc(6hk|SxBq`&)^N}V~ z_eBs%bb5!w*-MyzAL-m@HK#Ph0sLL%PhBcE^BiNeA{&R z8J)j5oPCJ-mqGN{l?%dY7eN*jrJ0vsC=#Ke4feGNaKc=VB=*7x;k$T-Ef4hw44G#0 z3C)7i4~vn#)G7JEiUFZfo7=Se!Vmk(prTo^Z7%u(NqlfsW-!JoyA{H;nM#ju9`u%c zYt7r=JSNXj?v%|~SFwte$IbJo9rb&k?Llt7sM_}r@J|4SN;J+uQA=|7w`7#VQ(>Dr zuo@>?T(`z-kmH?*l<&>|(w@YQP zEs`)^iMANoBWsBIQ@N*t0%sz)vS>Xpch$Q&fXurPy^#B$N+@kDe&y zaaTg+D_~4+lWSJ%f#6Va`>aktniz1sQy$bBFW7ngaw%-ANA0^CmFmVkZ4DdY(Fjnb z|9ac@+?Rgl9;>i4{%;|i6WXg56T9jAi0UW}cAz^}W}F6|G;_K2$W+QCfX2SMGZJ|B zm4M=KE;LR`LD5@le(b=AC3)|H=40;CY5R5& z8d~s+#)xEXSALk|5k}RLyvFwR=-}z{1*wz!gQo0b7xUuA$99tg=(Tu(Ts({-D(`?)CZ8NTB&+O^FK+hhjm7Jog3z5VDPy!%F<5Nsrm10oDC7 ze|`PLz?LR;^%e*Fff65qS-F@I$IL}bOW5~9kIcFQa`B=PEkNFbW|N7=sz!zHQdV11 z${kR=Sa!EP&gos{8Xx$x>|T9+ik5XqU9*w|>5<~RpEuv{)D(;VG~7B_+su_JLZDr9 zzO&v=Ns$NF%kH^H-{~4|wK{R{g-}jO`ZI%9?X<)6;Nnuse+KsNo*th$+GB+9!juy3 zx=fA$(+`^c17sr0QPm7!brcGpf02n?x8!}}Q>5Ie;hFZ#ekOVIGl_JMg)cHt>u+op zmZ+C}q-*c{+70QLy9iB6=6%1_xKK99M?k^Rxp$o=Re?B%fhQeYcGm*4YDX>~ytyRc z1jIcUc+z*7d|_argZIEm+6jnm;{o1Rty-7KFF-g5T$$m-v>|VZ{*ICxBT1%xTivE+ z|H-zZtk)Ti)29Wolc>;PQqZ3^Ncl^y`gsWUfEkK|%2cjF5|>!5cVRGKS(yeDq9F`! zBYspJIg||QXnQqptgncj|V?(XJoOU4O zJn@EZ@NH)wQM@&l%HWr44s{j`ZuCG%+2ewL@-bwmvnzB(w#&PXH$q1H13vC*6R%t5 zN%CrPE|jWoZ8cuh;NLWWp!gsQ*W|H215c#*98N!~FRqWYUzLM^v08A;5>G>-b8Bpu zTj`%7zKfE37c6?kZlSw?xH|(fhMZy||CXM#>ZUpys89spfDIEwVVQ}*#SU>E%tRIF7KOOD?{Fi1mMda)|n^u`T;|1LN>)Q&wKR##~TPNQjFn7AeJ2!Y`zP>q@-_WFV4a@BT_P zRjOLNori(q`{=fXSONLQ$?$!OqKG9&`|66x@T;Mx4 z}*=vdI>} zMTtG6!SSZXc&&GUBRlc>lgr$Y8XN?7KLZm@)$pxMD5b(ai}glt;mX@147u4+j~Wxz zzn{@vvm*(?KgtNr&Z_A`a29Q%rj0CsBN7r5*(h=V1*xpn`S!C2PB{k$1YtmOXBW)q z>h0T$-WvS)dL9W%qGrjqv-x)(sVidpk=zj65Wu_a%l(&a_XF0v`9B4V41j~wU?4M{ zf#CcvH$8Y-iB8kp7U+=@0dOqy^+r^`gM&8(p-04dcw^HWw%iXj?mal6+?S1|-Cx-Um?HO_zb@*Z@o$~35 zp}0%sk24kERHbJ05c^n_gu(3rxgWG0m%SA84dy}?7NdR6V(Xz0CZo5pK0dE#EspSq zx!$e;GHQo-JkH4b_E<;@QiJ70m1xg}bhMcH|1g6A?X@=4_v-m6 zdyvjzEZG0Dm}PE|W+Zj5T~|x(l(nNXyrvTHnh97jmD%-|7FwAlsX^pjW&$>3py13% z1*}JQs5{wb%%3XT4}v!M9DgiRxQ6_F>>UC&OrhYbBPd1SF6tYc68PuJXp>m6swWNM z)k`iyzO%bf@8NSdjCg*oA*8(YEH2;LyGfCsv)dHQkCnZtWo@rA+UypHUaV>L6^U*; z)#I?Zj0bDm9YO+y|Kyzon^s<^04;cBj{q6QPv=a^nxCe+WIpQel;|mPxx1|UxKOD24Or>9WpbLpE` zN?tF!XI1f?n>1uoq%^p&kt^r|>OX51))2SlpzSUj#SxBWT{lPl?w=<7P#_N@ICa-T z9U6fHTPsHFr_PJH@9TEHUxXwbfL}f0{Kg$--aCUETk1?_g`bZ#c;}NdXGe5}prV>0 zrBc_Ec~_J(igefHm?xPqWBGj0%x~Qd`1~qBA8QA0loBvr#`Qu(KF3*X(OybIs`*^ha8Y(I40XJF`$^fev4R>OQw>Rg-zar_qcg3D@IxbEAN(>QK;#nQAoy-q|CYuibhbaL9WjO# zAu|hUH14{LL+AeKHVZ?&bJfF!83ql6GxXNJZKS6KM>DU};9U$^%^)>^I2>(+`=?X> zbsl|1Kt2RNyy#!x_I;uIi`*C71)vy_3Yq)RFYtzfBfWbhfe@Ee7G8{QKKLfYZz-kh z;3xTH`qt$v=3F621?g-L{I8$tKl7fa$>r&aOh;P5P}fL@cT~O#rrFfbI9D4|5g_`X z;9#WDP3L*gCN0E6B+kTwXv`_jme-uwSE|8VaO_FdWdC@ZQ_(Ya_4~E^t-M5&9~#bhJDxpN7XaC_et!67N5=d@g87F zwBNT)bxzXyoqR55%2(rxjDZ3RXzPI;FZoE2;cVlZ#psGRyj3*#(N?ivqG^ufLI@RmOP$kM7ZV*#q)6O*7!hoPm$rlG#_t+9JqR@HME=Wy-L`E+&A&N@~~&dXRlVH#e$KW8}>=gj}}hUc$RV zTV$IzSSm|K+P zORhptr;9SBF7-UJSk2|TirjP(S)i!_u@Ka!jYi^Nr~}klFzU?|jNlt_4Pr6oC#E3> zt*F|#4nj5LGJD39LlNNqXfnc?OPu*-A_m8SmFSAGB=69!jNY&QC-LMgJ49hAkTzR5 zXWLZ8G@{Lz!x^!_nZ`t3X;_6)L!0+L9x4rKuiyNcs)=b4!0)&~6zD79jT&A{x^bo@ zEV%oeX$ahBbzc^0wDY91@67vCieq*!D=sh@R^6+$rzP9qXt1rk)o&;~0ThB1)vW!O zd#^R9qn91?Ol>dk&>;sNY+R`+Y;f=M2&tfDLYVIN4WOW zPg>NKNx-U5Lm|85N}%I>$e&*)S3?(E63hq=tUT)qb=>4KTx_uv-yhy&@o}ahsP$dW zBDg};s`9~H3tV0w7wdcl$CtoO_N=%ki1lQf$|t-s4MT^b3D!-#8N7u2LTkRbhx!Pr*ommh%RXggL{c}N?^E`_fI5f8h~xRlRq?*wtKe-t4#TI zZ_8b%MghiDU6<^?^2c?*E1!x4&Gq&Ocvrc!yBEy!*qH$*7Lj(V;dPHvD2jm8Ve(q0 ze2fI(jz9efd;?|o#>EC{QBplwcrf6%p)b|q*?3-3(Pgt2xk}&?kgs4t)o`>R%D+P> z?xDiWQ<@T3h;z_Of&Axf7dIVTv%k1prB1orysp&4aVn8dEbaWJFo_bn^I&!2It0R| z{&z3Hxm(1()I0*0v{Vl)k>7Ijo-$%2EEtAQ3~grvxcBJ7P8=*WoU`pjW~zI?^gbFSyg}Dj1%wrwf^kt5EH~h)g)Y8 zegWRASOMidUDEB)BP^eRdzkud6M&L8BJ&wy4lRsqe2CUzsKY%IWUt#T-*wcCTj) z1m7(23S{P&0v588)!rAD#9=oM=MOoK81Tg`eCx05CTulM27J|yw;!pvax+xT(wZ|Z z!cmsO2Igj45zw($ePWt8jHYC$vg`9VBySV2qwyJ|=1qImFN-KQ{r7d&c*vB{+J&z1 zta+{3jd?iL^v*0s57-_$lnFD8^a@ihW~cg77Z5c0d)gkWG}O}EIK3d`td63UFWDUd zk|f&SQB9w}ky*C_PG3+hu9MOLI;-yOevq?$1`LblWGk}Atfq33>u|CEEp4;W%!-FG zFeBEQ3==i%Efe>IuchYCj9F1ZjWmi_*MWEQzH25cB?k>}9MX5avp04P-vHY2wYnc= z6;5wGpEJ^eM+*mXnHy3q>lv1vpXu)=`d;z2np7bvJQ&^d4}avL73Y#jB7`|ZOzu8N z4YNDn^?)xljz-N0!n+oW#*El^g;fl_k@*m#>WgNdMeS_u-%2WD@G+U`dO=@PUie@0 znA1(0>h7yyPvck1@lAZ2^7sOiFp5Rq3H1>%*L7L>dX{uPl4*QiuILdud;eSU!5_By zp5$#qf&bvF{l$^dz>!UqOi#G*O zd$+T$F^%Rtuj83s5cty{Vh-mW!8{G+@4)d%!U^Cyfdz^&z44fQ|F_mZ#u=tzc1Bt z)cFIpW!A=>A$1<(bWp#s;m_*vgKgowx+v{K0hCaWt+1pWekjXh)a9|F1uQBchtz_p zHsv>oHfei26?weBotySsUml9vk3`kW39TT~)g+GMWGu*LA;xg=sujDeDc z_*N1s66Ke(rqRvq4MDAy6JzSKS&AS`Zktu7DXV?9+hQ+x2ot^*!l=QHW>_W#N(;(# zatmAVji8iJ-+ldGZmZ^aoCzrM5I%9KYroi|)ngmhF~Lq~QhgGPh1cE&e_ps|GW5u0 zq}i@u7}0n#u85K_Z6NwNaym=K>cKBUP&;kv8c(4y>rPIiG0kiZzBz)XwG;^0-EkiQ zOlDnWqbPn*nS&y;~1+P0}*K1ev z5-iEqVf#d56^)g@DngA;_EhIsuBe7y+Qa#kJzE{_(s12Y>BR+-eiU5VHn0Q(v$KJM zz>q2EUxj)z>3fXGDl(KZdavf{Ta&ptbwOjqQgo+!GICdFl|Fk#l>D2NBVrVH%7*^jdC z00lt}^K{};d-2ZLK+S~yg-w-@_H#C+6ewA5t`92TUNCMx+Qsivn|<6ft0V+ze=P2+ zh@EzBrl$6tt(>Xt!-qLnkM;z^P(Qy_#i24;{}rQ*hZ$D6>tmz=<*=kKk)JR^;9uU)1u*Rz|~v3Qi?@ zvMBDhUfPmvVFqtro@O%-Suc7y;?)H3+3=T(00A0DA1U$irN4_Y%`cwI&YQ1>D<|i7 zrYNkUDQ;egWOQO2UD~Xf+#~Quq1!hUepSzUUl7&UZ55yWf&@v-u++Pu^Y-l`YsEM} zKV^jj9ArD-08=lAP5t$*`HRtFR=#Vy29>DlY4foZ`(P__PzDD};=wnyiO9249=u;a z3^U^&{@xPD7w|C|&4tX=SR^AvFDX|($ie2A%!i@E{upwS*~wTBB%mGV+;gG##ufJb?JyJa*-`!xuADp791qJ?T=#og2T~-$?CcFx6S(Oh=C+T~FZZUwIT@*VZzgYxN z2ZncdF#~$aYMaW9)%usgF%)qfqluzqBuCML?6Q_PfZw7R7>JvPkf{^g$0n@Ke?kHm zq56vpN??rr52f>sbw`mJM=HHiFTN7RKk854{)oP0b^!8-az|)GX-U(NErx(awgr)* zz*fBjaHcf=`_H@Mb1(Cyl{ER~p`yE&80MBVEWv%RjGtB+4KltJ$Imxg5NjN`_|B@| zPp|0-Nx#SDnjo1{J9(r|l)C^|TGz{ik>AYE!q3fBB9 z%ig5eNbpLII*_F2_p+0X{Na}g9TIc%hXh?MmJ|1!A`fz9A7jYwB{l}M`t=Yw=<;U~ z6XNu3&LSw&8Q%n$1*Y@vja^M4u#D<`RVsF>jt^%% zz16cQNil2r-~tD^TEBd>#4m93UtUN~b46U1=b4ixrfO6?AAiEXDT5`a=HGZ< zq#CTtZnE|5jw}{8y2Wwpd>VxY48Lbs!af=b&*L!yTB@}?%S_PkJI50sVGG^yIR;w7 z&g`auXJRLbj2@HHKpsR|)1}BuSDSCp~fI&0=!c*;aU9&2Hq?Y7Rk%mz9 zlgR;E7%L5uJJZT{o0}D8>Ru;vX@@0>v$kxEFIYpvsGnhVzGv=Khbq zI9?MMCv{LdhKS4&3vzJYPnCw_VldjQ3x!3Y`0VT8WOC*!4Eek~9qzSJL?^R<`0EmZ zucUb`a6hk#|5gBY{uNDTR3~@K>|o)mD`g4W! zSx$b7GeOIO(9*wi;{a6NzQ?-%8>+8phA!1Rl)y9_wH0{g+NSO`2w zcmgP7@QAKH*Gm%WjfuzzPwQsJzkw+*2C({NG{h)p-|QziAFHUaQTV)@@_v?~o-s=! zz*uqoZoSM0ltmFXo5h)Jj$!q=QAJ*3Y*egtgu+apc+ZXAMo8xQr<>KDUo2lmwML)j zb?Nvmm2fw-2oPC#`{9F?ey@-Y=E>TrC*zjno8)n&pPLSuH2YbRg>`M;n;$9PSmAO$ z@64NT60e7@|15X4H#zMEE(gm^xe+%pwu|u2FCeCjLK@k?zcpmP5;UK{hdYrEnSi)GF3v? z;)5XBSVSa#vkz3|%ZGP7aDPa9Nwgw1>meuk5cf~@wb-&mWy!0G5mRk4$ra3<=@s_t zYqU{3h#OqoqW=HOuFd~!NB&>8b^b-*g6-qOXDzdCrLXu0o%shbrqjn%t{j}@_Py0= z{f$oCmF3@iPyL6w2l6&)b^~&4vZQMIcN)|tn?7JW_V^x2@k&v)(t52f);*i;Y}2+g zc?DaLm64+`*yyOFrTck^JPZy{`B)qHca5(=^}uoX#%MRa{JO#-5OLCnRcFJn5-)+}ro`g9P@{1fE z{zO9FWlot~7OH0%`FuQY2`;&g9zXvSR2ek*J3cvXlJ8a*x$J;!NoZOX;5}$uN;`rNM3*Zd&PgKz~y5T=+@8RLW?JmrzRW z$3*Hbs>@lJk7thOIl)8M(VUU4ow(!dP7f2s5t6u&U|k4YtbNVnniNbm_q|=o`RHwR z)69q{?jVbmmGHRPr`Z8_gMP<74-7vfW%UD_Pmi}E%%{R#n~5sC^o6jxI?j+}w0e0}u1r$xBflMeZ8wCYkd=V3 z&*=7YEmhweI;_}j+pPY6fX*E`{qBBfOdBvC?Onal1bZvb2v3sWyJZ_=bhwKsRh7Ct z?M%dc_;fTx0x!>qjZs%k-7E+>zcGuB{dbwr7=Lns-yQaBI!@AetB{0-cqj6Rb-ch_v^ak$#;VKe`)cEhieb=Gzt zt#`I2DMcwlr%#3i-3Xstx}RW$qw0-RqWbM_bq)n>zy0dszA^YnljF|V8E%ko9_=``NKyTLj fd!8k!hiZ)PHor>6;EezN_e4+ISgS as.line() + test_data <- rbind( + rvmf(100, mu = Line(90, 45), k = 10), + rvmf(50, mu = Line(0, 0), k = 20) + ) |> as.line() -ggstereo(grid = TRUE) + - ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) + ggstereo(grid = TRUE) + + ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) -ggstereo(earea = FALSE, centercross = TRUE) + - ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) - } + ggstereo(earea = FALSE, centercross = TRUE) + + ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y)) +} } diff --git a/man/ggstereocontour.Rd b/man/ggstereocontour.Rd index 4fb7c74..fd997e7 100644 --- a/man/ggstereocontour.Rd +++ b/man/ggstereocontour.Rd @@ -54,23 +54,23 @@ Stereonet contouring using ggplot } \examples{ if (require("mapproj")) { -test_data <- rbind( - rvmf(100, mu = Line(90, 45), k = 10), - rvmf(50, mu = Line(0, 0), k = 20) -) |> as.line() + test_data <- rbind( + rvmf(100, mu = Line(90, 45), k = 10), + rvmf(50, mu = Line(0, 0), k = 20) + ) |> as.line() -ggstereo() + - geom_contourf_stereo(gg(test_data)) + - ggplot2::scale_fill_viridis_d(option = "A") + - # guides(fill = guide_colorsteps(barheight = unit(8, "cm"), show.limits = TRUE)) + - geom_contour_stereo(gg(test_data), color = "grey") + - ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y), color = "lightgrey") + - ggframe() + ggstereo() + + geom_contourf_stereo(gg(test_data)) + + ggplot2::scale_fill_viridis_d(option = "A") + + # guides(fill = guide_colorsteps(barheight = unit(8, "cm"), show.limits = TRUE)) + + geom_contour_stereo(gg(test_data), color = "grey") + + ggplot2::geom_point(data = gg(test_data), ggplot2::aes(x = x, y = y), color = "lightgrey") + + ggframe() -ggstereo() + - geom_contourf_stereo(gg(test_data), norm = TRUE, bins = 50, threshold = .1) + - ggplot2::scale_fill_viridis_d(option = "A") - } + ggstereo() + + geom_contourf_stereo(gg(test_data), norm = TRUE, bins = 50, threshold = .1) + + ggplot2::scale_fill_viridis_d(option = "A") +} } \references{ Garcia Portugues, E. (2013). Exact risk improvement of diff --git a/man/or_eigen.Rd b/man/or_eigen.Rd index 885c045..fb97868 100644 --- a/man/or_eigen.Rd +++ b/man/or_eigen.Rd @@ -25,7 +25,7 @@ Decomposition of Orientation Tensor Eigenvectors and Eigenvalues } \examples{ set.seed(1) -mu <- rvmf(n = 1) |> vec2line() +mu <- rvmf(n = 1) |> vec2line() x <- rfb(100, mu = mu, k = 1, A = diag(c(10, 0, 0))) x_eigen <- or_eigen(x) x_eigen diff --git a/man/prepare-ggplot.Rd b/man/prepare-ggplot.Rd index f2909b7..27dcd0f 100644 --- a/man/prepare-ggplot.Rd +++ b/man/prepare-ggplot.Rd @@ -34,21 +34,30 @@ Prepare points and lines for ggplot } \examples{ if (require("mapproj")) { -x <- Plane(120, 85) -ggstereo() + - ggplot2::geom_point(data = gg(x), ggplot2::aes(x, y), color = "red") + - ggplot2::geom_path(data = ggl(x), ggplot2::aes(x, y), color = "red") + x <- Plane(120, 85) + ggstereo() + + ggplot2::geom_point(data = gg(x), ggplot2::aes(x, y), color = "red") + + ggplot2::geom_path(data = ggl(x), ggplot2::aes(x, y), color = "red") -x2 <- Line(120, 5) -ggstereo() + - ggplot2::geom_point(data = gg(x2), ggplot2::aes(x, y), color = "darkgreen") + - ggplot2::geom_path(data = ggl(x2, d = 8), ggplot2::aes(x, y, group = group), color = "darkgreen") - -x3 <- Plane(137, 71) -ggstereo() + - ggplot2::geom_point(data = gg(x3), ggplot2::aes(x, y), color = "darkgreen") + - ggplot2::geom_path(data = ggl(x3, d = 90), ggplot2::aes(x, y, group = group), color = "darkgreen", lwd = 1) + - ggplot2::geom_path(data = ggl(x3, d = 90 + 11), ggplot2::aes(x, y, group = group, color = 'sde <90')) + - ggplot2::geom_path(data = ggl(x3, d = 90 - 11), ggplot2::aes(x, y, group = group, color = 'sde >90')) - } + x2 <- Line(120, 5) + ggstereo() + + ggplot2::geom_point(data = gg(x2), ggplot2::aes(x, y), color = "darkgreen") + + ggplot2::geom_path(data = ggl(x2, d = 8), ggplot2::aes(x, y, group = group), color = "darkgreen") + + x3 <- Plane(137, 71) + ggstereo() + + ggplot2::geom_point(data = gg(x3), ggplot2::aes(x, y), color = "darkgreen") + + ggplot2::geom_path( + data = ggl(x3, d = 90), + ggplot2::aes(x, y, group = group), color = "darkgreen", lwd = 1 + ) + + ggplot2::geom_path( + data = ggl(x3, d = 90 + 11), + ggplot2::aes(x, y, group = group, color = "sde <90") + ) + + ggplot2::geom_path( + data = ggl(x3, d = 90 - 11), + ggplot2::aes(x, y, group = group, color = "sde >90") + ) +} } diff --git a/man/stats.Rd b/man/stats.Rd index 494b916..e95541e 100644 --- a/man/stats.Rd +++ b/man/stats.Rd @@ -49,7 +49,7 @@ Statistical estimators of the distribution of a set of vectors if otherwise). For enough large sample it approaches the angular standard deviation (\code{"csd"}) of the Fisher statistics. -\code{v_rdegree} returns the degree of preferred orientation of vectors (in \%). +\code{v_rdegree} returns the degree of preferred orientation of vectors, range: (0, 1). \code{v_sde} returns the spherical standard error (numeric). If the number of data is less than 25, if will print a additional message, that the output @@ -61,7 +61,9 @@ if otherwise). The \eqn{100(1-\alpha)\%} confidence interval is than given by \e \code{estimate_k} returns the estimated concentration of the von Mises-Fisher distribution \eqn{\kappa} (after Sra, 2011). } \examples{ +set.seed(1234) x <- rvmf(100, mu = Line(120, 50), k = 5) +v_mean(x) v_var(x) v_delta(x) v_rdegree(x) @@ -71,7 +73,7 @@ estimate_k(x) fisher_statistics(x) #' weights: -x2 <- Line(c(0, 0), c(0, 90)) +x2 <- Line(c(0, 0), c(0, 90)) v_mean(x2) v_mean(x2, w = c(1, 2)) v_var(x2) diff --git a/vignettes/Intro.Rmd b/vignettes/Intro.Rmd index eca0b0f..0a43fb4 100644 --- a/vignettes/Intro.Rmd +++ b/vignettes/Intro.Rmd @@ -55,12 +55,12 @@ ggstereo() + ```{r stat, message=FALSE,warning=FALSE} quality <- runif(nrow(lines), min = 1, max = 45) # assigning a random quality score to the data (can be replaced with real data) -lines_mean <- v_mean(lines, w = 1/quality) -lines_delta <- v_delta(lines, w = 1/quality) +lines_mean <- v_mean(lines, w = 1 / quality) +lines_delta <- v_delta(lines, w = 1 / quality) ggstereo() + geom_point(data = gg(lines, quality), aes(x, y, size = quality)) + - scale_size('Quality', range = c(3, .1)) + + scale_size("Quality", range = c(3, .1)) + geom_path(data = ggl(lines_mean, d = lines_delta), aes(x, y, color = "Std"), lwd = .1) + geom_point(data = gg(lines_mean), aes(x, y, color = "Mean"), size = 5, shape = 17) + labs(title = "Example data", color = NULL) @@ -86,23 +86,23 @@ VollmerPlot(lines, add = TRUE, col = "dodgerblue", pch = 19, cex = 2) ```{r density, message=FALSE,warning=FALSE} ggstereo() + geom_contourf_stereo(gg(planes), show.legend = TRUE, norm = TRUE) + - scale_fill_viridis_d('Density') + + scale_fill_viridis_d("Density") + geom_point(data = gg(planes), aes(x, y)) + labs(title = "Example data", color = NULL) ``` ## Facets -```{r density, message=FALSE,warning=FALSE} +```{r facets , message=FALSE,warning=FALSE} area_l <- LETTERS[sample.int(3, nrow(lines), replace = TRUE)] area_p <- LETTERS[sample.int(3, nrow(planes), replace = TRUE)] -lines_df <- gg(lines, area=area_l) +lines_df <- gg(lines, area = area_l) planes_df <- ggl(planes, area = area_p) -ggstereo(data = lines_df, aes(x=x, y=y, color = area)) + - geom_path(data = planes_df, aes(x=x, y=y, group = group), alpha = .25, color = 'grey') + +ggstereo(data = lines_df, aes(x = x, y = y, color = area)) + + geom_path(data = planes_df, aes(x = x, y = y, group = group), alpha = .25, color = "grey") + geom_point() + facet_wrap(vars(area)) + labs(title = "Example data", color = NULL) -``` \ No newline at end of file +```