From d151e6a8d6fcd3d80f9aed84dce3af9e621f0951 Mon Sep 17 00:00:00 2001 From: Andrew Barbour Date: Thu, 14 May 2015 12:17:56 -0700 Subject: [PATCH] Added mogi; appears to be functioning properly; needs testing though. --- DESCRIPTION | 11 +- NAMESPACE | 5 +- R/deform-utilities.R | 1 + R/{epsg.R => spatial.R} | 46 +++++++- R/{mag-dist.R => strain_magdist_scaling.R} | 0 R/volume-sources.R | 128 +++++++++++++++++++++ man/Simple-deformation.Rd | 2 +- man/cart2pol.Rd | 39 +++++++ man/deform.Rd | 2 +- man/deform_hull.Rd | 2 +- man/degreelen.Rd | 2 +- man/epsg.Rd | 4 +- man/erf.Rd | 2 +- man/hydraulic_diffusivity.Rd | 2 +- man/magnitude-distance.Rd | 4 +- man/mogi.Rd | 105 +++++++++++++++++ man/segall85.Rd | 2 +- 17 files changed, 340 insertions(+), 17 deletions(-) rename R/{epsg.R => spatial.R} (69%) rename R/{mag-dist.R => strain_magdist_scaling.R} (100%) create mode 100644 R/volume-sources.R create mode 100644 man/cart2pol.Rd create mode 100644 man/mogi.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 81168b7..b7322e8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,13 +1,16 @@ Package: deform Type: Package Title: Crustal deformation modeling tools in R -Version: 0.0-5 -Date: 2015-04-06 +Version: 0.0-6.9999 +Date: 2015-05-14 Author: Andrew J. Barbour; some of these tools were - adapted from Matlab codes by Brendan Crowell (UW), + adapted from Matlab(tm) codes by Brendan Crowell (UW), and Francois Beauducel (IPGP). Maintainer: A.J.Barbour -Description: Crustal deformation modeling tools. +Description: Crustal deformation modeling tools. Currently this + code supports calculation of surface deformation associated with + (1) a depleting layer using a Segall-type model, and + (2) a volume source using a Mogi-type model. License: GPL Depends: rgdal (>= 0.9-2) diff --git a/NAMESPACE b/NAMESPACE index eeb7aa3..fe5fb57 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -# Generated by roxygen2 (4.1.0): do not edit by hand +# Generated by roxygen2 (4.1.1): do not edit by hand export(.setleft) export(.surface_g) @@ -6,6 +6,7 @@ export(D_project) export(Tilt) export(Uniaxial_extension) export(agency.root) +export(cart2pol) export(deform_hull) export(degreelen) export(dynamic) @@ -16,6 +17,8 @@ export(getEPSG) export(hydraulic_diffusivity) export(ierf) export(ierfc) +export(mogi.pressure) +export(mogi.volume) export(rigid_hydraulic_diffusivity) export(segall85) export(static) diff --git a/R/deform-utilities.R b/R/deform-utilities.R index 033a36e..eb803d8 100644 --- a/R/deform-utilities.R +++ b/R/deform-utilities.R @@ -1,3 +1,4 @@ + #' Error functions #' @param x numeric; #' @export diff --git a/R/epsg.R b/R/spatial.R similarity index 69% rename from R/epsg.R rename to R/spatial.R index 2f86be8..c6d25da 100644 --- a/R/epsg.R +++ b/R/spatial.R @@ -1,3 +1,48 @@ +#' Transform cartesian coordinates to polar +#' @param x numeric vector, or numeric matrix (see \code{y}) +#' @param y numeric vector; if missing, then \code{x} must have at least +#' two columns: the first to represent \code{x} and the second to represent \code{y} +#' @param z numeric vector +#' @param dist.type character; the type of distance calculation +#' @param radians logical; should the output angle be in radians? (\code{FALSE} returns in degrees) +#' @author modified from function in beadarrayMSV +#' @export +#' @examples +#' +#' cart2pol(1:10, 1:10) +#' all.equal(cart2pol(1:10, 1:10), cart2pol(cbind(1:10, 1:10))) +#' +#' # Change radius computation +#' cart2pol(1:10, 1:10, 1:10) # assumes 'euclidean' +#' cart2pol(1:10, 1:10, 1:10, 'manhattan') +#' +#' # use degrees +#' cart2pol(1:10, 1:10, radians=FALSE) +#' +cart2pol <- function(x, y, z = NULL, dist.type = c('euclidean','manhattan'), radians=TRUE){ + pol <- list() + dist.type <- match.arg(dist.type) + if (missing(y)){ + X <- as.matrix(x) + stopifnot(ncol(X) >= 2) + y <- X[,2] + x <- X[,1] + } + x <- as.vector(x) + y <- as.vector(y) + r. <- switch(dist.type, + euclidean = sqrt(x^2 + y^2), + manhattan = x + y + ) + th. <- atan2(y, x) + if (!radians) th. <- th. * 180 / pi + rth <- cbind(Theta = th., Radius = r.) + if (!is.null(z)){ + rth <- cbind(rth, Z = as.vector(z)) + } + return(rth) +} + #' Functions to work with EPSG projection specifications #' @name epsg #' @aliases EPSG @@ -96,4 +141,3 @@ D_project <- function(Df, coords=names(Df)[1:2], p4, p4old, verbose=TRUE, ...){ stopifnot(!inherits(toret,"try-error")) return(toret) } -## \ No newline at end of file diff --git a/R/mag-dist.R b/R/strain_magdist_scaling.R similarity index 100% rename from R/mag-dist.R rename to R/strain_magdist_scaling.R diff --git a/R/volume-sources.R b/R/volume-sources.R new file mode 100644 index 0000000..91dc6f8 --- /dev/null +++ b/R/volume-sources.R @@ -0,0 +1,128 @@ +#' Surface deformation from Mogi-type volume sources +#' +#' Modified from [1]: Computes radial and vertical displacements, ground tilt, and +#' radial and tangential strain on the surface, at radial distances +#' from the top of the source, due to a hydrostatic pressure inside a +#' spherical cavity at some depth in a homogeneous elastic halfspace +#' (center of dilatation). See Mogi (1958). +#' +#' @details Strain is reckoned positive for extension. For insight into what a positive +#' tilt is, see \code{\link{Tilt}}. +#' +#' Units are shown next to the argument descriptions. +#' +#' Typical moduli are in the range of 1-10 GPa for shear +#' modulus, and 10-100 GPa for Young's modulus. +#' \emph{Note that one gigapascal (GPa) is +#' equal to 10^9 pascal (Pa).} +#' +#' @name mogi +#' @aliases Mogi +#' +#' @param surface.distance numeric; [m] the radial distance (positive away) from the center of the spherical cavity, projected onto the surface +#' @param sphere.depth numeric; [m] depth of the center of the spherical cavity from the surface (positive downward) +#' @param Volume.change numeric; [m^3] volumetric change of the spherical cavity (inflation positive) +#' @param sphere.radius numeric; [m] radius of the spherical cavity that's inflating or deflating +#' @param Pressure.change numeric; [Pa] hydrostatic pressure change in the spherical cavity (inflation positive) +#' @param shear.modulus numeric; [Pa] elastic shear modulus of the surrounding material, also known as 'rigidity' or Lame''s constant +#' @param Youngs.modulus numeric; [Pa] Young's modulus of elasticity of the surrounding material +#' @param nu numeric; [0-1] Poisson's ratio of the surrounding material +#' @param verbose logical; should messages be given? +#' +#' @references [1] \url{http://www.ipgp.fr/~beaudu/matlab.html#Mogi} +#' @references Mogi, K. (1958), Relations between the eruptions of various volcanoes and the +#' deformations of the ground surfaces around them, \emph{Bull. Earthquake Res. +#' Inst. Univ. Tokyo}, 36, 99-134. (\url{http://hdl.handle.net/2261/11909}) +#' +#' @author Original Matlab code 'mogi.m' written by Francois Beauducel; AJ Barbour ported to R code. +#' @family volume-sources +#' @seealso \code{\link{Uniaxial_extension}} and \code{\link{Tilt}} +#' @examples +#' library(RColorBrewer) +#' seq.palfun <- function(n) brewer.pal(n, 'Reds') +#' div.palfun <- function(n) brewer.pal(n, 'Spectral') +#' # http://www.ipgp.fr/~beaudu/images/mogi_example.png +#' +#' x <- y <- unique(sort(c(seq(-2,2,by=0.1), seq(-0.1,0.1,by=0.01)))) +#' theta_r <- cart2pol(expand.grid(x, y)) +#' r <- theta_r[,'Radius'] +#' +#' # calculate the deformation specifying a positive rate of volume-change (inflation) +#' MV <- mogi.volume(r, 1, 1e-5, verbose=TRUE) +#' +#' filled.contour(x, y, matrix((MV[,'Ur']), length(x)), asp=1, main='Ur', nlevels = 11, color.palette=div.palfun) +#' filled.contour(x, y, matrix((MV[,'Uz']), length(x)), asp=1, main='Uz', color.palette=div.palfun, levels=seq(0,2.5e-6, length.out = 11)) +#' +#' filled.contour(x, y, matrix((MV[,'Ett']), length(x)), levels=seq(0,2.5e-6,length.out=9), asp=1, main='Ett', color.palette=seq.palfun) +#' filled.contour(x, y, matrix((MV[,'Err']), length(x)), zlim=3e-6*c(-1,1), asp=1, main='Err', nlevels = 11) +#' +#' # There should be a null in the center of the tilt field +#' filled.contour(x, y, matrix((MV[,'Tilt']), length(x)), asp=1, main='Tilt', levels=seq(0,2.5e-6,length.out=9), color.palette=seq.palfun) +#' +#' # Calculate the undrained-to-drained effect: +#' # first calculate for nu=1/3 +#' MVud <- mogi.volume(r, 1, 1e-5, nu=1/3, verbose=TRUE) +#' # then subtract the 1/4 result +#' Resp.ud <- matrix((MV[,'Uz']), length(x)) - matrix((MVud[,'Uz']), length(x)) +#' filled.contour(x, y, log10(Resp.ud), asp=1, main='Undrained response: Uz', color.palette=seq.palfun, levels=seq(-8,-6.5,length.out=9)) +#' # +#' +#' Using values comparable to the previous example, but specifying pressure changes +#' MP <- mogi.pressure(r, 1, 0.01, 1e10, shear.modulus=3.08e9) +#' MP2 <- mogi.pressure(r, 1, 0.01, 1e10, Youngs.modulus=10e9) +#' +NULL + +#' @rdname mogi +#' @export +mogi.volume <- function(surface.distance, sphere.depth, Volume.change, ...){ + .mogi_check(surface.distance, sphere.depth) + Source <- Volume.change / pi + .mogi_calc(Source, surface.distance, sphere.depth, ...) +} + +#' @rdname mogi +#' @export +mogi.pressure <- function(surface.distance, sphere.depth, sphere.radius, Pressure.change, shear.modulus, Youngs.modulus, nu=1/4, ...){ + .mogi_check(surface.distance, sphere.depth, sphere.radius) + if (missing(shear.modulus)){ + if (missing(Youngs.modulus)) stop('Youngs.modulus cannot be missing if shear.modulus is missing') + shear.modulus <- Youngs.modulus / (2 * (1 + nu)) + } + Source <- sphere.radius ^ 3 * Pressure.change / shear.modulus + .mogi_calc(Source, surface.distance, sphere.depth, ...) +} + +#' @rdname mogi +.mogi_check <- function(surface.distance, sphere.depth, sphere.radius=NULL){ + stopifnot(all(surface.distance >= 0)) + stopifnot(all(sphere.depth >= 0)) + if (!is.null(sphere.radius)){ + stopifnot(all(sphere.radius >= 0)) + if (max(sphere.radius) / min(sphere.depth) > 1/10) warning('inaccurate results if sphere.depth is not much greater than sphere.radius') + } +} + +#' @rdname mogi +.mogi_calc <- function(src, surface.distance, sphere.depth, nu=1/4, verbose=FALSE){ + + stopifnot(all(nu>=0 & nu <=1)) + + # radial distance from source + radial.distance <- sqrt(sphere.depth^2 + surface.distance^2) + if (verbose) message("Poissons ratio: ", paste(nu)) + C. <- (1 - nu) * src + if (verbose) message("constant(s): ", paste(signif(C.))) + + # tangential horizontal (hoop) strain, radial horizontal displacement, and vertical displacement + Ett <- C. / radial.distance^3 + Ur <- surface.distance * Ett + Uz <- sphere.depth * Ett + + # tilt, and radial horizontal strain + gradUz <- 3 * C. * sphere.depth * surface.distance / radial.distance^5 + Err <- C. * (sphere.depth^2 - 2 * surface.distance^2) / radial.distance^5 + + data.frame(radial.distance, Ur, Uz, Ett, Err, Tilt=gradUz) + +} diff --git a/man/Simple-deformation.Rd b/man/Simple-deformation.Rd index bcc1c01..7167264 100644 --- a/man/Simple-deformation.Rd +++ b/man/Simple-deformation.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/segall85.R \name{Simple-deformation} \alias{.setleft} diff --git a/man/cart2pol.Rd b/man/cart2pol.Rd new file mode 100644 index 0000000..1f6ebc6 --- /dev/null +++ b/man/cart2pol.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/spatial.R +\name{cart2pol} +\alias{cart2pol} +\title{Transform cartesian coordinates to polar} +\usage{ +cart2pol(x, y, z = NULL, dist.type = c("euclidean", "manhattan"), + radians = TRUE) +} +\arguments{ +\item{x}{numeric vector, or numeric matrix (see \code{y})} + +\item{y}{numeric vector; if missing, then \code{x} must have at least +two columns: the first to represent \code{x} and the second to represent \code{y}} + +\item{z}{numeric vector} + +\item{dist.type}{character; the type of distance calculation} + +\item{radians}{logical; should the output angle be in radians? (\code{FALSE} returns in degrees)} +} +\description{ +Transform cartesian coordinates to polar +} +\examples{ +cart2pol(1:10, 1:10) +all.equal(cart2pol(1:10, 1:10), cart2pol(cbind(1:10, 1:10))) + +# Change radius computation +cart2pol(1:10, 1:10, 1:10) # assumes 'euclidean' +cart2pol(1:10, 1:10, 1:10, 'manhattan') + +# use degrees +cart2pol(1:10, 1:10, radians=FALSE) +} +\author{ +modified from function in beadarrayMSV +} + diff --git a/man/deform.Rd b/man/deform.Rd index 0630dd5..d75f092 100644 --- a/man/deform.Rd +++ b/man/deform.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/deform-package.R \docType{package} \name{deform} diff --git a/man/deform_hull.Rd b/man/deform_hull.Rd index 09d6027..a7f78bb 100644 --- a/man/deform_hull.Rd +++ b/man/deform_hull.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/seismicity.R \name{deform_hull} \alias{deform_hull} diff --git a/man/degreelen.Rd b/man/degreelen.Rd index 4d24bfe..c775846 100644 --- a/man/degreelen.Rd +++ b/man/degreelen.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/degreelen.R \name{degreelen} \alias{degreelen} diff --git a/man/epsg.Rd b/man/epsg.Rd index a1b34ab..5128842 100644 --- a/man/epsg.Rd +++ b/man/epsg.Rd @@ -1,5 +1,5 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand -% Please edit documentation in R/epsg.R +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/spatial.R \name{epsg} \alias{D_project} \alias{EPSG} diff --git a/man/erf.Rd b/man/erf.Rd index 1df0639..b18f258 100644 --- a/man/erf.Rd +++ b/man/erf.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/deform-utilities.R \name{erf} \alias{erf} diff --git a/man/hydraulic_diffusivity.Rd b/man/hydraulic_diffusivity.Rd index b998247..762e6db 100644 --- a/man/hydraulic_diffusivity.Rd +++ b/man/hydraulic_diffusivity.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/diffusivity.R \name{hydraulic_diffusivity} \alias{hydraulic_diffusivity} diff --git a/man/magnitude-distance.Rd b/man/magnitude-distance.Rd index 705dfda..3e86688 100644 --- a/man/magnitude-distance.Rd +++ b/man/magnitude-distance.Rd @@ -1,5 +1,5 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand -% Please edit documentation in R/mag-dist.R +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/strain_magdist_scaling.R \name{magnitude-distance} \alias{dynamic} \alias{magnitude-distance} diff --git a/man/mogi.Rd b/man/mogi.Rd new file mode 100644 index 0000000..77065c6 --- /dev/null +++ b/man/mogi.Rd @@ -0,0 +1,105 @@ +% Generated by roxygen2 (4.1.1): do not edit by hand +% Please edit documentation in R/volume-sources.R +\name{mogi} +\alias{.mogi_calc} +\alias{.mogi_check} +\alias{Mogi} +\alias{mogi} +\alias{mogi.pressure} +\alias{mogi.volume} +\title{Surface deformation from Mogi-type volume sources} +\usage{ +mogi.volume(surface.distance, sphere.depth, Volume.change, ...) + +mogi.pressure(surface.distance, sphere.depth, sphere.radius, Pressure.change, + shear.modulus, Youngs.modulus, nu = 1/4, ...) + +.mogi_check(surface.distance, sphere.depth, sphere.radius = NULL) + +.mogi_calc(src, surface.distance, sphere.depth, nu = 1/4, verbose = FALSE) +} +\arguments{ +\item{surface.distance}{numeric; [m] the radial distance (positive away) from the center of the spherical cavity, projected onto the surface} + +\item{sphere.depth}{numeric; [m] depth of the center of the spherical cavity from the surface (positive downward)} + +\item{Volume.change}{numeric; [m^3] volumetric change of the spherical cavity (inflation positive)} + +\item{sphere.radius}{numeric; [m] radius of the spherical cavity that's inflating or deflating} + +\item{Pressure.change}{numeric; [Pa] hydrostatic pressure change in the spherical cavity (inflation positive)} + +\item{shear.modulus}{numeric; [Pa] elastic shear modulus of the surrounding material, also known as 'rigidity' or Lame''s constant} + +\item{Youngs.modulus}{numeric; [Pa] Young's modulus of elasticity of the surrounding material} + +\item{nu}{numeric; [0-1] Poisson's ratio of the surrounding material} + +\item{verbose}{logical; should messages be given?} +} +\description{ +Modified from [1]: Computes radial and vertical displacements, ground tilt, and +radial and tangential strain on the surface, at radial distances +from the top of the source, due to a hydrostatic pressure inside a +spherical cavity at some depth in a homogeneous elastic halfspace +(center of dilatation). See Mogi (1958). +} +\details{ +Strain is reckoned positive for extension. For insight into what a positive +tilt is, see \code{\link{Tilt}}. + +Units are shown next to the argument descriptions. + +Typical moduli are in the range of 1-10 GPa for shear +modulus, and 10-100 GPa for Young's modulus. +\emph{Note that one gigapascal (GPa) is +equal to 10^9 pascal (Pa).} +} +\examples{ +library(RColorBrewer) +seq.palfun <- function(n) brewer.pal(n, 'Reds') +div.palfun <- function(n) brewer.pal(n, 'Spectral') +# http://www.ipgp.fr/~beaudu/images/mogi_example.png + +x <- y <- unique(sort(c(seq(-2,2,by=0.1), seq(-0.1,0.1,by=0.01)))) +theta_r <- cart2pol(expand.grid(x, y)) +r <- theta_r[,'Radius'] + +# calculate the deformation specifying a positive rate of volume-change (inflation) +MV <- mogi.volume(r, 1, 1e-5, verbose=TRUE) + +filled.contour(x, y, matrix((MV[,'Ur']), length(x)), asp=1, main='Ur', nlevels = 11, color.palette=div.palfun) +filled.contour(x, y, matrix((MV[,'Uz']), length(x)), asp=1, main='Uz', color.palette=div.palfun, levels=seq(0,2.5e-6, length.out = 11)) + +filled.contour(x, y, matrix((MV[,'Ett']), length(x)), levels=seq(0,2.5e-6,length.out=9), asp=1, main='Ett', color.palette=seq.palfun) +filled.contour(x, y, matrix((MV[,'Err']), length(x)), zlim=3e-6*c(-1,1), asp=1, main='Err', nlevels = 11) + +# There should be a null in the center of the tilt field +filled.contour(x, y, matrix((MV[,'Tilt']), length(x)), asp=1, main='Tilt', levels=seq(0,2.5e-6,length.out=9), color.palette=seq.palfun) + +# Calculate the undrained-to-drained effect: +# first calculate for nu=1/3 +MVud <- mogi.volume(r, 1, 1e-5, nu=1/3, verbose=TRUE) +# then subtract the 1/4 result +Resp.ud <- matrix((MV[,'Uz']), length(x)) - matrix((MVud[,'Uz']), length(x)) +filled.contour(x, y, log10(Resp.ud), asp=1, main='Undrained response: Uz', color.palette=seq.palfun, levels=seq(-8,-6.5,length.out=9)) +# + +Using values comparable to the previous example, but specifying pressure changes +MP <- mogi.pressure(r, 1, 0.01, 1e10, shear.modulus=3.08e9) +MP2 <- mogi.pressure(r, 1, 0.01, 1e10, Youngs.modulus=10e9) +} +\author{ +Original Matlab code 'mogi.m' written by Francois Beauducel; AJ Barbour ported to R code. +} +\references{ +[1] \url{http://www.ipgp.fr/~beaudu/matlab.html#Mogi} + +Mogi, K. (1958), Relations between the eruptions of various volcanoes and the + deformations of the ground surfaces around them, \emph{Bull. Earthquake Res. + Inst. Univ. Tokyo}, 36, 99-134. (\url{http://hdl.handle.net/2261/11909}) +} +\seealso{ +\code{\link{Uniaxial_extension}} and \code{\link{Tilt}} +} + diff --git a/man/segall85.Rd b/man/segall85.Rd index 0fa3495..e7f449b 100644 --- a/man/segall85.Rd +++ b/man/segall85.Rd @@ -1,4 +1,4 @@ -% Generated by roxygen2 (4.1.0): do not edit by hand +% Generated by roxygen2 (4.1.1): do not edit by hand % Please edit documentation in R/segall85.R \name{segall85} \alias{.surface_g}