StatGeochem

juliaatomicnumber::NamedTuple A named tuple containing the atomic number of each element.

juliaioniccharge::NamedTuple A named tuple containing the ionic charges corresponding to the ionic radii in ionicradius.

juliaionicradius::NamedTuple A named tuple containing the ionic radii, in picometers, for the naturally ocurring elements in their geologically common redox states. Where multiple redox states may be expected for a single element in nature, they are suffixed to the atomic symbol, e.g. Fe2 vs Fe3.

Radii are those reported as "crystal" ionic radii in 6-fold coordination by Shannon [1], as tabulated at https://en.wikipedia.org/wiki/Ionic_radius#Tables

[1] R. D. Shannon (1976). "Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides". Acta Crystallogr A. 32 (5): 751-767. doi:10.1107/S0567739476001551

TC = Ayers_tsphene(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5)

Calculate sphene saturation temperature in degrees Celsius Following the sphene saturation calibration of Ayers et al., 2022 (doi: 10.1007/s00410-022-01902-z)

TiO2Sat = Ayers_tspheneTiO2(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5, T)

Calculate sphene saturation TiO2 concentration (in wt. %) for a given temperature (in C) following the sphene saturation calibration of Ayers et al., 2022 (doi: 10.1007/s00410-022-01902-z)

TiO2Sat = Ayers_tspheneTiO2_18(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5, T)

Calculate sphene saturation TiO2 concentration (in wt. %) for a given temperature (in C) following the sphene saturation calibration of Ayers et al., 2018 (doi: 10.1130/abs/2018AM-320568)

TC = Ayers_tsphene_18(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5)

Calculate sphene saturation temperature in degrees Celsius Following the sphene saturation calibration of Ayers et al., 2018 (doi: 10.1130/abs/2018AM-320568)

T = Boehnke_tzirc(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5, Zr)

Calculate zircon saturation temperature in degrees Celsius Following the zircon saturation calibration of Boehnke, Watson, et al., 2013 (doi: 10.1016/j.chemgeo.2013.05.028)

M = Boehnke_tzircM(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5)

Calculate zircon saturation M-value based on major element concentrations Following the zircon saturation calibration of Boehnke, Watson, et al., 2013 (doi: 10.1016/j.chemgeo.2013.05.028)

ZrSat = Boehnke_tzircZr(SiO2, TiO2, Al2O3, FeOT, MnO, MgO, CaO, Na2O, K2O, P2O5, T)

Calculate zircon saturation Zr concentration for a given temperature (in C) Following the zircon saturation calibration of Boehnke, Watson, et al., 2013 (doi: 10.1016/j.chemgeo.2013.05.028)

Ti = Crisp_Ti_in_zircon(TC::Number, Pbar::Number, aSiO2::Number, aTiO2::Number)

Parts per million by weight of titanium in zircon at temperature TC degrees Celsius and pressure Pbar bar given aSiO2 silica activity and aTiO2 titanium activity, following the equations of Crisp et al., 2023. (doi: 10.1016/j.gca.2023.04.031)

Ti = Ferry_Ti_in_zircon(TC::Number, aSiO2::Number, aTiO2::Number)

Parts per million by weight of titanium in zircon at temperature TC degrees Celsius given aSiO2 silica activity and aTiO2 titanium activity, following the equations of Ferry and Watson, 2007. (doi: 10.1007/s00410-007-0201-0)

TC = Ferry_Ti_in_zirconT(TC::Number, aSiO2::Number, aTiO2::Number)

Calculate titanium-in-zircon temperature in degrees Celcius TC given Ti parts per million by weight of titanium in zircon, aSiO2 silica activity and aTiO2 titanium activity, following the equations of Ferry and Watson, 2007. (doi: 10.1007/s00410-007-0201-0)

Zr = Ferry_Zr_in_rutile(TC::Number, aSiO2::Number)

Parts per million by weight of zirconium in rutile at temperature TC degrees Celsius given aSiO2 silica activity, following the equations of Ferry and Watson, 2007. (doi: 10.1007/s00410-007-0201-0)

TC = Ferry_Zr_in_rutileT(Zr::Number, aSiO2::Number)

Calculate zirconium-in-rutile temperature in degrees Celcius given Zr parts per million by weight of zirconium in rutile and aSiO2 silica activity, following the equations of Ferry and Watson, 2007. (doi: 10.1007/s00410-007-0201-0)

TC = Harrison_tapatite(SiO2, P2O5)

Calculate apatite saturation temperature in degrees Celcius following the apatite saturation model of Harrison and Watson 1984 (doi: 10.1016/0016-7037(84)90403-4)

As Harrison_tapatiteP2O5, but returns saturation phosphorus concentration in PPM P

P2O5 = Harrison_tapatiteP2O5(SiO2, Al2O3, CaO, Na2O, K2O, T)

Calculate P2O5 concentration (in wt.%) required for apatite saturation at a given T (in C) following the apatite saturation model of Harrison and Watson 1984 (doi: 10.1016/0016-7037(84)90403-4) with the correction of Bea et al. 1992 (doi: 10.1016/0024-4937(92)90033-U) where applicable

TC = Hayden_trutile(SiO2, TiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, P2O5)

Calculate rutile saturation temperature in degrees Celcius following the rutile saturation model of Hayden and Watson, 2007 (doi: 10.1016/j.epsl.2007.04.020)

TC = Hayden_trutile(SiO2, TiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, P2O5, TC)

Calculate the TiO2 concentration in weight percent required for rutile saturation at temperature TC degrees Celcius, following the rutile saturation model of Hayden and Watson, 2007 (doi: 10.1016/j.epsl.2007.04.020)

LREEmolwt(La, Ce, Pr, Nd, Sm, Gd)

Returns the average molecular weight of the LREE considered in the REEt value from the monazite saturation model of Montel 1993 (doi: 10.1016/0009-2541(93)90250-M)

LREEt(La, Ce, Pr, Nd, Sm, Gd)

Returns the sum of the LREE concentrations divided by their respective molar masses. If REE are input in parts per million by weight (ppmw), the result is in units of moles per megagram. This is equivalent to the REEt value from the monazite saturation model of Montel 1993 (doi: 10.1016/0009-2541(93)90250-M)

TC = Montel_tmonazite(SiO2, TiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, Li2O, H2O, La, Ce, Pr, Nd, Sm, Gd)

Calculate monazite saturation temperature in degrees Celcius following the monazite saturation model of Montel 1993 (doi: 10.1016/0009-2541(93)90250-M)

REEt = Montel_tmonaziteREE(SiO2, TiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, Li2O, H2O, T)

Calculate monazite saturation REEt value (in [ppm/mol.wt.]) for a given temperature (in C) following the monazite saturation model of Montel 1993 (doi: 10.1016/0009-2541(93)90250-M), where:

D = (Na + K + Li + 2Ca) / Al * 1/(Al + Si)) # all as molar cation fractions (not at. %!) ln(REEt) = 9.50 + 2.34D + 0.3879√H2O - 13318/T # H2O as wt.% REEt = Σ REEᵢ(ppm) / at. weight (g/mol)

LREEt = Rusiecka_tmonaziteREE(P_ppm, TC)

Calculate the LREEt (mol/Megagram) value required for monazite saturation at a temperature of TC degrees celcius and P ppmw phosphorous present, following the solubility model of Rusiecka & Baker, 2019 (doi: 10.2138/am-2019-6931)

LREEt = Rusiecka_txenotimeY(P_ppm, TC)

Calculate the Y (ppmw) concentration required for xenotime saturation at a temperature of TC degrees celcius and P ppmw phosphorous present, following the solubility model of Rusiecka & Baker, 2019 (doi: 10.2138/am-2019-6931)

TC = Tollari_tapatite(SiO2, TiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, P2O5)

Calculate apatite saturation temperature in degrees Celcius following the apatite saturation model of Tollari et al. 2006 (doi: 10.1016/j.gca.2005.11.024)

P2O5 = Tollari_tapatiteP2O5(SiO2, CaO, T)

Calculate P2O5 concentration (in wt.%) required for apatite saturation at a given T (in C) following the apatite saturation model of Tollari et al. 2006 (doi: 10.1016/j.gca.2005.11.024)

bin_bsr([f!::Function=nanbinmean!], x::Vector, y::VecOrMat, xmin:step:xmax, [w];
    	x_sigma = 0.05x,
    	y_sigma = 0.05y,
    	nresamplings = 1000,
    	sem = :sigma,
    	p = 0.2
)
bin_bsr([f!::Function=nanbinmean!], x::Vector, y::VecOrMat, xmin, xmax, nbins, [w];
    	x_sigma = 0.05x,
    	y_sigma = 0.05y,
    	nresamplings = 1000,
    	sem = :sigma,
    	p = 0.2
)

Returns the bincenters c, means or medians m, and uncertainties of the mean or median for a variable y binned by independent variable x into nbins equal bins between xmin and xmax, after nresamplings boostrap resamplings with acceptance probability p.

If a 2-d array (matrix) of y values is provided, each column will be treated as a separate variable, means and uncertainties will be returned column-wise.

Optional keyword arguments and defaults:

x_sigma = 0.05x

A vector representing the absolute uncertainty (standard deviation) of each x value

y_sigma = 0.05y

A vector representing the absolute uncertainty (standard deviation) of each y value

nresamplings = 1000

The number of resamplings to conduct

sem = :sigma

Format of the uncertainty estimate of the distribution of the mean. If :sigma is chosen, a tuple of three vectors (c, m, e) will be returned, where e is the standard error of the mean. If :CI or :pctile is chosen, a tuple of four vectors (c, m, el, eu) will be returned, where el and eu are the lower and upper bounds of the 95% confidence interval.

p = 0.2

Resampling probabilities, either as a scalar or a vector of the same length as x

Examples:

(c,m,e) = bin_bsr(nanbinmedian!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:sigma)

(c,m,el,eu) = bin_bsr(nanbinmean!, x, y, 0, 4000, 40, x_sigma=0.05x, p=probability, sem=:pctile)

(c, m, el, eu) = bin_bsr_ratio_medians(x::Vector, num::Vector, denom::Vector, xmin, xmax, nbins, [w];
    	x_sigma = 0.05x,
    	num_sigma = zeros(size(num)),
    	denom_sigma = zeros(size(denom)),
    	nresamplings = 1000,
    	p::Union{Number,Vector} = 0.2
)

Equivalent to bin_bsr_ratios(nanbinmedian!, ...)

bin_bsr_ratios([f!::Function=nanbinmean!], x::Vector, num::Vector, denom::Vector, xmin:step:xmax, [w];
    	x_sigma = 0.05x,
    	num_sigma = zeros(size(num)),
    	denom_sigma = zeros(size(denom)),
    	nresamplings = 1000,
    	p::Union{Number,Vector} = 0.2
)
bin_bsr_ratios([f!::Function=nanbinmean!], x::Vector, num::Vector, denom::Vector, xmin, xmax, nbins, [w];
    	x_sigma = 0.05x,
    	num_sigma = zeros(size(num)),
    	denom_sigma = zeros(size(denom)),
    	nresamplings = 1000,
    	p::Union{Number,Vector} = 0.2
)

Returns the bincenters c, means m, as well as upper (el) and lower (eu) 95% CIs of the mean for a ratio num/den binned by x into nbins equal bins between xmin and xmax, after nresamplings boostrap resamplings with acceptance probability p.

Examples

julia> (c, m, el, eu) = bin_bsr_ratios(nanbinmean!, x, num, denom, xmin:step:xmax; x_sigma=0.05x)

julia> (c, m, el, eu) = bin_bsr_ratios(nanbinmean!, x, num, denom, xmin, xmax, nbins; x_sigma=0.05x)

(bincenters, N) = bincounts(x::AbstractArray, xmin::Number, xmax::Number, nbins::Integer;
    	relbinwidth::Number=1
)

Tally the number of samples that fall into each of nbins equally spaced x bins between xmin and xmax, aligned with bin edges as xmin:(xmax-xmin)/nbins:xmax

A relbinwidth of 1 represents normal space-filling bins, while a larger value implies bin overlap.

See also histcounts for a more efficient implementation without variable bin width.

binmeans(x, y, xmin:step:xmax, [weight]; resamplingratio=1, relbinwidth=1)
binmeans(x, y, xmin, xmax, nbins, [weight]; resamplingratio=1, relbinwidth=1)

The means (ignoring NaNs) of y values binned by x, into each of nbins equally spaced x bins between xmin and xmax, returning bincenters, means, and standard errors of the mean.

A relbinwidth of 1 represents normal space-filling bins, while a larger value implies bin overlap.

To more efficiently calculate binned means without variable bin widths (or suncertainties), see instead nanbinmean/nanbinmean!.

Examples

(c,m,e) = binmeans(x, y, 0, 4000, 40)

binmedians(x, y, xmin:step:xmax; resamplingratio=1 relbinwidth=1)
binmedians(x, y, xmin, xmax, nbins; resamplingratio=1 relbinwidth=1)

The medians (ignoring NaNs) of y values binned by x, into each of nbins equally spaced x bins between xmin and xmax, returning bincenters, medians, and equivalent standard errors of the mean (1.4828 * median abolute deviation).

A relbinwidth of 1 represents normal space-filling bins, while a larger value implies bin overlap.

To more efficiently calculate binned medians without variable bin widths (or suncertainties), see instead nanbinmedian/nanbinmedian!.

Examples

(c,m,e) = binmedians(x, y, 0:100:4000)
(c,m,e) = binmedians(x, y, 0, 4000, 40)

bsr!([f::Function=gaussian], resampled::Array, index::Vector{Int}, data, sigma, p;
    	rng::AbstractRNG=MersenneTwister()
)

Fill resampled with data boostrap resampled from a (sample-per-row / element-per-column) dataset data with uncertainties sigma and resampling probabilities p, optionally using random numbers generated by f where f is a function of the form f(rng, data[i], sigma[i])

resampled = bsresample(data::AbstractArray, sigma, nrows, [p];
    	 kernel = gaussian,
    	 rng = MersenneTwister(),
    	 return_index = false
)

Bootstrap resample a (sample-per-row / element-per-column) array of data with uncertainties sigma and resampling probabilities p

resampled = bsresample(dataset::Union{Dict,NamedTuple}, nrows, [elements], [p];
    	 kernel = gaussian,
    	 rng = MersenneTwister()
)

Bootstrap resample a dictionary-based dataset with uncertainties stored either in dataset["err"] or dataset["[variable]_sigma"]

carbonateconversion!(dataset::NamedTuple)

Convert carbonates (CaCO3, MgCO3) into corresponding metal oxides and CO2 if extant, in place, as well as synchonizing TIC, TOC, TC, C and CO2. All are assumed to be reported in the same units, (likely wt. %) except for C, which is assumed to be equivalent to unitratio * TC,

x, y, z = cartesian(ρ, φ, θ)

Convert from spherical coordinates (ρ,φ,θ) to cartesian coordinates (x,y,z).

centroid(lats, lons)

Return the centroid of a set of latitudes and longitudes on a sphere.

changepoint(data, [sigma], nsteps; np, npmin, npmax)

Given an ordered array of data points, optionally with uncertainties sigma, use a Markov chain Monte Carlo approach based on that of Gallagher et al., 2010 (10.1016/j.epsl.2011.09.015) to estimate the position (by index) and optionally number of changepoints that best explain the data. Will return the results for nsteps steps of the Markov chain.

Optional keyword arguments:

np

Specify an exact integer number of changepoints.

npmin
nmpax

Specify the maximum and minimum possible integer number of changepoints. If np is not specified, the number of changepoints will allowed to vary freely between these bounds, as in the "transdimensional" approach of Gallagher et al.

Examples

julia> A = [randn(100).-2; randn(100).+2];

julia> dist = changepoint(A, 10000; np=1);

julia> dist[9000:end] # after burnin
1001-element Vector{Int64}:
 101
 101
 101
 101
   ⋮
 101
 101
 101

 julia> dist = changepoint(A, ones(size(A)), 10000; np=1)
10000×1 Matrix{Int64}:
  61
  61
  61
  61
   ⋮
 101
 101
 101

cipw_norm(SiO2, TiO2, Al2O3, Fe2O3, FeO, MnO, MgO, CaO, Na2O, K2O, P2O5)

Returns

quartz, orthoclase, plagioclase, corundum, nepheline, diopside, orthopyroxene, olivine, magnetite, ilmenite, apatite

claiborne_zircon_kd(elem::String, T::Number)

Returns the temperature-dependent zircon/melt partition coefficients of Claiborne et al. (2017) for the given element elem at temperature T degrees C, using equations of the form a * exp(b/TK) where TK = T + 273.15

In addition to the elements for which Claiborne provides equations, we also extrapolate these partition coefficients to include "La" and "Pr" using a lattce strain model, following the equations of Blundy and Wood (2003).

Reference: Claiborne, L. L., Miller, C. F., Gualda, G. A., Carley, T. L., Covey, A. K., Wooden, J. L., & Fleming, M. A. (2018). Zircon as magma monitor: Robust, temperature‐dependent partition coefficients from glass and zircon surface and rim measurements from natural systems. in Microstructural geochronology: Planetary records down to atom scale, 1-33. https://doi.org/10.1002/9781119227250.ch1

Examples

julia> claiborne_zircon_kd("Yb", 600)
1016.9198328977473

(c, m, el, eu) = constproportion(binbsrfunction::Function, dataset::Dict, xelem, [yelem | numelem, denomelem], xmin::Number, xmax::Number, nbins::Integer; 
    	norm_by = dataset["SiO2"], 
    	norm_bins = [43,55,65,78], 
    	x_sigma = dataset[xelem*"_sigma"], 
    	[y_sigma = dataset[yelem*"_sigma"] | num_sigma = dataset[numelem*"_sigma"], denom_sigma = dataset[denomelem*"_sigma"]],
    	nresamplings = 1000,
    	p = 0.2,
)

Call binbsrfunction repeatedly for each interval in norm_bins and return the results weighted by the fraction of norm_by that falls in each bin implied by norm_bins.

By default, combines the results assuming constant proportions of mafic, intermediate, and felsic lithologies, following the approach of Keller and Harrison 2020 (doi: 10.1073/pnas.2009431117)

Examples

julia> c,m,el,eu = constproportion(bin_bsr_means, dataset, "Age", "K2O", 0, 4000, 40);

julia> c,m,el,eu = constproportion(bin_bsr_ratio_medians, dataset, "Age", "La", "Yb", 0, 4000, 40);

digitize_plotline(img, line_color, xlims, ylims; atol=0.16)

Calculate approximate x (horizontal) and y (vertical) positions for a colored line in an image

Examples

img = load("xysin.png") # using FileIO, ImageIO
C = eltype(img)
(x,y) = digitize_plotline(img, C(0,0.604,0.976,1), [0,2pi], [-1.1,1.1])

digitize_plotmarkers(img, marker_color, xlims, ylims; atol=0.16)

Calculate approximate x (horizontal) and y (vertical) positions and position uncertainties for distinct colored markers in an image.

Examples

img = load("xyscatter.png") # using FileIO, ImageIO
C = eltype(img)
(x,dx,y,dy) = digitize_plotmarkers(img, C(0,0.604,0.976,1), [0,10], [0,10])

dist_uncert(lats, lons)

Find the decimal degree center and associated uncertainty (in arc degrees) from lists lats and lons of decimal degree coordinates.

Examples

(lat_ctr, lon_ctr, uncertainty) = dist_uncert(lats, lons)

eustar(Nd::Number, Sm::Number, Gd::Number, Tb::Number)

Calculate expected europium concentration, Eu*, based on abundance of adjacent rare earths.

Full four-element log-linear interpolation, assuming 3+ ionic radii and the chondritic abundances of Sun and McDonough 1989 (doi: 10.1144/gsl.sp.1989.042.01.19)

eustar(Sm::Number, Gd::Number)

Calculate expected europium concentration, Eu*, based on abundance of adjacent rare earths.

Simple geometric mean interpolation from Sm and Gd alone, assuming the chondritic abundances of Sun and McDonough 1989 (doi: 10.1144/gsl.sp.1989.042.01.19), that is Eu* = 0.0580*sqrt(Sm/0.1530 * Gd/0.2055)

feoconversion(FeO::Number=NaN, Fe2O3::Number=NaN, FeOT::Number=NaN, Fe2O3T::Number=NaN)

Compiles data from FeO, Fe2O3, FeOT, and Fe2O3T into a single FeOT value.

find_crust1_base(lat,lon,layer)

Return elevation (relative to sea level) of the layer base for a crust 1.0 layer at a given latitude and longitude.

Accepts lat and lon both as Numbers and as AbstractArrays, but given the overhead of opening and reading the crust1 files, you should generally aim to provide large arrays with as many values in a single query as possible.

Available layers:

    1 | :water
    2 | :ice
    3 | :upper_sediments
    4 | :middle_sediments
    5 | :lower_sediments
    6 | :upper_crust
    7 | :middle_crust
    8 | :lower_crust
    9 | :mantle

Results are returned in form (Vp, Vs, Rho, thickness)

Examples

julia> find_crust1_base([43.702245], [-72.0929], 8)
1-element Vector{Float64}:
 -36.26

julia> find_crust1_base([43.702245], [-72.0929], :lower_crust)
1-element Vector{Float64}:
-36.26

find_crust1_layer(lat,lon,layer)

Return all point data (Vp, Vs, Rho, layer thickness) for a given latitude, longitude, and crustal layer.

Accepts lat and lon both as Numbers and as AbstractArrays, but given the overhead of opening and reading the crust1 files, you should generally aim to provide large arrays with as many values in a single query as possible.

Available layers:

    1 | :water
    2 | :ice
    3 | :upper_sediments
    4 | :middle_sediments
    5 | :lower_sediments
    6 | :upper_crust
    7 | :middle_crust
    8 | :lower_crust
    9 | :mantle

Results are returned in form (Vp, Vs, Rho, thickness)

Examples

julia> vp, vs, rho, thickness = find_crust1_layer([43.702245], [-72.0929], 8)
([7.0], [3.99], [2950.0], [7.699999999999999])

find_crust1_seismic(lat,lon,layer)

Return all seismic data (Vp, Vs, Rho) for a given latitude, longitude, and crustal layer.

Accepts lat and lon both as Numbers and as AbstractArrays, but given the overhead of opening and reading the crust1 files, you should generally aim to provide large arrays with as many values in a single query as possible.

Available layers:

    1 | :water
    2 | :ice
    3 | :upper_sediments
    4 | :middle_sediments
    5 | :lower_sediments
    6 | :upper_crust
    7 | :middle_crust
    8 | :lower_crust
    9 | :mantle

Results are returned in form (Vp, Vs, Rho, thickness)

Examples

julia> vp, vs, rho = find_crust1_seismic([43.702245], [-72.0929], 8)
([7.0], [3.99], [2950.0])

julia> vp, vs, rho = find_crust1_seismic([43.702245], [-72.0929], :lower_crust)
([7.0], [3.99], [2950.0])

find_crust1_thickness(lat,lon,layer)

Return layer thickness for a crust 1.0 layer at a given latitude and longitude.

Accepts lat and lon both as Numbers and as AbstractArrays, but given the overhead of opening and reading the crust1 files, you should generally aim to provide large arrays with as many values in a single query as possible.

Available layers:

    1 | :water
    2 | :ice
    3 | :upper_sediments
    4 | :middle_sediments
    5 | :lower_sediments
    6 | :upper_crust
    7 | :middle_crust
    8 | :lower_crust
    9 | :mantle

Results are returned in form (Vp, Vs, Rho, thickness)

Examples

julia> find_crust1_thickness([43.702245], [-72.0929], 8)
1-element Vector{Float64}:
 7.699999999999999

julia> find_crust1_thickness([43.702245], [-72.0929], :lower_crust)
1-element Vector{Float64}:
7.699999999999999

find_etopoelev([etopo], lat, lon, [T=Float64])

Find the elevation of points at position (lat, lon) on the surface of the Earth, using the ETOPO1 one-arc-degree elevation model.

Units are meters of elevation and decimal degrees of latitude and longitude.

Reference: Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M. http://www.ngdc.noaa.gov/mgg/global/global.html

See also: get_etopo.

Examples

julia> etopo = get_etopo("elevation")
10800×21600 Matrix{Float64}:
   -58.0    -58.0    -58.0  …    -58.0    -58.0    -58.0
   -61.0    -61.0    -61.0       -61.0    -61.0    -61.0
   -62.0    -63.0    -63.0       -63.0    -63.0    -62.0
   -61.0    -62.0    -62.0       -62.0    -62.0    -61.0
     ⋮                      ⋱
 -4226.0  -4226.0  -4227.0     -4227.0  -4227.0  -4227.0
 -4228.0  -4228.0  -4229.0     -4229.0  -4229.0  -4229.0
 -4229.0  -4229.0  -4229.0     -4229.0  -4229.0  -4229.0

julia> find_etopoelev(etopo, 43.702245, -72.0929)
0-dimensional Array{Float64, 0}:
294.0

find_geolcont(lat,lon)

Find which geographic continent a sample originates from.

Continents:

  1: "Africa"
  2: "Eurasia"
  3: "North America"
  4: "South America"
  5: "Australia"
  6: "Antarctica"
  7: "NA"

See also: continents, continentcolors.

Examples

julia> find_geolcont(43.702245, -72.0929)
0-dimensional Array{Int64, 0}:
3

julia> continents[find_geolcont(43.702245, -72.0929)]
0-dimensional Array{String, 0}:
"North America"

find_geolprov(lat,lon)

Find which tectonic setting a sample originates from, based on a modified version of the USGS map of tectonic provinces of the world (c.f. https://commons.wikimedia.org/wiki/File:Worldgeologicprovinces.jpg)

Settings:

  10: Accreted Arc
  11: Island Arc
  12: Continental Arc
  13: Collisional orogen
  20: Extensional
  21: Rift
  22: Plume
  31: Shield
  32: Platform
  33: Basin
  00: No data

Settings returned are most representative modern setting at a given location and may not represent the tectonic setting where rocks (especially older/Precambrian rocks) originally formed.

Examples

julia> find_geolprov(43.702245, -72.0929)
0-dimensional Array{Int64, 0}:
10

julia> lat = rand(4)*180 .- 90
4-element Vector{Float64}:
 -28.352224011759773
  14.521710123066882
  43.301961981794335
  79.26368353708557

julia> lon = rand(4)*360 .- 180
4-element Vector{Float64}:
   5.024149409750521
 161.04362679392233
 123.21726489255786
 -54.34797401313695

julia> find_geolprov(lat, lon)
4-element Vector{Int64}:
  0
  0
 32
  0

find_land(lat,lon)

Find whether or not a given set of latitude, longitude points on the globe is above sea level, based on the etopo bedrock elevation dataset

Examples

julia> find_land(43.702245, -72.0929)
0-dimensional Array{Bool, 0}:
1

find_litho1_property(lat, lon, layer::Symbol, property::Symbol)

Return values for a LITHO1.0 property of a given layer at one or more given latitudes and longitudes, to the nearest 0.5-arc-degree grid point.

Accepts lat and lon both as Numbers and as AbstractArrays, but given the overhead of opening and reading the LITHO1.0 files, you should generally aim to provide large arrays with as many values in a single query as possible.

Available properties include:

    :vp             | p-wave velocity [m/s]
    :vs             | s-wave velocity [m/s]
    :rho            | density [kg/m^3]
    :bottom         | depth to bottom of the layer [km] (above sea level = negative)
    :thickness      | layer thickness [km]

while avialble layers are:

    1 | :ice
    2 | :water
    3 | :upper_sediments
    4 | :middle_sediments
    5 | :lower_sediments
    6 | :upper_crust
    7 | :middle_crust
    8 | :lower_crust
    9 | :sclm (or :lithosphere)
    10 | :asthenosphere

Examples

julia> find_litho1_property([43.702245, 44], [-72.0929, -73], :upper_crust, :vp)
2-element Vector{Float64}:
 6219.99
 6253.16

find_seafloorage([sfdata], lat, lon)

Find the age of the seafloor at positions (lat, lon) on the ocean crust, using the Müller et al. (2008) dataset.

Units are millions of years for age and mm/yr for rate.

Reference: M⁠⁠üller, R. D., M. Sdrolias, C. Gaina, and W. R. Roest (2008). Age, spreading rates, and spreading asymmetry of the world's ocean crust, Geochem. Geophys. Geosyst., 9, Q04006, doi:10.1029/2007GC001743. ftp://ftp.es.usyd.edu.au/pub/agegrid/2008/Grids/

See also: get_seafloorage.

Examples

julia> sfdata = get_seafloorage("seafloorage")
8640×10800 Matrix{Float64}:
   8.02    8.13    8.24    8.36  …    7.6     7.7     7.8     7.91
   8.01    8.12    8.23    8.34       7.61    7.71    7.81    7.91
   8.01    8.11    8.22    8.32       7.62    7.72    7.81    7.91
   8.01    8.11    8.2     8.3        7.64    7.73    7.82    7.91
   ⋮                             ⋱
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN

julia> find_seafloorage(sfdata, 43.702245, -40)
0-dimensional Array{Float64, 0}:
75.95

find_srtm15plus([srtm], lat, lon, [T=Float64])

Find the elevation of points at position (lat, lon) on the surface of the Earth, using the SRTM15plus 15-arc-second elevation model.

Units are meters of elevation and decimal degrees of latitude and longitude.

Reference: https://doi.org/10.5069/G92R3PT9

See also: get_srtm15plus.

Examples

julia> srtm = get_srtm15plus("elevation")
43201×86401 Matrix{Int16}:
 -32768  -32768  -32768  -32768  …  -32768  -32768  -32768
   3124    3124    3124    3124       3113    3113    3124
   3123    3123    3123    3122       3111    3111    3123
   3121    3121    3121    3121       3110    3110    3121
      ⋮                          ⋱                       ⋮
  -4225   -4224   -4224   -4224      -4224   -4225   -4225
  -4223   -4222   -4222   -4223      -4223   -4223   -4223
  -4223   -4223   -4223   -4223      -4223   -4223   -4223
  -4230   -4230   -4230   -4230  …   -4230   -4230   -4230

julia> find_srtm15plus(srtm, 43.702245, -72.0929)
0-dimensional Array{Float64, 0}:
252.0

find_tc1_age(lat::Number,lon::Number)

Return a tuple (age, age_min, age_max) containing the nominal, upper, and lower tc1 age bounds for the 1x1 arc degree grid cell containing lat and lon

find_tc1_age(lat::AbstractArray,lon::AbstractArray)

Return a tuple (age, age_min, age_max) where age, age_min, and age_max are arrays containing the nominal, upper and lower tc1 age bounds for each location pair lat[i], lon[i]

Data is sourced from the global 1x1° "TC1" thermal model of Artemieva 2006, doi: 10.1016/j.tecto.2005.11.022

find_tc1_crust(lat::Number,lon::Number)

Find the depth to the 550C isotherm for the 1x1 arc degree grid cell containing lat and lon

find_tc1_crust(lat::AbstractArray,lon::AbstractArray)

For each pair of latitudes and longitudes given by lat and lon, find the depth to the 550C isotherm for the 1x1 arc degree grid cell containing lat[i] and lon[i]

Data is sourced from the global 1x1° "TC1" thermal model of Artemieva 2006, doi: 10.1016/j.tecto.2005.11.022

find_tc1_lith(lat::Number,lon::Number)

Find the depth to the 1300C isotherm for the 1x1 arc degree grid cell containing lat and lon

find_tc1_lith(lat::AbstractArray,lon::AbstractArray)

For each pair of latitudes and longitudes given by lat and lon, find the depth to the 1300C isotherm for the 1x1 arc degree grid cell containing lat[i] and lon[i]

Data is sourced from the global 1x1° "TC1" thermal model of Artemieva 2006, doi: 10.1016/j.tecto.2005.11.022

get_etopo([varname])

Read ETOPO1 (1 arc minute topography) file from HDF5 storage, downloading from cloud if necessary.

Available varnames (variable names) include:

  "elevation"
  "y_lat_cntr"
  "x_lon_cntr"
  "cellsize"
  "scalefactor"
  "reference"

Units are meters of elevation and decimal degrees of latitude and longitude.

Reference: Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M. http://www.ngdc.noaa.gov/mgg/global/global.html

See also: find_etopoelev.

Examples

julia> get_etopo()
Dict{String, Any} with 6 entries:
  "cellsize"    => 0.0166667
  "scalefactor" => 60
  "x_lon_cntr"  => [-179.992, -179.975, -179.958, -179.942, -179.925, -1…
  "reference"   => "Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minut…
  "y_lat_cntr"  => [-89.9917, -89.975, -89.9583, -89.9417, -89.925, -89.…
  "elevation"   => [-58.0 -58.0 … -58.0 -58.0; -61.0 -61.0 … -61.0 -61.0…

julia> get_etopo("elevation")
10800×21600 Matrix{Float64}:
  -58.0    -58.0    -58.0  …    -58.0    -58.0    -58.0
  -61.0    -61.0    -61.0       -61.0    -61.0    -61.0
  -62.0    -63.0    -63.0       -63.0    -63.0    -62.0
  -61.0    -62.0    -62.0       -62.0    -62.0    -61.0
    ⋮                      ⋱
 -4226.0  -4226.0  -4227.0     -4227.0  -4227.0  -4227.0
 -4228.0  -4228.0  -4229.0     -4229.0  -4229.0  -4229.0
 -4229.0  -4229.0  -4229.0     -4229.0  -4229.0  -4229.0

get_seafloorage(varname="")

Read seafloor age file from HDF5 storage, downloading from cloud if necessary.

Available varnames (variable names) include:

  "seafloorage"
  "seafloorage_sigma",
  "seafloorrate"
  "information"
  "reference"

Units are millions of years for age and mm/yr for rate.

Reference: M⁠⁠üller, R. D., M. Sdrolias, C. Gaina, and W. R. Roest (2008). Age, spreading rates, and spreading asymmetry of the world's ocean crust, Geochem. Geophys. Geosyst., 9, Q04006, doi:10.1029/2007GC001743. ftp://ftp.es.usyd.edu.au/pub/agegrid/2008/Grids/

See also: find_seafloorage.

Examples

julia> get_seafloorage()
Dict{String, Any} with 5 entries:
  "seafloorrate"      => [6.83 6.79 … 6.91 6.87; 6.82 6.78 … 6.91 6.…
  "reference"         => "Muller, R. D., M. Sdrolias, C. Gaina, and …
  "seafloorage"       => [8.02 8.13 … 7.8 7.91; 8.01 8.12 … 7.81 7.9…
  "information"       => "Mercator projection, from 80.738 to -80.73…
  "seafloorage_sigma" => [2.15 2.19 … 2.08 2.11; 2.15 2.19 … 2.08 2.…

julia> get_seafloorage("seafloorage")
8640×10800 Matrix{Float64}:
   8.02    8.13    8.24    8.36  …    7.6     7.7     7.8     7.91
   8.01    8.12    8.23    8.34       7.61    7.71    7.81    7.91
   8.01    8.11    8.22    8.32       7.62    7.72    7.81    7.91
   8.01    8.11    8.2     8.3        7.64    7.73    7.82    7.91
   ⋮                             ⋱
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN
 NaN     NaN     NaN     NaN        NaN     NaN     NaN     NaN

get_srtm15plus([varname])

Read SRTM15plus file from HDF5 storage (15 arc second topography from the Shuttle Radar Topography Mission), downloading from cloud if necessary.

Available varnames (variable names) include:

  "elevation"
  "y_lat_cntr"
  "x_lon_cntr"
  "cellsize"
  "scalefactor"
  "nanval"
  "reference"

Units are meters of elevation and decimal degrees of latitude and longitude.

Reference: https://doi.org/10.5069/G92R3PT9

See also: find_srtm15plus.

Examples

julia> get_srtm15plus()
Dict{String, Any} with 7 entries:
  "cellsize"    => 0.00416667
  "scalefactor" => 240
  "x_lon_cntr"  => [-180.0, -179.996, -179.992, -179.988, -179.983,…
  "reference"   => "http://topex.ucsd.edu/WWW_html/srtm30_plus.html"
  "y_lat_cntr"  => [-90.0, -89.9958, -89.9917, -89.9875, -89.9833, …
  "nanval"      => -32768
  "elevation"   => Int16[-32768 -32768 … -32768 -32768; 3124 3124 ……

julia> get_srtm15plus("elevation")
43201×86401 Matrix{Int16}:
 -32768  -32768  -32768  -32768  …  -32768  -32768  -32768
   3124    3124    3124    3124       3113    3113    3124
   3123    3123    3123    3122       3111    3111    3123
   3121    3121    3121    3121       3110    3110    3121
      ⋮                          ⋱                       ⋮
  -4225   -4224   -4224   -4224      -4224   -4225   -4225
  -4223   -4222   -4222   -4223      -4223   -4223   -4223
  -4223   -4223   -4223   -4223      -4223   -4223   -4223
  -4230   -4230   -4230   -4230  …   -4230   -4230   -4230

haversine(lat₁, lon₁, lat₂, lon₂)

Calculate the arc degree distance between two decimal degree points (lat₁, lon₁) and (lat₂, lon₂).

k = invweight(lat::AbstractArray, lon::AbstractArray, age::AbstractArray;
    	lp::Number=2,
    	spatialscale=1.8,
    	agescale=38.0
)

Find the inverse weights k (proportional to spatiotemporal sample density) for a set of geological samples with specified latitude (lat), logitude (lon), and age (of crystallization, deposition, etc.).

The default spatialscale and agescale are taken from Keller and Schoene 2012. However, alternative scalings can be supplied. If an array is supplied for either spatialscale, agescale, or both, a 3-d matrix of k values will be returned, with dimensions length(spatialscale)length(agescale)nrows.

k = invweight(x::AbstractArray, xscale::Number; lp=2)

Find the inverse weights for a single array x for a given xscale, and exponent lp (default lp = 2).

Returns an array k where k[i] is the "inverse weight" for element i of the input array.

k = invweight_age(age::AbstractArray;
    	lp::Number=2, 
    	agescale::Number=38.0
)

Find the inverse weights k (proportional to temporal sample density) for a set of geological samples with specified age (of crystallization, deposition, etc.).

k = invweight_location(lat::AbstractArray, lon::AbstractArray;
    	lp::Number=2,
    	spatialscale::Number=1.8
)

Find the inverse weights k (proportional to spatial sample density) for a set of geological samples with specified latitude (lat), and logitude (lon).

materialize(x)

Convert an array-like object to an materialized, allocated Array, and leaving other types unchanged.

Unlike collect, will merely pass through an existing Array without change, rather than allocating new memory and making a copy.

Examples

julia> StatGeochem.materialize(1:100)
100-element Vector{Int64}:
   1
   2
   3
   ⋮
  99
 100

julia> StatGeochem.materialize(5)
5

melts_clean_modes(scratchdir::String; index=1)

Read and parse / clean-up modal phase proportions from specified MELTS run directory Returns an elementified dictionary

melts_configure(meltspath::String, scratchdir::String, composition::Collection{Number},
    	elements::Collection{String},
    	T_range=(1400, 600),
    	P_range=(10000,10000);)

Configure and run a MELTS simulation using alphaMELTS. Optional keyword arguments and defaults include:

batchstring::String = "1

sc.melts 10 1 3 1 liquid 1 1.0 0 10 0 4 0 "

A string defining the sequence of options that would be entered to produce the desired calculation if running alphaMELTS at the command line. The default string specifies a batch calculation starting at the liquidus.

dT = -10

The temperature step, in degrees, between each step of the MELTS calculation

dP = 0

The pressure step, in bar, between each step of the MELTS calculation

index = 1

An optional variable used to specify a unique suffix for the run directory name

version::String = "pMELTS"

A string specifying the desired version of MELTS. Options include MELTS and pMELTS.

mode::String = "isobaric"

A string specifying the desired calculation mode for MELTS. Options include isothermal, isobaric, isentropic, isenthalpic, isochoric, geothermal and PTPath.

fo2path::String = "FMQ"

A string specifying the oxygen fugacity buffer to follow, e.g., FMQ or NNO+1. Available buffers include IW,COH,FMQ,NNO,HM, and None

fractionatesolids::Bool = false

Fractionate all solids? default is false

suppress::Collection{String} = String[]

Supress individual phases (specify as strings in array, i.e. ["leucite"])

verbose::Bool = true

Print verbose MELTS output to terminal (else, write it to melts.log)

melts_query(scratchdir::String; index=1)

Read all phase proportions from Phase_main_tbl.txt in specified MELTS run directory Returns an elementified dictionary

melts_query_liquid(scratchdir::String; index=1)

Read liquid composition from Liquid_comp_tbl.txt in specified MELTS run directory Returns an elementified dictionary

melts_query_modes(scratchdir::String; index=1)

Read modal phase proportions from Phase_mass_tbl.txt in specified MELTS run Returns an elementified dictionary

melts_query_solid(scratchdir::String; index=1)

Read solid composition from Solid_comp_tbl.txt in specified MELTS run directory Returns an elementified dictionary

melts_query_system(scratchdir::String; index=1, importas=:Dict)

Read system thermodynamic and composition data from System_main_tbl.txt in specified MELTS run directory. Returns an elementified dictionary or tuple.

dataset = metalconversion!(dataset::Union{Dict,NamedTuple}; unitratio::Number=10000)

Convert minor element oxides (MnO, Cr2O3, NiO, ...) into corresponding metals (Mn, Cr, Ni, ...) in place if extant.

If metals are expected as parts per million (ppm), set unitratio=10000 (default); if metals are as wt%, set unitratio = 1

See also metalconversion, c.f. oxideconversion!

converted_dataset = metalconversion(dataset::Union{Dict,NamedTuple}; unitratio::Number=10000)

As metalconversion!, but returning a copy rather than modifying in-place

oxideconversion!(dataset::Dict; unitratio::Number=10000)

Convert major elements (Ti, Al, etc.) into corresponding oxides (TiO2, Al2O3, ...) in place if extant.

If metals are expected as PPM, set unitratio=10000 (default); if metals are as wt%, set unitratio = 1

See also oxideconversion, c.f. metalconversion!

converted_dataset = oxideconversion(dataset::Union{Dict,NamedTuple}; unitratio::Number=10000)

As oxideconversion!, but returning a copy rather than modifying in-place

perplex_configure_geotherm(scratchdir::String, composition, [elements],
    	P_range=(280,28000), T_surf::Number=273.15, geotherm::Number=0.1;
    	dataset::String="hp633ver.dat",
    	index::Integer=1,
    	npoints::Integer=100,
    	solution_phases::String="melt(HGPH)\nPl(I1,HP)\nFsp(C1)\nSp(HGP)\nGt(HGP)\nO(HGP)\nOpx(HGP)\nCpx(HGP)\nCrd(HGP)\nBi(HGP)\nMica(W)\nEp(HP)\ncAmph(G)\nIlm(WPH)\nChl(W)\n",
    	excludes::String="ged\nfanth\ng\n",
    	mode_basis::String="wt",  #["vol", "wt", "mol"]
    	composition_basis::String="wt",  #["vol", "wt", "mol"]
    	fluid_eos::Integer=5)

Set up a PerpleX calculation for a single bulk composition along a specified geothermal gradient and pressure (depth) range. P specified in bar and T_surf in Kelvin, with geothermal gradient in units of Kelvin/bar

perplex_configure_isobar(scratchdir::String, composition, [elements],
    	P::Number=10000, T_range::NTuple{2,Number}=(500+273.15, 1500+273.15);
    	dataset::String="hp633ver.dat",
    	index::Integer=1,
    	npoints::Integer=100,
    	solution_phases::String="melt(HGPH)\nPl(I1,HP)\nFsp(C1)\nSp(HGP)\nGt(HGP)\nO(HGP)\nOpx(HGP)\nCpx(HGP)\nCrd(HGP)\nBi(HGP)\nMica(W)\nEp(HP)\ncAmph(G)\nIlm(WPH)\nChl(W)\n",
    	excludes::String="ged\nfanth\ng\n",
    	mode_basis::String="wt",  #["vol", "wt", "mol"]
    	composition_basis::String="wt",  #["wt", "mol"]
    	nonlinear_subdivision::Bool=false,
    	fluid_eos::Integer=5)

Set up a PerpleX calculation for a single bulk composition along a specified isobaric temperature gradient. P specified in bar and T_range in Kelvin

perplex_configure_path(scratchdir::String, composition::Collection{Number}, PTdir::String="", PTfilename::String="",
    	elements::String=("SiO2","TiO2","Al2O3","FeO","MgO","CaO","Na2O","K2O","H2O"),
    	T_range::NTuple{2,Number}=(500+273.15, 1050+273.15);
    	dataset::String="hp633ver.dat",
    	index::Integer=1,
    	solution_phases::String="melt(HGPH)\nPl(I1,HP)\nFsp(C1)\nSp(HGP)\nGt(HGP)\nO(HGP)\nOpx(HGP)\nCpx(HGP)\nCrd(HGP)\nBi(HGP)\nMica(W)\nEp(HP)\ncAmph(G)\nIlm(WPH)\nChl(W)\n",
    	excludes::String="ged\nfanth\ng\n",
    	mode_basis::String="wt",  #["vol", "wt", "mol"]
    	composition_basis::String="wt",  #["wt", "mol"]
    	nonlinear_subdivision::Bool=false,
    	fluid_eos::Integer=5,
    	fractionate::Integer=0)

Set up a PerpleX calculation for a single bulk composition along a specified pressure–temperature path with T as the independent variable.

P specified in bar and T_range in Kelvin

perplex_configure_pseudosection(scratchdir::String, composition, [elements::Collection{String}],
    	P::NTuple{2,Number}=(280, 28000), T::NTuple{2,Number}=(273.15, 1500+273.15);
    	dataset::String="hp633ver.dat",
    	index::Integer=1,
    	xnodes::Integer=42,
    	ynodes::Integer=42,
    	solution_phases::String="melt(HGPH)\nPl(I1,HP)\nFsp(C1)\nSp(HGP)\nGt(HGP)\nO(HGP)\nOpx(HGP)\nCpx(HGP)\nCrd(HGP)\nBi(HGP)\nMica(W)\nEp(HP)\ncAmph(G)\nIlm(WPH)\nChl(W)\n",
    	excludes::String="ged\nfanth\ng\n",
    	mode_basis::String="wt", #["vol", "wt", "mol"]
    	composition_basis::String="wt", #["wt", "mol"]
    	fluid_eos::Number=5)

Set up a PerpleX calculation for a single bulk composition across an entire 2d P-T space. P specified in bar and T in Kelvin.

perplex_query_modes(scratchdir::String, P::NTuple{2,Number}, T::NTuple{2,Number};
    	index::Integer=1, 
    	npoints::Integer=200, 
    	include_fluid="y", 
    	importas=:Dict,
)

perplex_query_modes(scratchdir::String, P::AbstractArray, T::AbstractArray;
    	index::Integer=1, 
    	include_fluid="y", 
    	importas=:Dict,
)

Query modal mineralogy (mass proportions) along a specified P-T path using a pre-computed pseudosection. Results are returned as a dictionary.

perplex_query_modes(scratchdir::String;
    	index::Integer=1, 
    	dof::Integer=1, 
    	include_fluid="y", 
    	npoints::Integer=0,
    	manual_grid::Bool=npoints>0,
    	importas=:Dict,
)

Query modal mineralogy (mass proportions) along a previously configured 1-d path (dof=1, isobar, geotherm, or P–T path) or 2-d grid / pseudosection (dof=2). Results are returned as a dictionary.

Currently returns wt%

perplex_query_phase(scratchdir::String, phase::String, P::NTuple{2,Number}, T::NTuple{2,Number};
    	index::Integer=1, 
    	npoints::Integer=200, 
    	include_fluid="y", 
    	clean_units::Bool=true, 
    	importas=:Dict,
)

perplex_query_phase(scratchdir::String, phase::String, P::AbstractArray, T::AbstractArray;
    	index::Integer=1, 
    	include_fluid="y", 
    	clean_units::Bool=true, 
    	importas=:Dict,
)

Query all perplex-calculated properties for a specified phase (e.g. "Melt(G)") along a specified P-T path using a pre-computed pseudosection. Results are returned as a dictionary.

perplex_query_phase(scratchdir::String, phase::String;
    	index::Integer=1, 
    	dof::Integer=1,
    	include_fluid="y", 
    	clean_units::Bool=true, 
    	npoints::Integer=0,
    	manual_grid::Bool=npoints>0, 
    	importas=:Dict,
)

Query all perplex-calculated properties for a specified phase (e.g. "Melt(G)") along a previously configured 1-d path (dof=1, isobar, geotherm, or P–T path) or 2-d grid / pseudosection (dof=2). Results are returned as a dictionary.

perplex_query_point(scratchdir::String, P::Number, T::Number; index::Integer=1)

Query perplex results at a single P,T point in a pseudosection. Results are returned as a string.

perplex_query_point(scratchdir::String, indvar::Number; index::Integer=1)

Query perplex results at a single temperature on an isobar or single pressure on a geotherm. Results are returned as a string.

perplex_query_seismic(scratchdir::String, P::NTuple{2,Number}, T::NTuple{2,Number};
    	index::Integer=1,
    	npoints::Integer=200, 
    	include_fluid="n", 
    	importas=:Dict,
)

perplex_query_seismic(scratchdir::String, P::AbstractArray, T::AbstractArray;
    	index::Integer=1, 
    	include_fluid="n", 
    	importas=:Dict,
)

Query perplex seismic results along a specified P-T path using a pre-computed pseudosection. Results are returned as a dictionary.

perplex_query_seismic(scratchdir::String;
    	index::Integer=1, 
    	dof::Integer=1, 
    	include_fluid::String="n", 
    	npoints::Integer=0
    	manual_grid::Bool=npoints>0, 
    	importas=:Dict,
)

Query perplex seismic results along a previously configured 1-d path (dof=1, isobar or geotherm) or 2-d grid / pseudosection (dof=2). Results are returned as a dictionary.

perplex_query_system(scratchdir::String, P::NTuple{2,Number}, T::NTuple{2,Number};
    	index::Integer=1, 
    	npoints::Integer=200, 
    	include_fluid="y", 
    	clean_units::Bool=true, 
    	importas=:Dict,
)

perplex_query_system(scratchdir::String, P::AbstractArray, T::AbstractArray;
    	index::Integer=1, 
    	include_fluid="y", 
    	clean_units::Bool=true, 
    	importas=:Dict,
)

Query all perplex-calculated properties for the system (with or without fluid) along a specified P-T path using a pre-computed pseudosection. Results are returned as a dictionary. Set include_fluid="n" to return solid+melt only.

perplex_query_system(scratchdir::String;
    	index::Integer=1, 
    	dof::Integer=1, 
    	include_fluid="y", 
    	clean_units::Bool=true,
    	npoints::Integer=0,
    	manual_grid::Bool=npoints>0, 
    	importas=:Dict,
)

Query all perplex-calculated properties for the system (with or without fluid) along a previously configured 1-d path (dof=1, isobar or geotherm) or 2-d grid / pseudosection (dof=2). Results are returned as a dictionary. Set include_fluid="n" to return solid+melt only.

perplextrace_query(scratchdir, composition::LinearTraceComposition, args...; 
    	apatite = :Harrison,                  # Harrison & Watson 1984 + Bea et al. 1992 (doi: 10.1016/0024-4937(92)90033-U)
    	zircon = :Boehnke,                    # Boehnke, Watson, et al., 2013 (doi: 10.1016/j.chemgeo.2013.05.028)
    	sphene = :Ayers,                      # Ayers et al., 2018 (doi: 10.1130/abs/2018AM-320568)
    	monazite = :Montel,                   # Montel 1993 (doi: 10.1016/0009-2541(93)90250-M)
    	zircon_kd_source = :average,          # [:Claiborne | :GERM | :average]
    	export_bulk_kds::Bool = false,
    	export_mineral_kds::Bool = false,
    	export_mineral_compositions::Bool = false,
    	export_melt_parameters::Bool = true,
    	export_empty_columns::Bool = false,
    	require_phase_for_export = "",
    	exportas = :Dict,
    	kwargs...,                            # Any other keyword arguments are forwarded to `perplex_query_` functions
)

Calculate trace element compositions of all phases, including accessory minerals, by

1. Calculating the phase assemblage and (if applicable) melt composition given the bulk 
    major element composition via equilibrium thermodynamics using Perple_X

2. Calculating equilibrium trace element concentrations in all phases (particularly
    including melt, if present) using GERM-based partition coefficients (`germ_kd`) 

3. Evaluating the saturation state of accessory minerals given the calculated 
    trace element abundances
    
4. Recalculating trace element partitioning in the presence of newly-saturated
    accessory phases

The results are returned as a dictionary.

perplextrace_query_melt(scratchdir, composition::LinearTraceComposition, args...; 
    	apatite = :Harrison,                  # Harrison & Watson 1984 + Bea et al. 1992 (doi: 10.1016/0024-4937(92)90033-U)
    	zircon = :Boehnke,                    # Boehnke, Watson, et al., 2013 (doi: 10.1016/j.chemgeo.2013.05.028)
    	sphene = :Ayers,                      # Ayers et al., 2018 (doi: 10.1130/abs/2018AM-320568)
    	monazite = :Montel,                   # Montel 1993 (doi: 10.1016/0009-2541(93)90250-M)
    	zircon_kd_source = :average,          # [:Claiborne | :GERM | :average]
    	export_melt_parameters::Bool = true,
    	export_empty_columns::Bool = false,
    	require_phase_for_export = "",
    	exportas = :Dict,
    	kwargs...,                            # Any other keyword arguments are forwarded to `perplex_query_` functions
)

As perplextrace_query, but only return melt properties

The results are returned as a dictionary.

randsample(data, nrows, [p])

Bootstrap resample (without uncertainty) a data array to length nrows. Optionally provide weights p either as a vector (one-weight-per-sample) or scalar.

randsample(dataset::Dict, nrows, [elements], [p])

Bootstrap resample (without uncertainty) a dataset dict to length nrows. Optionally provide weights p either as a vector (one-weight-per-sample) or scalar.

resamplingprobability(k; median_probability::Number=1/6)

Calculate scaled resampling probabilities p given a vector of inverse weights k and a desired median resampling probability (1/6 by default).

ρ, θ, φ = cartesian(x, y, z)

Convert from cartesian coordinates (x,y,z) to spherical coordinates (ρ,φ,θ).

system(cmdstr::AbstractString)

Direct access to the command line through C's system function – without stripping/sanitizing special characters, in contrast to Julia's safer run() function. This allows pipelining, etc. in shell commands. Returns 0 on success.

Examples

julia> system("pwd")
/Users/foo/code/StatGeochem.jl
0