Gravitational and magnetic attraction of 3-D vertical rectangular prisms
gravmagsubs is a software package for the R language (R Core Team 2021) for forward modeling of gravity and magnetic anomalies from collections of 3-D right rectangular prisms. The package implements computational approaches developed by (Plouff 1975a, 1975b, 1976) for forward modeling of gravity and magnetic anomalies.
This package is dependent on the following R packages:
Rcpp
- Seamless R and C++ Integration (Eddelbuettel and
Francois 2011)
citation("gravmagsubs")
##
## To cite gravmagsubs in publications, please use:
##
## Cronkite-Ratcliff, C., Phelps, G., Scheirer, D., 2022, gravmagsubs:
## Gravitational and magnetic attraction of 3-D vertical rectangular
## prisms, U.S. Geological Survey software release, R package, version
## 1.0, https://doi.org/10.5066/P9HBSPRM
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## author = {Collin Cronkite-Ratcliff and Geoffrey Phelps and Daniel Scheirer},
## title = {gravmagsubs: Gravitational and magnetic attraction of 3-D vertical rectangular prisms},
## publisher = {U.S. Geological Survey},
## address = {Reston, VA},
## version = {1.0},
## note = {R package},
## institution = {U.S. Geological Survey},
## year = {2022},
## doi = {10.5066/P9HBSPRM},
## url = {https://doi.org/10.5066/P9HBSPRM},
## }
It is assumed that users already have installed R (https://www.r-project.org) (R Core Team 2021).
Pandoc (https://pandoc.org) is required to create the package vignettes. Pandoc can be obtained by installing the Rstudio IDE (https://posit.co/), which comes with Pandoc. If you do not install the Rstudio IDE, you will need to install Pandoc individually (https://pandoc.org/installing.html).
Use the remotes
package (Csárdi et al. 2022) to install
the package directly from code.usgs.gov. When running
remotes::install_gitlab
, you may be asked:
# If you've never installed the "remotes" package:
install.packages("remotes")
#Installing the package from gitlab
::install_gitlab("gmegsc/gravmagsubs",
remoteshost = "code.usgs.gov",
build_opts = c("--no-resave-data",
"--no-manual"),
dependencies = TRUE,
build_vignettes = TRUE)
These packages have more recent versions available.
It is recommended to update all of them.
Which would you like to update?
1: All
2: CRAN packages only
3: None
: 3 Enter one or more numbers, or an empty line to skip updates
Choose “3. None” by typing the number 3, and update all packages after the installation.
Regularly, it is a good idea to update ALL your packages in R. If using RStudio, this is quite easy, there’s an Update button in the “Packages” tab. This checks CRAN for updates. It is a good idea to click this update regularly.
After installation, in order to use the package, load the package into the R session:
library(gravmagsubs)
For help on how to use the package, use the following command to view the documentation for the package:
?gravmagsubs
or
help(gravmagsubs)
To view the documentation for the functions
rectprismgrav()
or rectprismmag()
:
?rectprismgrav ?rectprismmag
or
help(rectprismgrav)
help(rectprismmag)
The gravmagsubs
package provides tools for computing the
gravitational and magnetic anomalies generated by 3-D vertical
rectangular prisms at specific observation points. The package consists
of two functions:
rectprismgrav()
: Computes the gravitational
attraction of 3-D right rectangular prisms.
rectprismmag()
: Computes the magnetic effect of 3-D
right rectangular prisms.
Each function can compute the total anomaly of a series of N prisms at M observation points, returned as a matrix of M rows and 1 column.
Each function also has a logical flag bycell
(default
FALSE
). If bycell=TRUE
, the function returns
the contribution from each individual prism. In this case, the function
will return a matrix of M rows by N columns.
The rectprismgrav()
function calculates the gravity
anomaly from a series of 3-D rectangular prisms with vertical sides
using the algorithm of Plouff (1975a).
To calculate the gravity anomaly of a single prism at a single point:
#########################################################
## gravity anomaly of a single prism at a single point ##
# location of the point where the gravity anomaly will be calculated
<- data.frame(x=0, y=0, z=0)
gravstation
# the rectangular prism is defined by its six edges (distance in km)
<- data.frame(xmin=-5, xmax=5,
prism1 ymin=-5, ymax=5,
zmin=-10, zmax=-5)
# density contrast in g/cc
<- 0.3
drho
<- rectprismgrav(gravstation$x, gravstation$y, gravstation$z,
gravanom $xmin, prism1$xmax,
prism1$ymin, prism1$ymax,
prism1$zmin, prism1$zmax, drho)
prism1
# gravity anomaly (mGal)
print(gravanom)
## [,1]
## [1,] 13.17689
To calculate the gravity anomaly at multiple stations along a 2-D cross-section:
#########################################################
## gravity anomaly from a 2-D section ##
# gravity station x-axis locations
<- data.frame(X = seq(-25, by=.1, length=500))
grav.calc
# density contrast in g/cc
<- -0.67
rho
# 2-D section consists of a layered model with density
# contrast extending to 1 km on the left side and 1.5 km
# on the right side
<- 1.5
h1 <- 1
h2
# create the 2-D section, with prisms measuring 1 km
# horizontally and 100 m vertically.
<- 1
prism.width.h <- 0.1
prism.width.v <- seq(-499.5, 499.5, by=prism.width.h)
X1 <- seq(-3.95, -0.05, by=prism.width.v)
Z1 <- expand.grid(xcenter=X1, Y=0, zcenter=Z1)
XZr
# define the density contrast as 0 or -0.67 g/cc.
$density <- 0
XZr$density[XZr$zcenter > -h1] <- rho
XZr$density[XZr$zcenter > -(h1+h2) & XZr$xcenter > 0] <- rho
XZr
# the y- and z-locations are both zero
$Y <- 0
grav.calc$Z <- 0
grav.calc
# Calculate the gravity anomaly.
# Extend the prisms +/-500 km in the Y direction
# to approximate a perfect 2-D model.
$Ggms <- rectprismgrav(grav.calc$X, grav.calc$Y, grav.calc$Z,
grav.calc$xcenter - prism.width.h/2,
XZr$xcenter + prism.width.h/2,
XZr$Y - 500, XZr$Y + 500,
XZr$zcenter - prism.width.v/2,
XZr$zcenter + prism.width.v/2, XZr$density)
XZr
library(ggplot2)
ggplot() +
geom_line(data = grav.calc, linewidth=2, aes(x = X, y = Ggms)) +
labs(title = "", x = "Distance [km]", y = "Gravity anomaly (mGal)") +
theme(panel.background=element_rect(colour="black", fill="white"),
panel.grid = element_blank(),
plot.title = element_text(face = "bold",hjust = 0.5, size=16 ),
aspect.ratio = 1/4)
Use the bycell=TRUE
option and unit density (1 g/cc) to
obtain the sensitivity matrix.
# Use unit density contrast to obtain sensitivity matrix
$density1 <- 1
XZr<- rectprismgrav(grav.calc$X, grav.calc$Y, grav.calc$Z,
XZr.sensmat $xcenter - prism.width.h/2,
XZr$xcenter + prism.width.h/2,
XZr$Y - 500, XZr$Y + 500,
XZr$zcenter - prism.width.v/2,
XZr$zcenter + prism.width.v/2, XZr$density1,
XZrbycell=TRUE)
Now define two new cross-sections with different density values. The
gravity anomaly of these new cross-sections can be obtained by
multiplying the new density values with the sensitivity matrix, rather
than re-running rectprismgrav()
.
# Define different density contrast values.
$density2 <- 0
XZr$density2[XZr$density == rho] <- -0.9
XZr$density3 <- 0
XZr$density3[XZr$density == rho] <- -0.4
XZr
# Use matrix multiplication to obtain new gravity anomaly.
<- XZr.sensmat %*% as.matrix(XZr[,c("density", "density2", "density3")])
gravanom.models
$gmod.orig <- gravanom.models[,1]
grav.calc$gmod.low <- gravanom.models[,2]
grav.calc$gmod.high <- gravanom.models[,3]
grav.calc
ggplot() +
geom_line(data = grav.calc, linewidth=2,
aes(x = X, y = gmod.low, colour = "-0.9 g/cc")) +
geom_line(data = grav.calc, linewidth=2,
aes(x = X, y = gmod.orig, colour = "-0.67 g/cc")) +
geom_line(data = grav.calc, linewidth=2,
aes(x = X, y = gmod.high, colour = "-0.4 g/cc")) +
scale_colour_manual(name="Density contrast", values=c("blue", "black", "red")) +
labs(title = "", x = "Distance [km]", y = "Gravity anomaly (mGal)") +
theme(panel.background=element_rect(colour="black", fill="white"),
legend.key = element_rect(fill = "white"),
legend.position = "right",
panel.grid = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5, size=16 ),
aspect.ratio = 1/4)
To calculate the gravity anomaly at a 2-D grid of stations located above a 3-D grid of prisms with different density values:
#########################################################
## gravity anomaly from 3-D gridded source model ##
# first, create an array of stations
<- seq(-5, 5, by=0.1)
xstations <- seq(-5, 5, by=0.1)
ystations <- expand.grid(xstations, ystations)
D3.stations $Zstation <- 0
D3.stationsnames(D3.stations) <- c("Xstation", "Ystation", "Zstation")
# second, create an array of prisms
<- 0.1
width <- width / 2
half.width
<- seq(-2 + half.width, 2 - half.width, by=width)
xsource <- seq(-2 + half.width, 2 - half.width, by=width)
ysource <- seq(-3 + half.width, -1 - half.width, by=width)
zsource
<- expand.grid(xsource, ysource, zsource)
D3.source names(D3.source) <- c("xcenter", "ycenter", "zcenter")
# density varies along the y-axis
$density <- 0.1 * D3.source$ycenter
D3.source
<- width / 2
half.width
<- rectprismgrav(D3.stations$Xstation,
D3.gravanom $Ystation,
D3.stations$Zstation,
D3.stations$xcenter - half.width,
D3.source$xcenter + half.width,
D3.source$ycenter - half.width,
D3.source$ycenter + half.width,
D3.source$zcenter - half.width,
D3.source$zcenter + half.width,
D3.source$density,
D3.sourcebycell=FALSE)
library(scales)
# create data frame
<- data.frame(x=D3.stations$Xstation, y=D3.stations$Ystation,
d3.st val1 = D3.gravanom)
# plot with custom color scale
ggplot() +
geom_raster(data = d3.st, aes(x = x, y = y, fill = val1)) +
scale_fill_gradientn(
colours = c("darkblue", "blue", "skyblue", "white", "tomato", "red", "darkred"),
values = rescale(c(-1.5, 0, 1.5)),
breaks = seq(-1.5, 1.5, 0.5),
labels = paste(seq(-1.5, 1.5, 0.5)),
limits = c(-1.6, 1.6),
name = "") +
labs(x = "Easting [km]",
y = "Northing [km]",
title = "Gravity anomaly (mGal)") +
theme(panel.background = element_rect(colour = "black", fill = "white"),
panel.grid = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5, size=16),
aspect.ratio = 1)
The rectprismmag()
function calculates the magnetic
field anomaly from a series of 3-D rectangular prisms with vertical
sides using the algorithm of Plouff (1975b).
To calculate the magnetic anomaly of a single prism at a single point:
#########################################################
## magnetic anomaly of single prism at a single point ##
# location of the point where the magnetic anomaly will be calculated
<- data.frame(x=-0.3, y=-1.5, z=0)
magstation
# the rectangular prism is defined by its six edges (distance in km)
<- data.frame(xmin=-2, xmax=2,
prism1 ymin=-2, ymax=2,
zmin=-3, zmax=-1)
<- 0.015 # susceptiblity (SI)
susc <- 0 # remanent magnetization (A/m)
mstr <- 0 # remanent inclination (deg)
mincl <- 0 # remanent declination (deg)
mdecl <- 48800 # Earth's field intensity (nT)
ftotal <- 60 # field inclination (deg)
fincl <- 12 # field declination (deg)
fdecl
<- rectprismmag(magstation$x, magstation$y, magstation$z,
maganom $xmin, prism1$xmax,
prism1$ymin, prism1$ymax,
prism1$zmin, prism1$zmax, susc,
prism1
mstr, mincl, mdecl,
ftotal, fincl, fdecl)
# magnetic anomaly (nT)
print(maganom)
## [,1]
## [1,] 131.161
Calculate the magnetic anomaly at a 2-D grid of stations located above a 3-D grid of prisms. In this example, susceptibility values are generated randomly from a normal distribution with a mean of 0.015 SI units and a standard deviation of 0.001 SI units. No remanent magnetization is present in this example, so the magnetic anomaly is generated entirely from induced magnetization.
#########################################################
## magnetic anomaly from 3-D gridded source model ##
# first, create an array of stations
<- seq(-5, 5, by=0.1)
xstations <- seq(-5, 5, by=0.1)
ystations <- expand.grid(xstations, ystations)
D3.stations $Zstation <- 0
D3.stationsnames(D3.stations) <- c("Xstation", "Ystation", "Zstation")
# second, create an array of prisms
<- 0.1
width <- width / 2
half.width
<- seq(-2 + half.width, 2 - half.width, by=width)
xsource <- seq(-2 + half.width, 2 - half.width, by=width)
ysource <- seq(-3 + half.width, -1 - half.width, by=width)
zsource
<- expand.grid(xsource, ysource, zsource)
D3.source names(D3.source) <- c("xcenter", "ycenter", "zcenter")
# normally distributed, spatially uncorrelated, magnetic susceptibility
set.seed(1)
$suscnorm <- rnorm(n=length(D3.source[,1]), mean=0.015, sd=0.001)
D3.source
$nrmstr <- 0
D3.source$nrmdecl <- 0
D3.source$nrmincl <- 0
D3.source
# define the ambient field
$fieldtotal <- 48800
D3.source$fielddecl <- 12
D3.source$fieldincl <- 60
D3.source
<- rectprismmag(D3.stations$Xstation,
D3.maganomind $Ystation,
D3.stations$Zstation,
D3.stations$xcenter - half.width,
D3.source$xcenter + half.width,
D3.source$ycenter - half.width,
D3.source$ycenter + half.width,
D3.source$zcenter - half.width,
D3.source$zcenter + half.width,
D3.sourcesuscvolsi = D3.source$suscnorm,
nrmstr = D3.source$nrmstr,
nrmdecl = D3.source$nrmdecl,
nrmincl = D3.source$nrmincl,
fieldtotal = D3.source$fieldtotal,
fielddecl = D3.source$fielddecl,
fieldincl = D3.source$fieldincl,
bycell=FALSE)
# create data frame
<- D3.stations
d3.st $val <- D3.maganomind
d3.st
# plot with custom color scale
ggplot() +
geom_raster(data = as.data.frame(d3.st),
aes(x = Xstation, y = Ystation, fill = val)) +
scale_fill_gradientn(
colours = c("darkblue", "blue", "lightblue", "white", "pink", "red", "darkred"),
values = rescale(c(min(d3.st$val), 0, max(d3.st$val))),
breaks = c(-50, 0, 50, 100),
labels = paste(c("-50", 0, "50", "100")),
limits = c(-50, max(d3.st$val)),
name = "") +
labs(x = "Distance [km]",
y = "Distance [km]",
title = "Induced magnetic anomaly (nT)") +
theme(panel.background = element_rect(colour="black", fill="white"),
panel.grid = element_blank(),
plot.title = element_text(face = "bold", hjust = 0.5, size=16),
aspect.ratio = 1)
See the current DISCLAIMER