An open source radiocarbon (14C) calibration software. Written in Python http://c14.iosa.it/
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 

339 lines
10 KiB

#! /usr/bin/env python
# -*- coding: utf-8 -*-
# filename: plot.py
# Copyright 2018 Stefano Costa <steko@iosa.it>
# Copyright 2020 Karl Håkansson <kalle@jomenvisst.de>
#
# This file is part of IOSACal, the IOSA Radiocarbon Calibration Library.
# IOSACal is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# IOSACal is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with IOSACal. If not, see <http://www.gnu.org/licenses/>.
import matplotlib.pyplot as plt
import iosacal
from scipy.stats import norm
COLORS = {
'bgcolor': '#e5e4e5',
}
def single_plot(calibrated_age, oxcal=False, output=None, BP='bp'):
calibrated_age = calibrated_age
f_m = calibrated_age.radiocarbon_sample.date
sigma_m = calibrated_age.radiocarbon_sample.sigma
radiocarbon_sample_id = calibrated_age.radiocarbon_sample.id
calibration_curve = calibrated_age.calibration_curve
intervals = calibrated_age.intervals
sample_interval = calibration_curve[:,0].copy() # for determination plot
# adjust plot bounds
min_year, max_year = (50000, -50000)
min_x = min(calibrated_age[:,0])
max_x = max(calibrated_age[:,0])
if min_year > min_x:
min_year = min_x
if max_year < max_x:
max_year = max_x
# do not plot the part of calibration curve that is not visible
# greatly reduces execution time \o/
cutmin = calibration_curve[calibration_curve[:,0]>min_x]
cutmax = cutmin[cutmin[:,0]<max_x]
calibration_curve = cutmax
if BP != 'bp':
if min_year < 0 and max_year > 0:
ad_bp_label = "BC/AD"
elif min_year < 0 and max_year < 0:
ad_bp_label = "BC"
elif min_year > 0 and max_year > 0:
ad_bp_label = "AD"
else:
ad_bp_label = "BP"
string68 = '{:{fmt}}'.format(calibrated_age.intervals[68], fmt=BP)
string95 = '{:{fmt}}'.format(calibrated_age.intervals[95], fmt=BP)
fig = plt.figure(figsize=(12,8))
fig.clear()
ax1 = plt.subplot(111)
ax1.set_facecolor(COLORS['bgcolor'])
plt.xlabel("Calibrated age ({})".format(ad_bp_label))
plt.ylabel("Radiocarbon determination (BP)")
plt.text(0.5, 0.95,
r'{}: {:.0f} ± {:.0f} BP'.format(radiocarbon_sample_id, f_m, sigma_m),
horizontalalignment='center',
verticalalignment='center',
transform = ax1.transAxes,
bbox=dict(facecolor='white', alpha=0.9, lw=0))
plt.text(0.75, 0.80,
'68.2% probability\n{}\n\n95.4% probability\n{}'.format(string68, string95),
horizontalalignment='left',
verticalalignment='center',
transform = ax1.transAxes,
bbox=dict(facecolor='white', alpha=0.9, lw=0))
plt.text(0.0, 1.0,'IOSACal v{}; {}'.format(iosacal.__VERSION__, calibration_curve.title),
horizontalalignment='left',
verticalalignment='bottom',
transform = ax1.transAxes,
size=10,
bbox=dict(facecolor='white', alpha=0.9, lw=0))
# Calendar Age
ax2 = plt.twinx()
if oxcal is True:
# imitate OxCal
ax1.set_facecolor('white')
ax2.fill_between(
calibrated_age[:,0],
calibrated_age[:,1]*0,
calibrated_age[:,1],
facecolor='k',
alpha=0.3,
label='Calendar Age'
)
ax2.plot(
calibrated_age[:,0],
calibrated_age[:,1],
'k',
alpha=0
)
else:
ax2.fill_between(
calibrated_age[:,0],
calibrated_age[:,1]*0,
calibrated_age[:,1],
facecolor='k',
alpha=0.3,
label='Calendar Age'
)
ax2.plot(
calibrated_age[:,0],
calibrated_age[:,1],
'k',
alpha=0
)
ax2.set_ybound(min(calibrated_age[:,1]),max(calibrated_age[:,1])*3)
ax2.set_xbound(min(calibrated_age[:,0]),max(calibrated_age[:,0]))
ax2.set_axis_off()
# Radiocarbon Age
sample_curve = norm.pdf(sample_interval, f_m, sigma_m)
if oxcal is True:
sample_fill_color = '#fac5cd'
else:
sample_fill_color = 'w'
ax3 = plt.twiny(ax1)
ax3.fill(
sample_curve,
sample_interval,
'1.0',
facecolor=sample_fill_color,
alpha=0.8
)
ax3.set_xbound(0,max(sample_curve)*4)
ax3.set_axis_off()
# Calibration Curve
curve_low = calibration_curve[:,1] - calibration_curve[:,2]
curve_high = calibration_curve[:,1] + calibration_curve[:,2]
ax1.fill_between(calibration_curve[:,0],
curve_low,
curve_high,
facecolor='#000000',
edgecolor='none',
alpha=0.15)
ax1.plot(calibration_curve[:,0], calibration_curve[:,1], 'k', lw=0.5, alpha=0.7)
# Confidence intervals
if oxcal is True:
for i in intervals[68]:
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0.05,
ymax=0.07,
facecolor='none',
alpha=0.8)
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0.068,
ymax=0.072,
facecolor='w',
edgecolor='w',
lw=2)
for i in intervals[95]:
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0.025,
ymax=0.045,
facecolor='none',
alpha=0.8)
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0.043,
ymax=0.047,
facecolor='w',
edgecolor='w',
lw=2)
else:
for i in intervals[68]:
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0,
ymax=0.02,
facecolor='k',
alpha=0.5)
for i in intervals[95]:
ax1.axvspan(
i.from_year,
i.to_year,
ymin=0,
ymax=0.02,
facecolor='k',
alpha=0.5)
# FIXME the following values 10 and 5 are arbitrary and could be probably
# drawn from the f_m value itself, while preserving their ratio
ax1.set_ybound(f_m - sigma_m * 15, f_m + sigma_m * 5)
ax1.set_xbound(min(calibrated_age[:,0]),max(calibrated_age[:,0]))
ax1.invert_xaxis()
if BP != 'bp':
ax1.set_xticklabels([int(item*-1+1950) for item in ax1.get_xticks().tolist()])
if output:
plt.savefig(output)
def stacked_plot(calibrated_ages,name='Stacked plot',oxcal=False, BP='bp', output=None):
'''Plot multiple calibrated ages, vertically stacked.
``calibrated_ages`` is an iterable of CalAge objects.'''
# Define the legend and descriptive text
min_year, max_year = (50000, -50000)
for calibrated_age in calibrated_ages:
radiocarbon_sample_id = calibrated_age.radiocarbon_sample.id
calibration_curve = calibrated_age.calibration_curve
calibration_curve_title = calibration_curve.title
intervals = calibrated_age.intervals
if min_year > min(calibrated_age[:,0]):
min_year = min(calibrated_age[:,0])
if max_year < max(calibrated_age[:,0]):
max_year = max(calibrated_age[:,0])
if BP != 'bp':
if min_year < 0 and max_year > 0:
ad_bp_label = "BC/AD"
elif min_year < 0 and max_year < 0:
ad_bp_label = "BC"
elif min_year > 0 and max_year > 0:
ad_bp_label = "AD"
else:
ad_bp_label = "BP"
numrows = len(calibrated_ages)
fig, axs = plt.subplots(numrows, 1, sharex=True, figsize=(12, 2*numrows))
plt.suptitle("{}".format(name))
if BP == 'bp':
axs[0].invert_xaxis() # just once, because the axis is shared
for n, calibrated_age in enumerate(calibrated_ages):
ax = axs[n]
# Calendar Age
ax.fill_between(
calibrated_age[:,0],
calibrated_age[:,1]*0,
calibrated_age[:,1],
facecolor='k',
alpha=0.3,
label='Calendar Age'
)
ax.plot(
calibrated_age[:,0],
calibrated_age[:,1],
'k',
alpha=0.6
)
ax.set_ybound(
min(calibrated_age[:,1]),
max(calibrated_age[:,1])*2
)
ax.set_xbound(min_year, max_year)
# remove labels for Y axis - values are meaningless
ax.get_yaxis().set_ticklabels('')
# Legend
plt.text(0.95, 0.85,'{0.radiocarbon_sample.id!s}'.format(calibrated_age),
horizontalalignment='right',
verticalalignment='center',
transform = ax.transAxes,
bbox=dict(facecolor='white', alpha=0.9, lw=0))
# Confidence intervals
for i in calibrated_age.intervals[95]:
ax.axvspan(
i.from_year,
i.to_year,
ymin=0.6,
ymax=0.7,
facecolor='k',
alpha=0.4)
for i in calibrated_age.intervals[68]:
ax.axvspan(
i.from_year,
i.to_year,
ymin=0.6,
ymax=0.7,
facecolor='k',
alpha=1.0)
if n + 1 == numrows:
if BP != 'bp':
plt.xlabel("Calibrated age ({})".format(ad_bp_label), y = 0.05) # out of plot boundaries :(
else:
plt.xlabel("Calibrated age (BP)", y = 0.05)
if output:
plt.savefig(output)
def iplot(calibrated_ages, **kwds):
'''A generic function for plotting in the IPython Notebook.'''
try:
# ugly, ugly hack
calibrated_ages.intervals
except AttributeError:
stacked_plot(calibrated_ages, **kwds)
else:
single_plot(calibrated_ages, **kwds)