![\"Launch](\"https://mybinder.org/badge_logo.svg\"){kind=icon}
\n",
"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!
\n", "\n", "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.\n", "***\n", "### Load required Julia packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": false }, "outputs": [], "source": [ "## --- Load (and install if neccesary) the StatGeochem package which has the resampling functions we'll want\n", "try\n", " using StatGeochem\n", "catch\n", " using Pkg\n", " Pkg.add(PackageSpec(url=\"https://github.com/brenhinkeller/StatGeochem.jl\"))\n", " using StatGeochem\n", "end\n", "\n", "using Statistics, StatsBase, DelimitedFiles\n", "using Plots; gr();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load dataset" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "Dict{String,Any} with 117 entries:\n", " \"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.…\n", " \"Pd\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 0.01, N…\n", " \"MgO\" => [0.89, 0.5, 0.69, 2.83, 3.5, 2.81, 2.97, 2.32, 2.33, 1.74 … …\n", " \"C\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, Na…\n", " \"Nb\" => [3.5, 1.3, 3.4, 11.6, 10.4, 12.2, 11.5, 8.3, 7.0, 4.9 … 11.…\n", " \"Ag\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 0.66, N…\n", " \"Gd\" => [NaN, NaN, NaN, 7.99, 5.92, 6.12, 5.94, 5.46, 2.01, 2.7 … N…\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.…\n", " \"Age\" => [1.0, 1.0, 1.0, 1650.0, 1650.0, 1650.0, 1650.0, 1650.0, 1650.…\n", " \"Sb\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, 1.…\n", " \"TiO2\" => [0.27, 0.25, 0.24, 0.6, 0.63, 0.63, 0.59, 0.39, 0.48, 0.38 ……\n", " \"Cs\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, 7.…\n", " \"Be\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 2.7, Na…\n", " \"Locality\" => [\"Adachi, Japan tonalite \", \"Adachi, Japan tonalite \", \"Adach…\n", " \"Sr\" => [312.0, 396.0, 272.0, 887.4, 798.0, 761.4, 698.7, 602.0, 790.…\n", " \"Material\" => [\"IGNEOUS\", \"IGNEOUS\", \"IGNEOUS\", \"IGNEOUS\", \"IGNEOUS\", \"IGNE…\n", " \"Ga\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 29.0, N…\n", " \"Ta\" => [NaN, NaN, NaN, 0.62, 0.59, 0.98, 0.67, 0.38, 0.41, 0.21 … …\n", " \"Te\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … NaN, Na…\n", " \"Eustar\" => [NaN, NaN, NaN, 3.18574, 2.14821, 2.45012, 2.3435, 2.07734, 0…\n", " \"Al2O3\" => [19.19, 14.9, 15.77, 16.26, 17.65, 16.87, 16.21, 15.67, 15.38…\n", " \"CO2\" => [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN … 7.3, 0.…\n", " \"Freeair\" => [31.0, 31.0, 31.0, 140.0, 140.0, 140.0, 140.0, 140.0, 140.0, …\n", " \"Y\" => [22.4, 7.5, 17.6, 29.1, 21.9, 22.0, 23.4, 19.3, 6.3, 8.2 … …\n", " \"U\" => [NaN, NaN, NaN, 1.18, 1.24, 3.33, 1.25, 1.09, 0.84, 1.42 … …\n", " ⋮ => ⋮" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## --- Download and unzip Keller and Schoene (2012) dataset\n", "\n", "if ~isfile(\"ign.h5\") # Unless it already exists\n", " download(\"https://storage.googleapis.com/statgeochem/ign.h5.gz\",\"./ign.h5.gz\")\n", " download(\"https://storage.googleapis.com/statgeochem/err2srel.csv\",\"./err2srel.csv\")\n", " run(`gunzip -f ign.h5.gz`) # Unzip file\n", "end\n", "\n", "# Read HDF5 file\n", "using HDF5\n", "ign = h5read(\"ign.h5\",\"vars\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "## --- Compute proximity coefficients (inverse weights)\n", "\n", "# # Compute inverse weights\n", "# k = invweight(ign[\"Latitude\"] .|> Float32, ign[\"Longitude\"] .|> Float32, ign[\"Age\"] .|> Float32)\n", "\n", "# Since this is somewhat computatually intensive, let's load a precomputed version instead\n", "k = ign[\"k\"]\n", "\n", "# Probability of keeping a given data point when sampling\n", "p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the data in each resampling\n", "p[vec(ign[\"Elevation\"].<-100)] .= 0 # Consider only continental crust\n", "\n", "# Set absolute uncertainties for each variable where possible, using errors defined inerr2srel.csv\n", "err2srel = importdataset(\"err2srel.csv\", ',')\n", "for e in ign[\"elements\"]\n", " # If there's an err2srel for this variable, create a \"_sigma\" if possible\n", " if haskey(err2srel, e) && !haskey(ign, e*\"_sigma\")\n", " ign[e*\"_sigma\"] = ign[e] .* (err2srel[e] / 2);\n", " end\n", "end\n", "\n", "# Special cases: age uncertainty\n", "ign[\"Age_sigma\"] = (ign[\"Age_Max\"]-ign[\"Age_Min\"])/2;\n", "t = (ign[\"Age_sigma\"] .< 50) .| isnan.(ign[\"Age_sigma\"]) # Find points with < 50 Ma absolute uncertainty\n", "ign[\"Age_sigma\"][t] .= 50 # Set 50 Ma minimum age uncertainty (1-sigma)\n", "\n", "# Special cases: location uncertainty\n", "ign[\"Latitude_sigma\"] = ign[\"Loc_Prec\"]\n", "ign[\"Longitude_sigma\"] = ign[\"Loc_Prec\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Introduction: major and trace-element variability at constant silica\n", "To demonstrate the motivation for this reference model, we first consider the variability in major and trace element concentrations _at constant silica_ over time," ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Single element differentiation example\n", "xelem = \"SiO2\"\n", "xmin = 45\n", "xmax = 75\n", "nbins = 8\n", "elem = \"MgO\"\n", "\n", "h = plot(xlabel=xelem, ylabel=\"$(elem)\",xlims=(xmin,xmax),framestyle=:box,grid=:off,fg_color_legend=:white) # Format plot\n", "\n", "rt = [0,1,2,3,4]\n", "colors = reverse(resize_colormap(viridis[1:end-20],length(rt)-1))\n", "for i=1:length(rt)-1\n", " t = (ign[\"Age\"].>rt[i]*1000) .& (ign[\"Age\"].