This Jupyter notebook uses the StatGeochem package to calculate and plot the geochemical evolution of igneous whole-rock geochemical compositions preserved in the extant continental crust, for any element or ratio included in the dataset of Keller & Schoene 2012 and 2018, using the weighted boostrap resampling algorithm therefrom, implemented in the Julia language.
If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!
Hint: shift
-enter
to run a single cell, or from the Cell
menu select Run All
to run the whole file. Any code from this notebook can be copied and pasted into the Julia REPL or a .jl
script.
## --- Load the StatGeochem package which has the resampling functions we'll want
using StatGeochem
using Plots
## --- Download and unzip the Keller and Schoene (2012) dataset
if ~isfile("ign.h5") # Unless it already exists
download("https://storage.googleapis.com/statgeochem/ign.h5.gz","./ign.h5.gz")
download("https://storage.googleapis.com/statgeochem/err2srel.csv","./err2srel.csv")
run(`gunzip -f ign.h5.gz`) # Unzip file
end
# Read HDF5 file
using HDF5
ign = h5read("ign.h5","vars")
Dict{String, Any} with 117 entries: "Lower_Vp" => [7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2, 7.2 … 7.2, 7.… "Pd" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 0.01, N… "MgO" => [0.89, 0.5, 0.69, 2.83, 3.5, 2.81, 2.97, 2.32, 2.33, 1.74 … … "C" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, Na… "Nb" => [3.5, 1.3, 3.4, 11.6, 10.4, 12.2, 11.5, 8.3, 7.0, 4.9 … 11.… "Ag" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 0.66, N… "Gd" => [NaN, NaN, NaN, 7.99, 5.92, 6.12, 5.94, 5.46, 2.01, 2.7 … N… "Middle_Rho" => [2.9, 2.9, 2.9, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8, 2.8 … 2.9, 2.… "Age" => [1.0, 1.0, 1.0, 1650.0, 1650.0, 1650.0, 1650.0, 1650.0, 1650.… "Sb" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, 1.… "TiO2" => [0.27, 0.25, 0.24, 0.6, 0.63, 0.63, 0.59, 0.39, 0.48, 0.38 …… "Cs" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, 7.… "Be" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 2.7, Na… "Locality" => ["Adachi, Japan tonalite ", "Adachi, Japan tonalite ", "Adach… "Sr" => [312.0, 396.0, 272.0, 887.4, 798.0, 761.4, 698.7, 602.0, 790.… "Material" => ["IGNEOUS", "IGNEOUS", "IGNEOUS", "IGNEOUS", "IGNEOUS", "IGNE… "Ga" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 29.0, N… "Ta" => [NaN, NaN, NaN, 0.62, 0.59, 0.98, 0.67, 0.38, 0.41, 0.21 … … "Te" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, Na… "Eustar" => [NaN, NaN, NaN, 3.18574, 2.14821, 2.45012, 2.3435, 2.07734, 0… "Al2O3" => [19.19, 14.9, 15.77, 16.26, 17.65, 16.87, 16.21, 15.67, 15.38… "CO2" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 7.3, 0.… "Freeair" => [31.0, 31.0, 31.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, … "Y" => [22.4, 7.5, 17.6, 29.1, 21.9, 22.0, 23.4, 19.3, 6.3, 8.2 … … "U" => [NaN, NaN, NaN, 1.18, 1.24, 3.33, 1.25, 1.09, 0.84, 1.42 … … ⋮ => ⋮
## --- Compute proximity coefficients (inverse weights)
# # Calculate inverse weights based on spatiotemporal sample density
# k = invweight(ign["Latitude"] .|> Float32, ign["Longitude"] .|> Float32, ign["Age"] .|> Float32)
# Since this is pretty computatually intensive, let's load a precomputed version instead
k = ign["k"]
# Probability of keeping a given data point when sampling
p = 1.0./((k.*nanmedian(5.0./k)) .+ 1.0); # Keep roughly one-fith of the data in each resampling
p[vec(ign["Elevation"].<-100)] .= 0 # Consider only continental crust
# Set absolute uncertainties for each element where possible, using errors defined inerr2srel.csv
err2srel = importdataset("err2srel.csv", ',', importas=:Dict)
for e in ign["elements"]
# If there's an err2srel for this variable, create a "_sigma" if possible
if haskey(err2srel, e) && !haskey(ign, e*"_sigma")
ign[e*"_sigma"] = ign[e] .* (err2srel[e] / 2);
end
end
# Special cases: age uncertainty
ign["Age_sigma"] = (ign["Age_Max"]-ign["Age_Min"])/2;
t = (ign["Age_sigma"] .< 50) .| isnan.(ign["Age_sigma"]) # Find points with < 50 Ma absolute uncertainty
ign["Age_sigma"][t] .= 50 # Set 50 Ma minimum age uncertainty (1-sigma)
# Special cases: location uncertainty
ign["Latitude_sigma"] = ign["Loc_Prec"]
ign["Longitude_sigma"] = ign["Loc_Prec"]
68696-element Vector{Float64}: 0.5 0.5 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ⋮ 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
## --- Resample a single variable
tmin = 50 # Minimum age
tmax = 3850 # Maximum age
nbins = 38
elem = "MgO"
# Look only at samples in the basaltic silica range
t = 43 .< ign["SiO2"] .< 51 # Mafic
# Resample, returning binned means and uncertainties
# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
(c,m,el,eu) = bin_bsr_means(ign["Age"][t], ign[elem][t], tmin,tmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t])
# Plot results
plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label="")
plot!(xlabel="Age (Ma)", ylabel="$elem (wt. %)",framestyle=:box,grid=:off,xflip=true) # Format plot
## --- Resample a ratio
tmin = 50 # Minimum age
tmax = 3850 # Maximum age
nbins = 38
num = "La" # Numerator
denom = "Yb" # Denominator
# Look only at samples in the basaltic silica range
t = 43 .< ign["SiO2"] .< 51 # Mafic
# Exclude outliers
t = t .& inpctile(ign[num], 99) .& inpctile(ign[denom], 99)
# Resample, returning binned means and uncertainties
# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
(c,m,el,eu) = bin_bsr_ratios(ign["Age"][t], ign[num][t], ign[denom][t], tmin,tmax,nbins, p=p[t], x_sigma=ign["Age_sigma"][t], num_sigma=ign[num][t]*0.05, denom_sigma=ign[denom][t]*0.05)
# Plot results
h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkred,markerstrokecolor=:darkred,label="")
plot!(h, xlabel="Age (Ma)", ylabel="$(num) / $(denom)",framestyle=:box,grid=:off,xflip=true) # Format plot
display(h)
## --- Ratio differentiation
xelem = "SiO2"
xmin = 40 # Minimum age
xmax = 80 # Maximum age
nbins = 20
num = "Sc" # Numerator
denom = "Yb" # Denominator
# Exclude outliers
t = inpctile(ign[num], 99) .& inpctile(ign[denom], 99)
# Resample, returning binned means and uncertainties
# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
(c,m,el,eu) = bin_bsr_ratios(ign[xelem][t], ign[num][t], ign[denom][t], xmin,xmax,nbins, p=p[t], x_sigma=ign[xelem][t]*0.01, num_sigma=ign[num][t]*0.05, denom_sigma=ign[denom][t]*0.05)
# Plot results
h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label="")
plot!(h, xlabel=xelem, ylabel="$(num) / $(denom)",xlims=(xmin,xmax),framestyle=:box,grid=:off) # Format plot
display(h)