![\"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": {}, "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": "code", "execution_count": 2, "metadata": {}, "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 the Keller and Schoene (2012) dataset\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": [ { "data": { "text/plain": [ "68696-element Array{Float64,1}:\n", " 0.5\n", " 0.5\n", " 0.5\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " 0.2\n", " ⋮\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01\n", " 0.01" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## --- Compute proximity coefficients (inverse weights)\n", "\n", "# # Calculate inverse weights based on spatiotemporal sample density\n", "# k = invweight(ign[\"Latitude\"] .|> Float32, ign[\"Longitude\"] .|> Float32, ign[\"Age\"] .|> Float32)\n", "\n", "# Since this is pretty 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 element 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": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "## --- Resample a single variable\n", "tmin = 50 # Minimum age\n", "tmax = 3850 # Maximum age\n", "nbins = 38\n", "elem = \"MgO\"\n", "\n", "# Look only at samples in the basaltic silica range\n", "t = 43 .< ign[\"SiO2\"] .< 51 # Mafic\n", "\n", "# Resample, returning binned means and uncertainties\n", "# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)\n", "(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])\n", "\n", "# Plot results\n", "plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label=\"\")\n", "plot!(xlabel=\"Age (Ma)\", ylabel=\"$elem (wt. %)\",framestyle=:box,grid=:off,xflip=true) # Format plot" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Resample a ratio\n", "tmin = 50 # Minimum age\n", "tmax = 3850 # Maximum age\n", "nbins = 38\n", "num = \"La\" # Numerator\n", "denom = \"Yb\" # Denominator\n", "\n", "# Look only at samples in the basaltic silica range\n", "t = 43 .< ign[\"SiO2\"] .< 51 # Mafic\n", "\n", "# Exclude outliers\n", "t = t .& inpctile(ign[num], 99) .& inpctile(ign[denom], 99)\n", "\n", "# Resample, returning binned means and uncertainties \n", "# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)\n", "(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)\n", "\n", "# Plot results\n", "h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkred,markerstrokecolor=:darkred,label=\"\")\n", "plot!(h, xlabel=\"Age (Ma)\", ylabel=\"$(num) / $(denom)\",framestyle=:box,grid=:off,xflip=true) # Format plot\n", "display(h)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "## --- Ratio differentiation\n", "xelem = \"SiO2\"\n", "xmin = 40 # Minimum age\n", "xmax = 80 # Maximum age\n", "nbins = 20\n", "num = \"Sc\" # Numerator\n", "denom = \"Yb\" # Denominator\n", "\n", "# Exclude outliers\n", "t = inpctile(ign[num], 99) .& inpctile(ign[denom], 99)\n", "\n", "# Resample, returning binned means and uncertainties\n", "# (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)\n", "(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)\n", "\n", "# Plot results\n", "h = plot(c,m,yerror=(el,eu),seriestype=:scatter,color=:darkblue,markerstrokecolor=:darkblue,label=\"\")\n", "plot!(h, xlabel=xelem, ylabel=\"$(num) / $(denom)\",xlims=(xmin,xmax),framestyle=:box,grid=:off) # Format plot\n", "display(h)\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.0", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.0" } }, "nbformat": 4, "nbformat_minor": 2 }