Skip to content

Develop #30

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
davewhipp merged 2 commits into from
Feb 27, 2023
Merged

Develop #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e1c4ee6
Update handling of measured age data file to provide better user feed…
davewhipp Feb 27, 2023
3b6a855
Update version to 0.3.2
davewhipp Feb 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from setuptools import setup

setup(name='tcplotter',
version='0.3.1',
version='0.3.2',
description='Plots thermochronometer ages and closure temperatures',
url='https://github.com/HUGG/tcplotter',
author='David Whipp',
Expand Down
186 changes: 109 additions & 77 deletions tcplotter/tcplotter.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,8 @@ def get_tc_exec(command):


# Function for reading age data file
def read_age_data(file):
"""Reads in age data from a csv file"""
def read_age_data(file, ap_x, ap_y, zr_x, zr_y):
"""Reads in age data from a csv file and stores plotable data"""
# Make empty lists for column values
ahe_age = []
ahe_uncertainty = []
Expand All @@ -54,35 +54,72 @@ def read_age_data(file):
if (len(data[i][3]) > 0) and (len(data[i][4]) > 0):
# Append AHe data if the age type is AHe
if data[i][0].lower() == "ahe":
ahe_age.append(float(data[i][1]))
ahe_uncertainty.append(float(data[i][2]))
ahe_eu.append(float(data[i][3]))
ahe_radius.append(float(data[i][4]))
# Check if value is not within eU range
if not (min(ap_x) <= float(data[i][3]) <= max(ap_x)):
print(
f"Warning: eU concentration for {data[i][0].lower()} age on line {i + 1} not within modelled range."
)
print(
f" This age will be excluded from plotting and the misfit calculation."
)
# Check if value is not within radius range
elif not (min(ap_y) <= float(data[i][4]) <= max(ap_y)):
print(
f"Warning: Grain radius for {data[i][0].lower()} age on line {i + 1} not within modelled range."
)
print(
f" This age will be excluded from plotting and the misfit calculation."
)
else:
ahe_age.append(float(data[i][1]))
ahe_uncertainty.append(float(data[i][2]))
ahe_eu.append(float(data[i][3]))
ahe_radius.append(float(data[i][4]))

# Append ZHe data if the age type is ZHe
elif data[i][0].lower() == "zhe":
zhe_age.append(float(data[i][1]))
zhe_uncertainty.append(float(data[i][2]))
zhe_eu.append(float(data[i][3]))
zhe_radius.append(float(data[i][4]))
# Check if value is not within eU range
if not (min(zr_x) <= float(data[i][3]) <= max(zr_x)):
print(
f"Warning: eU concentration for {data[i][0].lower()} age {i + 1} not within modelled range."
)
print(
f" This age will be excluded from plotting and the misfit calculation."
)
# Check if value is not within radius range
elif not (min(zr_y) <= float(data[i][4]) <= max(zr_y)):
print(
f"Warning: Grain radius for {data[i][0].lower()} age {i + 1} not within modelled range."
)
print(
f" This age will be excluded from plotting and the misfit calculation."
)
else:
zhe_age.append(float(data[i][1]))
zhe_uncertainty.append(float(data[i][2]))
zhe_eu.append(float(data[i][3]))
zhe_radius.append(float(data[i][4]))
else:
if (data[i][0].lower() != "ahe") and (data[i][0].lower() != "zhe"):
print(
f"Warning: Unsupported age type ({data[i][0].lower()}) on line {i + 1}."
)
elif len(data[i][3]) == 0:
print(
f"Warning: {data[i][0].lower()} age on line {i + 1} of age file missing eU value."
f"Warning: {data[i][0].lower()} age on line {i + 1} missing eU value."
)
elif len(data[i][4]) == 0:
print(
f"Warning: {data[i][0].lower()} age on line {i + 1} of age file missing radius value."
f"Warning: {data[i][0].lower()} age on line {i + 1} missing radius value."
)
print(f" Age will not be plotted.")
# Create new lists with data file values
n_ahe = len(ahe_age)
n_zhe = len(zhe_age)
ahe_data = [ahe_age, ahe_uncertainty, ahe_eu, ahe_radius]
zhe_data = [zhe_age, zhe_uncertainty, zhe_eu, zhe_radius]

return ahe_data, zhe_data
return n_ahe, ahe_data, n_zhe, zhe_data


def chi_squared(observed, expected, std):
Expand All @@ -94,7 +131,7 @@ def chi_squared(observed, expected, std):
return misfit / len(observed)


def calculate_misfit(age_data, age_type, age_list, param_x, param_y):
def calculate_misfit(age_data, age_list, param_x, param_y):
"""Calculates misfit between measured and predicted ages."""
predicted_ages = []
measured_ages = []
Expand All @@ -103,27 +140,13 @@ def calculate_misfit(age_data, age_type, age_list, param_x, param_y):
# Create interpolation function
age_grid = np.array(age_list).reshape((len(param_x), len(param_y)))
age_interp = RegularGridInterpolator((param_x, param_y), age_grid)
# Append interpolated age if within the eU and radius ranges
# Check if value is not within eU range
if not (min(param_x) <= age_data[2][i] <= max(param_x)):
print(
f"Warning: eU concentration for {age_type} age {i + 1} not within modelled range."
)
print(f" This age will be excluded from the misfit calculation.")
elif not (min(param_y) <= age_data[3][i] <= max(param_y)):
print(
f"Warning: Grain radius for {age_type} age {i + 1} not within modelled range."
)
print(f" This age will be excluded from the misfit calculation.")
else:
predicted_ages.append(age_interp([age_data[2][i], age_data[3][i]]))
# Include only measured ages within the eU/radius ranges
measured_ages.append(age_data[0][i])
std_dev.append(age_data[1][i])
predicted_ages.append(age_interp([age_data[2][i], age_data[3][i]]))
measured_ages.append(age_data[0][i])
std_dev.append(age_data[1][i])
# Calculate misfit
misfit = chi_squared(measured_ages, predicted_ages, std_dev)
n_ages = len(measured_ages)
return misfit[0], n_ages
return misfit[0]


# Define function for creating plot of cooling rates
Expand Down Expand Up @@ -400,13 +423,6 @@ def eu_vs_radius(
)
use_widget = False

# Read in measured ages from file, if age_data_file is defined
if len(age_data_file) > 0:
ahe_age_data, zhe_age_data = read_age_data(age_data_file)
else:
ahe_age_data = []
zhe_age_data = []

# Ensure relative paths work by setting working dir to dir containing this script file
wd_orig = os.getcwd()
script_path = os.path.abspath(__file__)
Expand Down Expand Up @@ -462,35 +478,46 @@ def eu_vs_radius(
# Set plot style
plt.style.use(plot_style)

# Set plot loop variables
ap_x = ap_eu
ap_y = ap_rad
zr_x = zr_eu
zr_y = zr_rad

# Read in measured ages from file, if age_data_file is defined
if len(age_data_file) > 0:
print(f"Reading age data file {age_data_file}...")
fp = os.path.join(wd_orig, age_data_file)
num_ahe_age_data, ahe_age_data, num_zhe_age_data, zhe_age_data = read_age_data(
fp, ap_x, ap_y, zr_x, zr_y
)
print(
f"Found {num_ahe_age_data} usable AHe age(s) and {num_zhe_age_data} usable ZHe age(s)."
)

# Define plot size and number of subplots
fig_width = 10
if plot_type < 3:
fig_height = 5
# Make plot longer if plotting data
if plot_type == 1:
if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
fig_height += 1.5
else:
if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
fig_height += 1.5
# Create figure and axes
fig, ax = plt.subplots(1, 2, figsize=(fig_width, fig_height))
else:
fig_height = 10
# Make plot longer if plotting data
if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
fig_height += 1.5
if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
fig_height += 1.5
# Create figure and axes
fig, ax = plt.subplots(2, 2, figsize=(fig_width, fig_height))

# Set plot loop variables
ap_x = ap_eu
ap_y = ap_rad
zr_x = zr_eu
zr_y = zr_rad

# Create lists for storing closure temperatures, ages
ahe_tc_list = []
ahe_age_list = []
Expand Down Expand Up @@ -596,20 +623,25 @@ def eu_vs_radius(
os.remove(tt_file)

# Calculate age misfits if age data file is used and option is selected
if calc_misfit and len(age_data_file) > 0:
total_misfit = 0
if len(ahe_age_data) > 0:
ahe_misfit, ahe_n_ages = calculate_misfit(
ahe_age_data, "AHe", ahe_age_list, ap_x, ap_y
if calc_misfit:
if (num_ahe_age_data > 0) or (num_zhe_age_data > 0):
total_misfit = 0
if num_ahe_age_data > 0:
ahe_misfit = calculate_misfit(ahe_age_data, ahe_age_list, ap_x, ap_y)
print(f"AHe misfit (n = {num_ahe_age_data} ages): {ahe_misfit}")
total_misfit += ahe_misfit
if num_zhe_age_data > 0:
zhe_misfit = calculate_misfit(zhe_age_data, zhe_age_list, zr_x, zr_y)
print(f"ZHe misfit (n = {num_zhe_age_data} ages): {zhe_misfit}")
total_misfit += zhe_misfit
print(
f"Total misfit (n = {num_ahe_age_data + num_zhe_age_data} ages): {total_misfit}"
)
print(f"AHe misfit (n = {ahe_n_ages} ages): {ahe_misfit}")
if len(zhe_age_data) > 0:
zhe_misfit, zhe_n_ages = calculate_misfit(
zhe_age_data, "ZHe", zhe_age_list, zr_x, zr_y
else:
print(
f"Warning: Misfit calculation requested, but no AHe or ZHe ages found in data file."
)
print(f"ZHe misfit (n = {zhe_n_ages} ages): {zhe_misfit}")
total_misfit += ahe_misfit + zhe_misfit
print(f"Total misfit (n = {ahe_n_ages + zhe_n_ages} ages): {total_misfit}")
print(f" Misfit will not be calculated.")

# Apatite eU versus radius
if plot_type == 1:
Expand All @@ -625,7 +657,7 @@ def eu_vs_radius(
# Add age contour labels
ax[0].clabel(ap_contours_age)
# Determine bounds for contour colors if plotting age data
if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
age_min = min(min(ahe_age_list), min(ahe_age_data[0]))
age_max = max(max(ahe_age_list), max(ahe_age_data[0]))
else:
Expand All @@ -642,8 +674,8 @@ def eu_vs_radius(
vmax=age_max,
alpha=plot_alpha,
)
if len(ahe_age_data) > 0:
ahe_label = f"AHe data (n = {ahe_n_ages}"
if num_ahe_age_data > 0:
ahe_label = f"AHe data (n = {num_ahe_age_data}"
if calc_misfit:
ahe_label += f"; misfit: {ahe_misfit:.2f}"
ahe_label += ")"
Expand Down Expand Up @@ -694,7 +726,7 @@ def eu_vs_radius(
alpha=plot_alpha,
)

if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
# Plot the colorbar if plotting data
fig.colorbar(
ap_contourf_tc,
Expand All @@ -721,7 +753,7 @@ def eu_vs_radius(
# Add age contour labels
ax[0].clabel(zr_contours_age)
# Determine bounds for contour colors if plotting age data
if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
age_min = min(min(zhe_age_list), min(zhe_age_data[0]))
age_max = max(max(zhe_age_list), max(zhe_age_data[0]))
else:
Expand All @@ -738,8 +770,8 @@ def eu_vs_radius(
vmax=age_max,
alpha=plot_alpha,
)
if len(zhe_age_data) > 0:
zhe_label = f"ZHe data (n = {zhe_n_ages}"
if num_zhe_age_data > 0:
zhe_label = f"ZHe data (n = {num_zhe_age_data}"
if calc_misfit:
zhe_label += f"; misfit: {zhe_misfit:.2f}"
zhe_label += ")"
Expand Down Expand Up @@ -790,7 +822,7 @@ def eu_vs_radius(
alpha=plot_alpha,
)

if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
# Plot the colorbar if plotting data
fig.colorbar(
zr_contourf_tc,
Expand All @@ -817,7 +849,7 @@ def eu_vs_radius(
# Add age contour labels
ax[0][0].clabel(ap_contours_age)
# Determine bounds for contour colors if plotting age data
if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
age_min = min(min(ahe_age_list), min(ahe_age_data[0]))
age_max = max(max(ahe_age_list), max(ahe_age_data[0]))
else:
Expand All @@ -834,8 +866,8 @@ def eu_vs_radius(
vmax=age_max,
alpha=plot_alpha,
)
if len(ahe_age_data) > 0:
ahe_label = f"AHe data (n = {ahe_n_ages}"
if num_ahe_age_data > 0:
ahe_label = f"AHe data (n = {num_ahe_age_data}"
if calc_misfit:
ahe_label += f"; misfit: {ahe_misfit:.2f}"
ahe_label += ")"
Expand Down Expand Up @@ -886,7 +918,7 @@ def eu_vs_radius(
alpha=plot_alpha,
)

if len(ahe_age_data) > 0:
if num_ahe_age_data > 0:
# Plot the colorbar if plotting data
fig.colorbar(
ap_contourf_tc,
Expand All @@ -911,7 +943,7 @@ def eu_vs_radius(
# Add age contour labels
ax[1][0].clabel(zr_contours_age)
# Determine bounds for contour colors if plotting age data
if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
age_min = min(min(zhe_age_list), min(zhe_age_data[0]))
age_max = max(max(zhe_age_list), max(zhe_age_data[0]))
else:
Expand All @@ -928,8 +960,8 @@ def eu_vs_radius(
vmax=age_max,
alpha=plot_alpha,
)
if len(zhe_age_data) > 0:
zhe_label = f"ZHe data (n = {zhe_n_ages}"
if num_zhe_age_data > 0:
zhe_label = f"ZHe data (n = {num_zhe_age_data}"
if calc_misfit:
zhe_label += f"; misfit: {zhe_misfit:.2f}"
zhe_label += ")"
Expand Down Expand Up @@ -980,7 +1012,7 @@ def eu_vs_radius(
alpha=plot_alpha,
)

if len(zhe_age_data) > 0:
if num_zhe_age_data > 0:
# Plot the colorbar if plotting data
fig.colorbar(
zr_contourf_tc,
Expand Down