/
SPD_curve.py
95 lines (73 loc) · 3.47 KB
/
SPD_curve.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
# -*- coding: utf-8 -*-
"""
File name: SPD_curve.py
Author: Mario Gutiérrez-Roig
Date created: 14/07/2016
Python Version: 3.4.3
Description: This code generates a plot of the Sum of Probabilities
Distribution (SPD) of a collection of dates stored in datasets called as
an input. It onky takes into account filtered dates according an
archaeological criteria. Dates with a certain percentage of marine source
carbon are calibrated with a mixture of atmospherical-marine curve.
"""
import numpy as np
import iosacal
import pkg_resources
import pickle
# Input data file path
input_file="Data/Iberia.csv"
# Output file name
output_file="Plots/SPD_Iberia_no_taph.pdf"
output_file_2="Simulations/SPD_Iberia_no_taph.pkl"
# Raw data taking some colums of the input file (fields should be delimited by ";")
raw = np.genfromtxt(input_file, delimiter=";", names=True, usecols=("ID_date", "C14Age", "Std", "Type_site", "Marine", "Marine_err", "Reservoir", "Reservoir_err", "Filtered"), dtype=('|S16', float, float, '|S16', float, float, float, float, int), skip_footer=0)
# Consider only the filtered dates
filtered = [[x[0], x[1], x[2], x[3]] for x in raw if x[4] == 0 and x[8] == 1]
marinedates = [[x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7]] for x in raw if x[4] > 0 and x[8] == 1]
# Radiocarbon dates seggregated by site type
radiocarbon_dates = [iosacal.core.RadiocarbonDetermination(x[1], x[2], x[0]) for x in filtered]
# Calibrate dates
cal_dates = [x.calibrate('intcal13', norm=True) for x in radiocarbon_dates]
# Performs the same analysis for marine dates
if len(marinedates) > 0:
# Calibrate Marine dates
calmarine = []
for x in marinedates:
r = iosacal.core.RadiocarbonDetermination(x[1], x[2], x[0])
# Define the main calibration curve
curve_name = "intcal13"
# Determine the path to the calibration curve file
curve_path = pkg_resources.resource_filename("iosacal", "data/%s.14c" % curve_name)
# CalibrationCurve format of iosacal
calibration_curve = iosacal.core.CalibrationCurve(curve_path)
# Mix it with secondary curve with a certain percentage P +/- D and Rerservoir effect deltaR +/- err_deltaR
calibration_curve.mixing("marine13", P=x[4], D=x[5], deltaR=x[6], err_deltaR=x[7])
cal_r = r.calibrate(calibration_curve)
cal_dates.append(cal_r)
# Calculate SPDs
SPD = iosacal.core.SPD(cal_dates, norm=False)
# x-range limits
maxx = 18000
minx = 7500
# Save final SPD curve into pickle file
with open(output_file_2, 'wb') as output:
pickle.dump(SPD, output, -1)
# Plot script #
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
minorLocator = AutoMinorLocator()
fig, ax = plt.subplots(figsize=(12,6))
plt.tick_params(labelsize=14)
plt.xlim(maxx,minx)
plt.xlabel("Calendar Age (BP)", fontsize=14)
plt.ylabel("Probability", fontsize=14)
ax.xaxis.set_minor_locator(minorLocator)
# Uncomment in case of log-scale plot
#plt.ylim(5e-6,5e-4)
#plt.yscale('log')
#This plots red and green vertical lines at 8240 and 9300
#plt.plot([8240,8240], [0.0006,0], '-', linewidth=1, color='red')
#plt.plot([9300,9300], [0.0006,0], '-', linewidth=4, color='green')
plt.fill_between(SPD.T[0], SPD.T[1], facecolor="blue", alpha="0.25", edgecolor="none")
plt.plot(SPD.T[0], SPD.T[1], '-', linewidth=1, color='black')
plt.savefig(output_file, format='pdf')