Fork me on GitHub

Plot Thermocalc avPT results with Python

Plot Thermocalc avPT results, including error ellipses with Python.

png

The figure above plots the following PT results saved as file called avPT.csv. The data were taken from Dave Waters' Practical Aspects of Thermobarometry.

Sample avP sdP avT sdT cor
No 1 8.5 1.2 751 80 0.922
No 2 5.2 1.0 721 42 -0.500
No 3 9.2 0.9 650 60 0.883
No 4 2.5 0.9 686 135 0.869
No 5 5.6 0.5 480 60 0.750
No 6 3.5 1.0 520 50 0.400

The Python code to plot Thermocalc results as shown above.

# setup python environment
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# define ellipse functions
def eigsorted(cov):
        vals, vecs = np.linalg.eigh(cov)
        order = vals.argsort()[::-1]
        return vals[order], vecs[:,order]

def get_ellipse_coords(a=0.0, b=0.0, x=0.0, y=0.0, angle=0.0, k=1):
    """ Draws an ellipse using (360*k + 1) discrete points
    # Source: http://central.scipy.org/item/23/2/plot-an-ellipse
    k = 1 means 361 points (degree by degree)
    a = major axis distance,
    b = minor axis distance,
    x = offset along the x-axis
    y = offset along the y-axis
    angle = clockwise rotation [in degrees] of the ellipse;
        * angle=0  : the ellipse is aligned with the positive x-axis
        * angle=30 : rotated 30 degrees clockwise from positive x-axis
    """
    pts = np.zeros((360*k+1, 2))

    beta = -angle * np.pi/180.0
    sin_beta = np.sin(beta)
    cos_beta = np.cos(beta)
    alpha = np.radians(np.r_[0.:360.:1j*(360*k+1)])

    sin_alpha = np.sin(alpha)
    cos_alpha = np.cos(alpha)

    pts[:, 0] = x + (a * cos_alpha * cos_beta - b * sin_alpha * sin_beta)
    pts[:, 1] = y + (a * cos_alpha * sin_beta + b * sin_alpha * cos_beta)

    return pts

# import data
# data from https://www.earth.ox.ac.uk/~davewa/pt/th_tools.html avpt.xls
df = pd.read_csv('avPT.csv')

# aluminosilicate reactions
p_sk = (4.4, 6.0, 8.0, 10.0, 12.0, 14.0)
t_sk = (550, 623, 715, 806, 896, 984)
p_ak = (2.0, 4.0, 4.4)                                                     
t_ak = (352, 514, 550)
p_as = (2.0, 4.4)
t_as = (717, 550)
plt.plot(t_sk, p_sk, 'k-', lw=1)
plt.plot(t_ak, p_ak, 'k-', lw=1)
plt.plot(t_as, p_as, 'k-', lw=1)

# plot pt ellipses
for i in range(0, len(df.avT)):
    sxy = df.cor[i] * df.sdT[i] * df.sdP[i]
    cov = np.array([[df.sdT[i]**2, sxy], [sxy, df.sdP[i]**2]])
    width = np.sqrt(max(np.linalg.eigvalsh(cov)))
    height = np.sqrt(min(np.linalg.eigvalsh(cov)))
    vals, vecs = eigsorted(cov)
    theta = -1*np.degrees(np.arctan2(*vecs[:,0][::-1]))
    ellip = get_ellipse_coords(a=width, b=height, x=df.avT[i], y=df.avP[i], angle=theta, k=2)

    plt.plot(ellip[:,0], ellip[:, 1], lw=1)
    plt.scatter(df.avT[i], df.avP[i], label=df.Sample[i])

# add plot labels
plt.xlim(300, 900)
plt.ylim(0, 14)
plt.xlabel('Temperature (C)')
plt.ylabel('Pressure (kbar)')
plt.legend(title='Sample')

Comments !

social