API tools faq
paste
Guest User

Untitled

a guest
Mar 17th, 2023
546
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 14.67 KB | None | 0 0
raw download clone embed print report
  1. #!/usr/bin/env python3
  2. # -*- coding: utf-8 -*-
  3. """
  4. This script determines what proxy data was used in the AD 1400–1449
  5. portion of the MBH98 temperature reconstruction by regressing it on
  6.  
  7. * the proxy data used in Wahl & Ammann 2007 (WA), and
  8. * an extended dataset from the Climategate archive.
  9.  
  10. Results:
  11.  
  12. * MBH98 and WA used different datasets.
  13. * The MBH98 proxy data is in the Climategate archive.
  14.  
  15. NOTE: It may take a while to download the data (about 150 MB).
  16.  
  17. Tested on Linux with Python 3.8.10 and these library versions:
  18.  
  19. NumPy         1.17.4
  20. SciPy         1.3.3
  21. pandas        0.25.3
  22. Matplotlib    3.1.2
  23. """
  24.  
  25. import os.path
  26. import tarfile
  27. import requests
  28. import numpy as np
  29. import scipy.signal as sps
  30. import pandas as pd
  31. import matplotlib.pyplot as plt
  32.  
  33.  
  34. # Get submatrix of x where the first column of x is in range t.
  35. def submatrix(x, t):
  36.     return x[(x[:, 0] >= t[0]) & (x[:, 0] <= t[-1]), :]
  37.  
  38.  
  39. # Fill in missing data in the last rows of matrix x.
  40. def fill_in_missing(x):
  41.     for i in range(1, x.shape[1]):
  42.         if np.isnan(x[-1, i]):
  43.             index = (~np.isnan(x[:, i])).nonzero()[0]
  44.             x[index[-1]:, i] = x[index[-1], i]
  45.  
  46.  
  47. # Download data from url.
  48. def get_data_from(url):
  49.     filename = url.split("/")[-1]
  50.     if not os.path.isfile(filename):
  51.         r = requests.get(url, stream=True)
  52.         with open(filename, "wb") as fd:
  53.             for chunk in r.iter_content(chunk_size=1024**2):
  54.                 fd.write(chunk)
  55.     return filename
  56.  
  57.  
  58. # Untar tar archive in directory.
  59. def untar(tar, directory):
  60.     if not os.path.isdir(directory):
  61.         with tarfile.TarFile(tar, "r") as td:
  62.             td.extractall(directory)
  63.  
  64.  
  65. # Get MBH98 temperature reconstruction.
  66. def get_mbh98_reconstruction():
  67.     # Get data.
  68.     url = "http://www.meteo.psu.edu/holocene/public_html/shared/research/MANNETAL98/nhmean.dat"
  69.     filename = get_data_from(url)
  70.  
  71.     # Load data.
  72.     mbh = np.genfromtxt(filename, skip_header=1)[:, [0, 1]]
  73.     return mbh
  74.  
  75.  
  76. # Get Wahl & Ammann 2007 emulation of MBH98.
  77. def get_wahl_and_ammann_reconstruction():
  78.     # Get archive.
  79.     url = "https://projects.ral.ucar.edu/rc4a/millennium/WA_MBH-emulation.tar"
  80.     tarname = get_data_from(url)
  81.     directory = tarname.split(".tar")[0]
  82.     untar(tarname, directory)
  83.    
  84.     # Load data.
  85.     filename = "WAscenarios/WA_original/WA_original/NHem.reconstruction.dat"
  86.     path = os.path.join(directory, filename)
  87.     wa = np.genfromtxt(path)
  88.     wa = wa[:, [0, 2]]
  89.     return wa
  90.  
  91.  
  92. # Get proxy data used in Wahl & Ammann 2007.
  93. def get_wahl_and_ammann_proxy_data():
  94.     # Get proxy list.
  95.     url_datalist = "http://www.meteo.psu.edu/holocene/public_html/shared/research/MANNETAL98/PROXY/datalist1400.dat"
  96.     filename_datalist = get_data_from(url_datalist)
  97.    
  98.     # Get proxy data.
  99.     url_data = "http://www.meteo.psu.edu/holocene/public_html/shared/research/MANNETAL98/PROXY/data1400.dat"
  100.     filename_data = get_data_from(url_data)
  101.    
  102.     # Load proxy list and edit PC names.
  103.     datalist = pd.read_table(filename_datalist)
  104.     datalist = datalist.iloc[:, -1]
  105.     for i, item in enumerate(datalist):
  106.         if not item.startswith("pc"):
  107.             line = item
  108.         else:
  109.             datalist[i] = f"({line}) {item}"
  110.     datalist = [x for x in datalist
  111.                 for y in [".dat", ".out", ".txt"]
  112.                 if y in x]
  113.    
  114.     # Load proxy data and fill in missing values.
  115.     proxy = np.genfromtxt(filename_data)
  116.     fill_in_missing(proxy)    
  117.  
  118.     return proxy, datalist
  119.  
  120.  
  121. # Get MBH98 proxy data from Climategate archive.
  122. def get_climategate_proxy_data(t):
  123.     # Get Climategate archive.
  124.     url = "https://sealevel.info/FOIA/2009/FOIA/documents/mbh98-osborn/mbh98.tar"
  125.     filename = get_data_from(url)
  126.     directory = filename.split(".tar")[0]
  127.     untar(filename, directory)
  128.    
  129.     # Extended proxy list from Climategate dataset.
  130.     datalist = [
  131.         "ICE/MANNETAL97/quelc1-o18.dat",
  132.         "ICE/MANNETAL97/quelc1-accum.dat",
  133.         "ICE/MANNETAL97/quelc2-o18.dat",
  134.         "ICE/MANNETAL97/quelc2-accum.dat",
  135.         "ICE/MANNETAL97/svalbard.dat",
  136.         "ICE/FISHER/westgreen-o18.dat",
  137.         "TREE/MANNETAL97/tasmania.dat",
  138.         "TREE/MANNETAL97/npatagonia.dat",
  139.         "TREE/MANNETAL97/treeline11.dat",
  140.         "TREE/MANNETAL97/seprecip-nc.dat",
  141.         "TREE/MANNETAL97/seprecip-sc.dat",
  142.         "TREE/MANNETAL97/seprecip-ga.dat",
  143.         "TREE/MISC/yakut.dat",
  144.         "TREE/MISC/fenno-new.dat",
  145.         "TREE/MISC/urals-new.dat",
  146.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc01.out",
  147.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc02.out",
  148.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc03.out",
  149.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc04.out",
  150.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc05.out",
  151.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc06.out",
  152.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc07.out",
  153.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc08.out",
  154.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc09.out",
  155.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc10.out",
  156.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc11.out",
  157.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc12.out",
  158.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc13.out",
  159.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc14.out",
  160.         "TREE/ITRDB/NOAMER/BACKTO_1400/pc15.out",
  161.         "TREE/ITRDB/MISC/fran009.txt",
  162.         "TREE/ITRDB/MISC/fran010.txt",
  163.         "TREE/ITRDB/MISC/moro007.txt",
  164.         "TREE/ITRDB/MISC/moro008.txt",
  165.         "TREE/STAHLE/SWM/BACKTO_1400/pc01.out",
  166.         ]
  167.    
  168.     # Load proxy data and fill in missing values.
  169.     proxy = np.full((t.size, len(datalist) + 1), np.nan)
  170.     proxy[:, 0] = t
  171.     for i, filename in enumerate(datalist):
  172.         path = os.path.join(directory, filename)
  173.         p = pd.read_table(path, header=None)
  174.         p = p[0].str.split(expand=True).astype(float).to_numpy()
  175.         p = submatrix(p, t)
  176.         index = p[:, 0] - t[0]
  177.         proxy[index.astype(int), i + 1] = p[:, 1]
  178.     fill_in_missing(proxy)
  179.    
  180.     return proxy, datalist
  181.  
  182.  
  183. # Standardize data matrix x (time in first column) over interval t using MBH98
  184. # or WA standardization (detrended versus nondetrended standard deviation).
  185. def standardize(x, t, std):
  186.     x_z = x.copy()
  187.     x_t = submatrix(x, t)
  188.     mu = x_t[:, 1:].mean(axis=0)
  189.     x_z[:, 1:] -= mu
  190.     if std == "MBH":
  191.         sigma = sps.detrend(x_t[:, 1:], axis=0).std(axis=0, ddof=0)
  192.         x_z[:, 1:] /= sigma
  193.     if std == "WA":
  194.         sigma = x_t[:, 1:].std(axis=0, ddof=1)
  195.         x_z[:, 1:] /= sigma
  196.     return x_z
  197.  
  198.  
  199. # Regress the second column of y on all except the first column of x
  200. # over interval t.
  201. def regress(y, x, t):
  202.     y_t = submatrix(y, t)
  203.     x_t = submatrix(x, t)
  204.     beta_ols = np.linalg.lstsq(x_t[:, 1:], y_t[:, 1], rcond=None)[0]
  205.     y_ols = np.full((x.shape[0], 2), np.nan)
  206.     y_ols[:, 0] = x[:, 0]
  207.     y_ols[:, 1] = x[:, 1:] @ beta_ols
  208.     return beta_ols, y_ols
  209.  
  210.  
  211. # Make lineplot.
  212. def lineplot(x, y, fmt, ms, labels, title, figname):
  213.     fig, ax = plt.subplots(figsize=(8, 4.5), dpi=200)
  214.     for i in range(len(y)):
  215.         ax.plot(x, submatrix(y[i], x)[:, 1],
  216.                 fmt[i], ms=ms[i], label=labels[i])
  217.     ax.set_xlabel("Year", fontsize=10, labelpad=5)
  218.     ax.set_ylabel("Temperature anomaly ($\degree$C)", fontsize=10)
  219.     ax.set_title(title, loc="left", fontsize=12)
  220.     ax.tick_params(axis="both", labelsize=10)
  221.     ax.grid(b=True)
  222.     ax.legend(fontsize=10)
  223.     ax.set_facecolor((0.95, 0.95, 0.95))
  224.     fig.tight_layout(rect=(0, 0, 0.95, 1))
  225.     fig.savefig(figname)
  226.     plt.close(fig)
  227.  
  228.  
  229. # Make barplot.
  230. def barplot(beta, datalist, title, figname):
  231.     fig, ax = plt.subplots(figsize=(8, 10), dpi=200)
  232.     y = np.arange(len(datalist))
  233.     ax.barh(y, beta, height=0.8)
  234.     ax.axvline(color="black")
  235.     ax.invert_yaxis()
  236.     ax.set_yticks(y)
  237.     ax.set_yticklabels(datalist)
  238.     ax.set_ylim([len(y), -1])
  239.     ax.set_xlabel("Coefficient", fontsize=12, labelpad=5)
  240.     ax.set_title(title, fontsize=18, loc="left")
  241.     ax.tick_params(axis="both", labelsize=10)
  242.     ax.grid(b=True)
  243.     ax.set_axisbelow(b=True)
  244.     ax.set_facecolor((0.95, 0.95, 0.95))
  245.     fig.tight_layout()
  246.     fig.savefig(figname)
  247.     plt.close()
  248.  
  249.  
  250. # Make two barplots.
  251. def two_barplots(betas, datalists, titles, suptitle, figname):
  252.     fig, ax = plt.subplots(1, 2, figsize=(8, 5), dpi=200, sharex=True)
  253.     h = []
  254.     for i in range(2):
  255.         y = np.arange(len(datalists[i]))
  256.         h.append(ax[i].barh(y, betas[i], height=0.8))
  257.         ax[i].axvline(color="black")
  258.         ax[i].invert_yaxis()
  259.         ax[i].set_yticks(y)
  260.         ax[i].set_yticklabels(datalists[i])
  261.         ax[i].set_ylim([len(y), -1])
  262.         ax[i].set_xlabel("Coefficient", fontsize=8, labelpad=5)
  263.         ax[i].set_title(titles[i], fontsize=12, loc="left")
  264.         ax[i].tick_params(axis="both", labelsize=6)
  265.         ax[i].grid(b=True)
  266.         ax[i].set_axisbelow(b=True)
  267.         ax[i].set_facecolor((0.95, 0.95, 0.95))
  268.     ax[1].yaxis.tick_right()
  269.     for i, proxy in enumerate(zip(datalists[0], datalists[1])):
  270.         if not proxy[1].endswith(proxy[0].split("/")[-1]):
  271.             h[0][i].set_color("r")
  272.             h[1][i].set_color("r")
  273.     fig.tight_layout(rect=(0, 0, 1, 0.94))
  274.     fig.subplots_adjust(wspace=0.1)
  275.     x = np.mean([ax[0].get_position().xmax, ax[1].get_position().xmin])
  276.     fig.suptitle(suptitle, x=x, ha="center", fontsize=14)
  277.     fig.savefig(figname)
  278.     plt.close()
  279.    
  280.  
  281. # Perform the regression analysis.
  282. def main():
  283.     # Time intervals for reconstructions, calibrations and regressions.
  284.     t = np.arange(1400, 1981)
  285.     t_cal = np.arange(1902, 1981)
  286.     t_reg = np.arange(1400, 1450)
  287.        
  288.     # Get data.
  289.     mbh = get_mbh98_reconstruction()
  290.     wa = get_wahl_and_ammann_reconstruction()
  291.     proxy_wa, datalist_wa = get_wahl_and_ammann_proxy_data()
  292.     proxy_cg, datalist_cg = get_climategate_proxy_data(t)
  293.    
  294.     # MBH98 reconstruction regressed on WA proxy data.
  295.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  296.     _, mbh_reg_on_wa = regress(mbh, proxy_wa_zscores, t_reg)
  297.     x = t_reg
  298.     y = [mbh, mbh_reg_on_wa]
  299.     fmt = ["-"] * 2
  300.     ms = [None] * 2
  301.     labels = ["MBH98", "Regression"]
  302.     title = "MBH98 reconstruction regressed on WA proxy data"
  303.     figname = "fig1_mbh98_regressed_on_wa_data.png"
  304.     lineplot(x, y, fmt, ms, labels, title, figname)
  305.    
  306.     # MBH98 reconstruction regressed on WA proxy data,
  307.     # with WA reconstruction included for comparison.
  308.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  309.     _, mbh_reg_on_wa = regress(mbh, proxy_wa_zscores, t_reg)
  310.     x = t_reg
  311.     y = [mbh, mbh_reg_on_wa, wa]
  312.     fmt = ["-"] * 3
  313.     ms = [None] * 3
  314.     labels = ["MBH98", "MBH98 regressed on WA proxy data", "WA"]
  315.     title = "MBH98 and WA reconstructions"
  316.     figname = "fig2_mbh98_and_wa.png"
  317.     lineplot(x, y, fmt, ms, labels, title, figname)
  318.    
  319.     # WA reconstruction regressed on WA proxy data.
  320.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  321.     _, wa_reg_on_wa = regress(wa, proxy_wa_zscores, t_reg)
  322.     x = t_reg
  323.     y = [wa, wa_reg_on_wa]
  324.     fmt = ["-o", "-D"]
  325.     ms = [6, 3]
  326.     labels = ["WA", "Regression"]
  327.     title = "WA reconstruction regressed on WA proxy data"
  328.     figname = "fig3_wa_regressed_on_wa_data.png"
  329.     lineplot(x, y, fmt, ms, labels, title, figname)
  330.    
  331.     # MBH98 reconstruction regressed on Climategate proxy data.
  332.     proxy_cg_zscores = standardize(proxy_cg, t_cal, "MBH")
  333.     _, mbh_reg_on_cg = regress(mbh, proxy_cg_zscores, t_reg)
  334.     x = t_reg
  335.     y = [mbh, mbh_reg_on_cg]
  336.     fmt = ["-o", "-D"]
  337.     ms = [6, 3]
  338.     labels = ["MBH98", "Regression"]
  339.     title = "MBH98 reconstruction regressed on Climategate proxy data"
  340.     figname = "fig4_mbh98_regressed_on_climategate_data.png"
  341.     lineplot(x, y, fmt, ms, labels, title, figname)
  342.    
  343.     # Coefficients for WA regressed on WA proxy data.
  344.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  345.     beta_wa, _ = regress(wa, proxy_wa_zscores, t_reg)
  346.     title = "WA regression coefficients\nAD 1400–1449 interval"
  347.     figname = "fig5_coefficients_wa.png"
  348.     barplot(beta_wa, datalist_wa, title, figname)
  349.    
  350.     # Coefficients for MBH98 regressed on Climategate data.
  351.     proxy_cg_zscores = standardize(proxy_cg, t_cal, "MBH")
  352.     beta_mbh, _ = regress(mbh, proxy_cg_zscores, t_reg)
  353.     title = "MBH98 regression coefficients\nAD 1400–1449 interval"
  354.     figname = "fig6_coefficients_mbh98.png"
  355.     barplot(beta_mbh, datalist_cg, title, figname)
  356.    
  357.     # Regression coefficients for MBH98 and WA.
  358.     proxy_cg_zscores = standardize(proxy_cg, t_cal, "MBH")
  359.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  360.     beta_mbh, _ = regress(mbh, proxy_cg_zscores, t_reg)
  361.     beta_wa, _ = regress(wa, proxy_wa_zscores, t_reg)
  362.     zip_mbh = zip(datalist_cg, beta_mbh)
  363.     zip_mbh_nonzero = [x for x in zip_mbh if abs(x[1]) > 1e-6]
  364.     datalist_mbh, beta_mbh = zip(*zip_mbh_nonzero)
  365.     betas = [beta_mbh, beta_wa]
  366.     datalists = [datalist_mbh, datalist_wa]
  367.     titles = ["MBH98", "WA"]
  368.     suptitle = "AD 1400–1449 regression coefficients"
  369.     figname = "fig7_coefficients_mbh98_and_wa.png"
  370.     two_barplots(betas, datalists, titles, suptitle, figname)
  371.    
  372.     # Regression coefficients for MBH98 and WA using MBH98 standardization.
  373.     proxy_cg_zscores = standardize(proxy_cg, t_cal, "MBH")
  374.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "MBH")
  375.     beta_mbh, _ = regress(mbh, proxy_cg_zscores, t_reg)
  376.     beta_wa, _ = regress(wa, proxy_wa_zscores, t_reg)
  377.     zip_mbh = zip(datalist_cg, beta_mbh)
  378.     zip_mbh_nonzero = [x for x in zip_mbh if abs(x[1]) > 1e-6]
  379.     datalist_mbh, beta_mbh = zip(*zip_mbh_nonzero)
  380.     betas = [beta_mbh, beta_wa]
  381.     datalists = [datalist_mbh, datalist_wa]
  382.     titles = ["MBH98", "WA"]
  383.     suptitle = "AD 1400–1449 regression coefficients (MBH98 standardization)"
  384.     figname = "fig8_coefficients_mbh98_standardization.png"
  385.     two_barplots(betas, datalists, titles, suptitle, figname)
  386.    
  387.     # Regression coefficients for MBH98 and WA using WA standardization.
  388.     proxy_cg_zscores = standardize(proxy_cg, t_cal, "WA")
  389.     proxy_wa_zscores = standardize(proxy_wa, t_cal, "WA")
  390.     beta_mbh, _ = regress(mbh, proxy_cg_zscores, t_reg)
  391.     beta_wa, _ = regress(wa, proxy_wa_zscores, t_reg)
  392.     zip_mbh = zip(datalist_cg, beta_mbh)
  393.     zip_mbh_nonzero = [x for x in zip_mbh if abs(x[1]) > 1e-6]
  394.     datalist_mbh, beta_mbh = zip(*zip_mbh_nonzero)
  395.     betas = [beta_mbh, beta_wa]
  396.     datalists = [datalist_mbh, datalist_wa]
  397.     titles = ["MBH98", "WA"]
  398.     suptitle = "AD 1400–1449 regression coefficients (WA standardization)"
  399.     figname = "fig9_coefficients_wa_standardization.png"
  400.     two_barplots(betas, datalists, titles, suptitle, figname)
  401.  
  402.  
  403. if __name__ == '__main__':
  404.     main()
  405.  
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!