This Jupyter notebook demonstrates eruption (/deposition) age estimation component of Chron.jl
, based on the approach of Keller, Schoene, and Samperton (2018). For more information on Chron.jl
and its other capabilities (largely, age-depth modelling), see github.com/brenhinkeller/Chron.jl and and doi.org/10.17605/osf.io/TQX3F
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 Chron.jl package
using Chron
using Plots, DelimitedFiles
First, let's enter some basic information about your samples. How many are there, what are the sample names (needs to be a valid filename, BTW), and what are the stratigraphic heights and uncertainties? Then, we'll enter the ages as .csv files.
nSamples = 5 # The number of samples you have data for
smpl = ChronAgeData(nSamples)
smpl.Name = ("KJ08-157", "KJ04-75", "KJ09-66", "KJ04-72", "KJ04-70",)
smpl.Path = "MyData/" # Where are the data files? Must match where you put the csv files below.
smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.
smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files
smpl.Height_Unit = "cm"; # Unit of measurement for Height and Height_sigma
In this case, there is no need to specify sample heights, since we are not doing any age-depth modelling
Now let's see what's in our current directory (we'll use ;
to activate Julia's command-line shell mode, followed by a unix command, in this case ls
;ls
Chron1.0Coupled.ipynb Chron1.0Coupled.jl Chron1.0CoupledConcordia.ipynb Chron1.0CoupledConcordia.jl Chron1.0CoupledSystematic.jl Chron1.0DistOnly.ipynb Chron1.0DistOnly.jl Chron1.0Radiocarbon.ipynb Chron1.0Radiocarbon.jl Chron1.0StratOnly.ipynb Chron1.0StratOnly.jl ConcordiaExampleData DenverUPbExampleData EruptionDepositionAgeDemonstration.ipynb EruptionDepositionAgeDemonstration.jl FCT Manifest.toml MyData PlutonEmplacement.ipynb Project.toml archive.tar.gz ffsend
Equivalently, we can also run unix commands using the run()
function, i.e.
run(`ls`);
Now that we know how to access the command line, let's make a new folder for our example data. This can be called whatever you want, just make sure it matches smpl.Path
above
;mkdir -p MyData/
Now, let's make some files and paste in csv data for each one. For now, I'm pasting in example data from Clyde et al. (2016), 10.1016/j.epsl.2016.07.041
# You can just paste your data in here, in the following two-column format.
# The first column is age and the second column is two-sigma analytical uncertainty.
# You should generally exclude all systematic uncertainties here.
data = [
66.12 0.14
66.115 0.048
66.11 0.1
66.11 0.17
66.096 0.056
66.088 0.081
66.085 0.076
66.073 0.084
66.07 0.11
66.055 0.043
66.05 0.16
65.97 0.12
]
# Now, let's write this data to a file, delimited by commas (',')
# In this example the filename is KJ08-157.csv, in the folder called MyData
writedlm("MyData/KJ08-157.csv", data, ',')
data = [
66.24 0.25
66.232 0.046
66.112 0.085
66.09 0.1
66.04 0.18
66.03 0.12
66.016 0.08
66.003 0.038
65.982 0.071
65.98 0.19
65.977 0.042
65.975 0.066
65.971 0.082
65.963 0.074
65.92 0.12
65.916 0.088
]
writedlm("MyData/KJ04-75.csv",data,',')
data = [
66.13 0.15
66.066 0.052
65.999 0.045
65.989 0.057
65.98 0.11
65.961 0.057
65.957 0.091
65.951 0.066
65.95 0.11
65.929 0.059
]
writedlm("MyData/KJ09-66.csv",data,',')
data = [
66.11 0.2
66.003 0.063
66.003 0.058
65.98 0.06
65.976 0.089
65.973 0.084
65.97 0.15
65.963 0.055
65.959 0.049
65.94 0.18
65.928 0.066
65.92 0.057
65.91 0.14
]
writedlm("MyData/KJ04-72.csv",data,',')
data = [
66.22 0.27
66.06 0.11
65.933 0.066
65.918 0.087
65.92 0.34
65.916 0.067
65.91 0.18
65.892 0.09
65.89 0.063
65.89 0.15
65.88 0.2
65.812 0.069
65.76 0.15
]
writedlm("MyData/KJ04-70.csv",data,',')
Alternatively, you could download .csv files that you have posted somewhere online using the Julia download()
function, or using the unix command curl
throught the command-line interface
For each sample in smpl.Name
, we must have a .csv
file in smpl.Path
which contains each individual mineral age and uncertainty.
To learn more about the eruption/deposition age estimation model, see also Keller, Schoene, and Sameperton (2018) and the BayeZirChron demo notebook. It is important to note that this model (like most if not all others) has no knowledge of open-system behaviour, so e.g., Pb-loss will lead to erroneous results.
# Bootstrap a KDE of the pre-eruptive (or pre-deposition) zircon distribution
# shape from individual sample datafiles using a KDE of stacked sample data
BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
h = plot(range(0,1,length=length(BootstrappedDistribution)), BootstrappedDistribution,
label="Bootstrapped distribution", xlabel="Time (arbitrary units)", ylabel="Probability Density",
fg_color_legend=:white, framestyle=:box)
savefig(h, joinpath(smpl.Path,"BootstrappedDistribution.pdf"))
display(h)
┌ Info: Interpreting first two columns of KJ08-157.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: Interpreting first two columns of KJ04-75.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: Interpreting first two columns of KJ09-66.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: Interpreting first two columns of KJ04-72.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: Interpreting first two columns of KJ04-70.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 |
# Number of steps to run in distribution MCMC
distSteps = 10^6
distBurnin = distSteps÷100
# Choose the form of the prior closure/crystallization distribution to use
dist = BootstrappedDistribution
## You might alternatively consider:
# dist = UniformDistribution # A reasonable default
# dist = MeltsVolcanicZirconDistribution # A single magmatic pulse, truncated by eruption
# dist = ExponentialDistribution # Applicable for survivorship processes, potentially including inheritance/dispersion in Ar-Ar dates
# Run MCMC to estimate saturation and eruption/deposition age distributions
smpl = tMinDistMetropolis(smpl,distSteps,distBurnin,dist)
results = vcat(["Sample" "Age" "2.5% CI" "97.5% CI" "sigma"], hcat(collect(smpl.Name),smpl.Age,smpl.Age_025CI,smpl.Age_975CI,smpl.Age_sigma))
writedlm(smpl.Path*"results.csv", results, ',');
[ Info: Estimating eruption/deposition age distributions... ┌ Info: 1: KJ08-157 │ Interpreting first two columns of KJ08-157.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: 2: KJ04-75 │ Interpreting first two columns of KJ04-75.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: 3: KJ09-66 │ Interpreting first two columns of KJ09-66.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: 4: KJ04-72 │ Interpreting first two columns of KJ04-72.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 | ┌ Info: 5: KJ04-70 │ Interpreting first two columns of KJ04-70.csv as └ | Age | Age 2-sigma 𝑎𝑏𝑠𝑜𝑙𝑢𝑡𝑒 |
# Let's see what that did
run(`ls $(smpl.Path)`);
results = readdlm(smpl.Path*"results.csv",',')
AgeDepthModel.pdf AgeDepthModel.svg BootstrappedDistribution.pdf DepositionRateModelCI.pdf DepositionRateModelCI.svg KJ04-70.csv KJ04-70_distribution.pdf KJ04-70_distribution.svg KJ04-70_rankorder.pdf KJ04-70_rankorder.svg KJ04-72.csv KJ04-72_distribution.pdf KJ04-72_distribution.svg KJ04-72_rankorder.pdf KJ04-72_rankorder.svg KJ04-75.csv KJ04-75_distribution.pdf KJ04-75_distribution.svg KJ04-75_rankorder.pdf KJ04-75_rankorder.svg KJ08-157.csv KJ08-157_distribution.pdf KJ08-157_distribution.svg KJ08-157_rankorder.pdf KJ08-157_rankorder.svg KJ09-66.csv KJ09-66_distribution.pdf KJ09-66_distribution.svg KJ09-66_rankorder.pdf KJ09-66_rankorder.svg KR18-01.csv KR18-01_Concordia.pdf KR18-01_Concordia.svg KR18-01_Pbloss.pdf KR18-01_Pbloss.svg KR18-01_distribution.pdf KR18-01_distribution.svg KR18-04.csv KR18-04_Concordia.pdf KR18-04_Concordia.svg KR18-04_Pbloss.pdf KR18-04_Pbloss.svg KR18-04_distribution.pdf KR18-04_distribution.svg KR18-05.csv KR18-05_Concordia.pdf KR18-05_Concordia.svg KR18-05_Pbloss.pdf KR18-05_Pbloss.svg KR18-05_distribution.pdf KR18-05_distribution.svg agedist.csv distresults.csv height.csv lldist.csv results.csv
6×5 Matrix{Any}: "Sample" "Age" "2.5% CI" "97.5% CI" "sigma" "KJ08-157" 66.0717 66.038 66.0938 0.0140167 "KJ04-75" 65.9319 65.888 65.9693 0.0209099 "KJ09-66" 65.9367 65.8943 65.9763 0.0205472 "KJ04-72" 65.9565 65.9258 65.9768 0.0128253 "KJ04-70" 65.8314 65.7577 65.8938 0.0347877
To see what one of these plots looks like, try pasting this into a markdown cell like the one below
<img src="MyData/KJ04-75_rankorder.svg" align="left" width="600"/>
For each sample, the eruption/deposition age distribution is inherently asymmetric, because of the one-sided relationship between mineral closure and eruption/deposition. For example (KJ04-70):
(if no figure appears, double-click to enter this markdown cell and re-evaluate (shift
-enter
) after running the model above
As before, we can use the unix command ls
to see all the files we have written. Actually getting them out of here may be harder though.
;ls MyData
AgeDepthModel.pdf AgeDepthModel.svg BootstrappedDistribution.pdf DepositionRateModelCI.pdf DepositionRateModelCI.svg KJ04-70.csv KJ04-70_distribution.pdf KJ04-70_distribution.svg KJ04-70_rankorder.pdf KJ04-70_rankorder.svg KJ04-72.csv KJ04-72_distribution.pdf KJ04-72_distribution.svg KJ04-72_rankorder.pdf KJ04-72_rankorder.svg KJ04-75.csv KJ04-75_distribution.pdf KJ04-75_distribution.svg KJ04-75_rankorder.pdf KJ04-75_rankorder.svg KJ08-157.csv KJ08-157_distribution.pdf KJ08-157_distribution.svg KJ08-157_rankorder.pdf KJ08-157_rankorder.svg KJ09-66.csv KJ09-66_distribution.pdf KJ09-66_distribution.svg KJ09-66_rankorder.pdf KJ09-66_rankorder.svg KR18-01.csv KR18-01_Concordia.pdf KR18-01_Concordia.svg KR18-01_Pbloss.pdf KR18-01_Pbloss.svg KR18-01_distribution.pdf KR18-01_distribution.svg KR18-04.csv KR18-04_Concordia.pdf KR18-04_Concordia.svg KR18-04_Pbloss.pdf KR18-04_Pbloss.svg KR18-04_distribution.pdf KR18-04_distribution.svg KR18-05.csv KR18-05_Concordia.pdf KR18-05_Concordia.svg KR18-05_Pbloss.pdf KR18-05_Pbloss.svg KR18-05_distribution.pdf KR18-05_distribution.svg agedist.csv distresults.csv height.csv lldist.csv results.csv
We could use the trick we learned before to view the SVG files in markdown cells, which you should then be able to right click and download as real vector graphics. e.g. pasting something like
<img src="AgeDepthModel.svg" align="left" width="600"/>
in a markdown cell such as this one
Meanwhile, for the csv files we could try something like ; cat agedist.csv
, but agedist is probably too big to print. Let's try using ffsend instead, which should give you a download link. In fact, while we're at it, we might as well archive and zip the entire directory!
# Make gzipped tar archive of the the whole MyData directory
run(`tar -zcf archive.tar.gz ./MyData`);
# Download prebuilt ffsend linux binary
isfile("ffsend") || download("https://github.com/timvisee/ffsend/releases/download/v0.2.65/ffsend-v0.2.65-linux-x64-static", "ffsend")
# Make ffsend executable
run(`chmod +x ffsend`);
; ./ffsend upload archive.tar.gz
/bin/bash: ./ffsend: cannot execute binary file
You could alternatively use the ffsend command in this way to transfer individual files
Keep in mind that any changes you make to this online notebook won't persist after you close the tab (or after it times out) even if you save your changes! You have to either copy-paste or file
>Download as
a copy.