Figures_1-2-3-a.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45a361d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "This Source Code Form is subject to the terms of the Mozilla Public\n",
    "License, v. 2.0. If a copy of the MPL was not distributed with this\n",
    "file, You can obtain one at http://mozilla.org/MPL/2.0/.\n",
    "\n",
    "\n",
    "Authors:\n",
    "Roland Baatz roland.baatz@zalf.de\n",
    "\n",
    "Maintainers and contact:\n",
    "Currently maintained by the authors.\n",
    "\n",
    "This file requires output of the CRNS_*.ipynb\n",
    "\n",
    "Copyright (C) Leibniz Centre for Agricultural Landscape Research (ZALF)\n",
    "\"\"\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d75c1888",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "current_dir = os.getcwd()\n",
    "lib_dir = os.path.join(current_dir, 'lib')\n",
    "sys.path.append(lib_dir)\n",
    "\n",
    "import numpy as np\n",
    "import CRNS_library as CRNS_lib\n",
    "import time\n",
    "import scipy.io as sio\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import CRNS_optimizer as CRNS_opt\n",
    "from scipy.optimize import curve_fit\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "def func(x, p1, p2):\n",
    "    return p1 * np.cos(p2 * x) + p2 * np.sin(p1 * x)\n",
    "print('\\nBlock completed at time:                       '+ time.strftime(\"%H:%M:%S on %Y-%m-%d\")) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5060329c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fefe8fd2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95f2ecf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "### Iterate over EU Study sites\n",
    "#Read CRNS and calibrate\n",
    "plog_flag_id=0\n",
    "src_path = './data'\n",
    "suffix_csv = \"rmse_u\"\n",
    "csv_filename=src_path+\"my_dataframe_moderated_\"+suffix_csv+'.csv'\n",
    "\n",
    "\n",
    "if os.path.exists(csv_filename):\n",
    "    dfi = pd.read_csv(csv_filename, sep=';', index_col=0)  # Set index_col to 0 for row names\n",
    "    #LOAD DATAFRAME FROM CSV\n",
    "else:\n",
    "    print(\"ERROR File not present\")\n",
    "    exithere()\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5210b17",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4102bb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65772c97",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d87b180",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0790e65d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c0e1cb5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "269a2b62",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de703a44",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f99ba242",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f687b39",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d1f9171",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3d89afd",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "#df.drop(index=31, inplace=True) #ID: 31 Reading: LEC001\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(1, 3, figsize=(12, 4))\n",
    "\n",
    "ax=plt.subplot(1, 3, 1) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi_apty'],  df['psi_oulu'], color='k', marker='*', label='Apty vs Oulu')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "\n",
    "plt.xlabel('psi Apty (CR = 0.65 GV)')\n",
    "plt.ylabel('psi Oulu (CR = 0.8 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax=plt.subplot(1, 3,2) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi'],  df['psi_lmks'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Jungfraujoch (CR = 4.5 GV)')\n",
    "plt.ylabel('psi Lmks (CR = 3.8 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "ax=plt.subplot(1, 3,3) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi_mxco'],  df['psi_athens'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Mexico (CR = 8.3 GV)')\n",
    "plt.ylabel('psi Athens (CR = 9.5 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/3-subplots-psi-only{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8aebd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "#df.drop(index=31, inplace=True) #ID: 31 Reading: LEC001\n",
    "\n",
    "# Initialize subplots\n",
    "fig, axes = plt.subplots(2, 3, figsize=(12, 8))\n",
    "\n",
    "ax=plt.subplot(2, 3, 1) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi_apty'],  df['psi_oulu'], color='k', marker='*', label='Apty vs Oulu')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "\n",
    "plt.xlabel('psi Apty (CR = 0.65 GV)')\n",
    "plt.ylabel('psi Oulu (CR = 0.8 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax=plt.subplot(2, 3,2) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi'],  df['psi_lmks'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Jungfraujoch (CR = 4.5 GV)')\n",
    "plt.ylabel('psi Lmks (CR = 3.8 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "ax=plt.subplot(2, 3,3) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi_mxco'],  df['psi_athens'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Mexico (CR = 8.3 GV)')\n",
    "plt.ylabel('psi Athens (CR = 9.5 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "\n",
    "ax=plt.subplot(2, 3, 4) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi'],  df['psi_oulu'], color='k', marker='*', label='Apty vs Oulu')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "\n",
    "plt.xlabel('psi Jungfraujoch (CR = 4.5 GV)')\n",
    "plt.ylabel('psi  Oulu (CR = 0.65 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax=plt.subplot(2, 3,5) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi'],  df['psi_apty'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Jungfraujoch (CR = 4.5 GV)')\n",
    "plt.ylabel('psi Apti (CR = 0.65 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "ax=plt.subplot(2, 3,6) \n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "#fig=plt.figure(figsize=(3, 3))\n",
    "plt.scatter(df['psi'],  df['psi_athens'], color='k', marker='*', label='Apty vs Oulu')\n",
    "plt.xlabel('psi Jungfraujoch (CR = 4.5 GV)')\n",
    "plt.ylabel('psi Athens (CR = 9.5 GV)')\n",
    "plt.title('psi vs. Cutoff Rigidity')\n",
    "\"\"\"\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\"\"\"\n",
    "ax.plot([-0.5, 2], [-0.5,2], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/6-subplots-psi-only{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78053a79",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "634beeb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "df.drop(index=31, inplace=True)\n",
    "#ID: 31 Reading: LEC001\n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "\n",
    "# Create a new figure\n",
    "fig, axes = plt.subplots(figsize=(12, 5))\n",
    "# Plot error bars\n",
    "plt.subplot(1, 3, 1) \n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(df['CutoffRigidity'].iloc[ii], df['beta'].iloc[ii], yerr=df['beta_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "\n",
    "plt.scatter(df['CutoffRigidity'], df['beta'], color='k', marker='.', alpha=0.99)\n",
    "\n",
    "# Adding legend for blue stars (value) and grey lines (uncertainty)\n",
    "plt.plot([], [], 'b*', label='Value')  # Blue star legend\n",
    "plt.plot([], [], color='grey', linewidth=2, label='Uncertainty')  # Grey line legend\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['beta']\n",
    "\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 2)\n",
    "\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "\n",
    "curve_y = p(xaxis)\n",
    "plt.scatter(xdata, ydata, color='black', marker='*')\n",
    "plt.plot(xaxis, curve_y, color='black', label='Fitted line')\n",
    "\n",
    "plt.xlabel('Cutoff Rigidity [GV]')\n",
    "plt.ylabel('Beta (empirical)')\n",
    "plt.title('a)',loc='left')\n",
    "plt.grid(True)\n",
    "\n",
    "# Adding legend entries manually without plotting additional data\n",
    "plt.legend(loc=\"upper left\")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ax = plt.subplot(1, 3, 2) \n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(df['beta_McJD'].iloc[ii], df['beta'].iloc[ii], yerr=df['beta_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "# Set labels\n",
    "ax.set_xlabel('Beta (McJannet and Desilets (2023))')\n",
    "ax.set_ylabel('Beta (empirical)')\n",
    "\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "# Set the limits to be the same for both axes\n",
    "ax.set_xlim(min_val, max_val)\n",
    "ax.set_ylim(min_val, max_val)\n",
    "\n",
    "# Set the limits to be the same for both axes\n",
    "ax.set_xlim(x_min, x_max)\n",
    "ax.set_ylim(y_min, y_max)\n",
    "plt.xticks(rotation=45)\n",
    "\n",
    "# Add grid and legend\n",
    "ax.grid(True)\n",
    "ax.legend()\n",
    "ax.set_title('b)',loc='left')\n",
    "\n",
    "\n",
    "\n",
    "ax3 = plt.subplot(1, 3, 3) \n",
    "ax3.boxplot([df['beta_McJD'], df['beta']], vert=True, patch_artist=True)\n",
    "ax3.set_xticklabels(['Beta (McJannet and Desilets)', 'Beta (empirical)'])\n",
    "\n",
    "# Rotating the x-axis labels by 90 degrees\n",
    "ax3.set_xticklabels(['McJ. & D.', 'empiric'], rotation=45)\n",
    "\n",
    "# Annotate boxplots with axis labels\n",
    "ax3.set_ylabel('Beta values')\n",
    "ax3.set_title('c)', loc='left')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/Figure_1a-c_{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2efc3b4f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba750558",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "df.drop(index=31, inplace=True)\n",
    "#ID: 31 Reading: LEC001\n",
    "# Scatterplot: Beta vs. Cutoff Rigidity\n",
    "\n",
    "# Create a new figure\n",
    "fig, axes = plt.subplots(figsize=(8, 4))\n",
    "# Plot error bars\n",
    "plt.subplot(1, 2, 1) \n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(df['CutoffRigidity'].iloc[ii], df['beta'].iloc[ii], yerr=df['beta_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "\n",
    "plt.scatter(df['CutoffRigidity'], df['beta'], color='k', marker='.', alpha=0.7)\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['beta']\n",
    "\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 2)\n",
    "\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "\n",
    "curve_y = p(xaxis)\n",
    "plt.scatter(xdata, ydata, color='black', marker='*')\n",
    "plt.plot(xaxis, curve_y, color='black', label='Fitted Line')\n",
    "\n",
    "\n",
    "plt.xlabel('Cutoff Rigidity')\n",
    "plt.ylabel('Beta (empirical)')\n",
    "plt.title('a)',loc='left')\n",
    "plt.grid(True)\n",
    "\n",
    "\n",
    "ax = plt.subplot(1, 2, 2) \n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(df['beta_McJD'].iloc[ii], df['beta'].iloc[ii], yerr=df['beta_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "# Set labels\n",
    "ax.set_xlabel('Beta (McJannet and Desilets (2023))')\n",
    "ax.set_ylabel('Beta (empirical)')\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_max = ax.get_xlim()[1]\n",
    "y_max = ax.get_ylim()[1]\n",
    "max_val = max(x_max, y_max)\n",
    "\n",
    "# Get the maximum of the current axes\n",
    "x_min = ax.get_xlim()[0]\n",
    "y_min = ax.get_ylim()[0]\n",
    "min_val = max(x_min, y_min)\n",
    "\n",
    "min_val=np.min(df['beta'])\n",
    "max_val=np.max(df['beta'])\n",
    "\n",
    "# Plot the 1:1 line\n",
    "ax.plot([min_val, max_val], [min_val, max_val], color='red', linestyle='--', label='1-to-1 line')\n",
    "\n",
    "# Set the limits to be the same for both axes\n",
    "ax.set_xlim(min_val, max_val)\n",
    "ax.set_ylim(min_val, max_val)\n",
    "\n",
    "# Set the limits to be the same for both axes\n",
    "ax.set_xlim(x_min, x_max)\n",
    "ax.set_ylim(y_min, y_max)\n",
    "\n",
    "# Add grid and legend\n",
    "ax.grid(True)\n",
    "ax.legend()\n",
    "ax.set_title('b)',loc='left')\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/Figure_1a{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "df=dfi.copy()\n",
    "#df.drop(index=31, inplace=True) #ID: 31 Reading: LEC001s\n",
    "\n",
    "# Initialize subplots\n",
    "# Create a figure with a specific size\n",
    "fig = plt.figure(figsize=(8, 8))  # Adjust the size to make the plots approximately quadratic\n",
    "\n",
    "# Define a grid with 2 rows and 2 columns using GridSpec\n",
    "gs = gridspec.GridSpec(2, 2, height_ratios=[0.8, 1.2], width_ratios=[1, 1])\n",
    "\n",
    "# Create the subplots\n",
    "ax1 = fig.add_subplot(gs[0, 0])  # Top left subplot\n",
    "ax2 = fig.add_subplot(gs[0, 1])  # Top right subplot\n",
    "ax3 = fig.add_subplot(gs[1, :])  # Bottom subplot spanning both columns\n",
    "\n",
    "for ii in range(len(df)):\n",
    "    ax1.errorbar(df['CutoffRigidity'].iloc[ii], df['psi'].iloc[ii], yerr=df['psi_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "\n",
    "ax1.scatter(df['CutoffRigidity'], df['psi'], color='k', marker='.', alpha=0.7)\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['psi']\n",
    "\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "\n",
    "curve_y = p(xaxis)\n",
    "ax1.scatter(xdata, ydata, color='black', marker='*')\n",
    "ax1.plot(xaxis, curve_y, color='black', label='Fitted Line')\n",
    "uncertainty = np.sqrt(np.diag(pcov))\n",
    "#plt.fill_between(xaxis, curve_y - uncertainty, curve_y + uncertainty, color='grey', alpha=0.3)\n",
    "#plt.text(xaxis[10], curve_y[10], f'$f(x) = {popt[0]:.2f} \\\\cos({popt[1]:.2f}x) + {popt[1]:.2f} \\\\sin({popt[0]:.2f}x)$', fontsize=12)\n",
    "# Display the function as text on the plot\n",
    "#plt.text(xaxis[10], curve_y[10], f'$f(x) = {coefficients[0]:.2f} \\\\cos({coefficients[1]:.2f}x) + {coefficients[1]:.2f} \\\\sin({coefficients[0]:.2f}x)$', fontsize=12)\n",
    "function_text = f'$f(x) = {coefficients[0]:.2f}x^3 + {coefficients[1]:.2f}x^2 + {coefficients[2]:.2f}x + {coefficients[3]:.2f}$'\n",
    "\n",
    "ax1.set_xlabel('Cutoff Rigidity')\n",
    "ax1.set_ylabel('psi Jungfraujoch (empirical)')\n",
    "ax1.set_title('a)',loc='left')\n",
    "ax1.grid(True)\n",
    "\n",
    "ax2.scatter(df['psi_apty'],  df['psi_oulu'], color='k', marker='*', label='Apty vs Oulu')\n",
    "ax2.set_xlabel('psi Apty (CR=0.65 GV)')\n",
    "ax2.set_ylabel('psi Oulu (CR=0.8 GV)')\n",
    "ax2.set_title('b)',loc='left')\n",
    "ax2.grid(True)\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['psi_apty']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='grey', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='grey', label='Apty 0.7')\n",
    "\n",
    "\n",
    "ydata = df['psi_oulu']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='blue', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='blue', label='Oulu 0.8', linestyle='--')\n",
    "\n",
    "ydata = df['psi_lmks']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='m', marker='.')\n",
    "ax3.plot(xaxis, curve_y, color='m', label='Lmks 3.8')\n",
    "\n",
    "\n",
    "ydata = df['psi']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='black', marker='*')\n",
    "ax3.plot(xaxis, curve_y, color='black', label='Jungfraujoch 4.5', linestyle='--')\n",
    "\n",
    "ydata = df['psi_mxco']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='red', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='red', label='Mexico 8.3')\n",
    "\n",
    "\n",
    "\n",
    "ydata = df['psi_athens']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='y', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='y', label='Athens 8.5', linestyle='--')\n",
    "\n",
    "\n",
    "\n",
    "ax3.set_xlabel('Cutoff Rigidity')\n",
    "ax3.set_ylabel('psi Jungfraujoch (empirical)')\n",
    "ax3.set_title('c)',loc='left')\n",
    "ax3.grid(True)\n",
    "ax3.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/Figure_2_psi{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "df=dfi.copy()\n",
    "#df.drop(index=31, inplace=True)\n",
    "\n",
    "fig = plt.figure(figsize=(4,4))\n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(ii, df['omega'].iloc[ii], yerr=df['omega_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "# Adding legend for blue stars (value) and grey lines (uncertainty)\n",
    "plt.plot([], [], 'b*', label='Value')  # Blue star legend\n",
    "plt.plot([], [], color='grey', linewidth=2, label='Uncertainty')  # Grey line legend\n",
    "\n",
    "# Adding legend entries manually without plotting additional data\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.xlabel('Site ID')\n",
    "plt.ylabel('Omega')\n",
    "# Adjust layout to not overlap plots\n",
    "plt.tight_layout()\n",
    "pngname=f\"./plots/Figure_3_omegaa{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "\n",
    "\n",
    "df=dfi.copy()\n",
    "df.drop(index=59, inplace=True)\n",
    "\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "for ii in range(len(df)):\n",
    "    plt.errorbar(ii, df['omega'].iloc[ii], yerr=df['omega_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "# Adding legend for blue stars (value) and grey lines (uncertainty)\n",
    "plt.plot([], [], 'b*', label='Value')  # Blue star legend\n",
    "plt.plot([], [], color='grey', linewidth=2, label='Uncertainty')  # Grey line legend\n",
    "\n",
    "# Adding legend entries manually without plotting additional data\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.xlabel('Site ID')\n",
    "plt.ylabel('Omega')\n",
    "# Adjust layout to not overlap plots\n",
    "plt.tight_layout()\n",
    "pngname=f\"./plots/Figure_3_omega_{suffix_csv}_outlier-rm.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "plt.close(\"all\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef3c1db",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "from matplotlib.lines import Line2D\n",
    "\n",
    "\n",
    "\n",
    "df=dfi.copy()\n",
    "#df.drop(index=31, inplace=True) #ID: 31 Reading: LEC001\n",
    "\n",
    "# Initialize subplots\n",
    "# Create a figure with a specific size\n",
    "fig = plt.figure(figsize=(8, 8))  # Adjust the size to make the plots approximately quadratic\n",
    "\n",
    "# Define a grid with 2 rows and 2 columns using GridSpec\n",
    "gs = gridspec.GridSpec(2, 2, height_ratios=[0.8, 1.2], width_ratios=[1, 1])\n",
    "\n",
    "# Create the subplots\n",
    "ax1 = fig.add_subplot(gs[0, 0])  # Top left subplot\n",
    "ax2 = fig.add_subplot(gs[0, 1])  # Top right subplot\n",
    "ax3 = fig.add_subplot(gs[1, :])  # Bottom subplot spanning both columns\n",
    "\n",
    "for ii in range(len(df)):\n",
    "    ax1.errorbar(df['CutoffRigidity'].iloc[ii], df['psi'].iloc[ii], yerr=df['psi_sd'].iloc[ii],\n",
    "                 fmt='*', capsize=0, color='blue', ecolor='grey', elinewidth=2, capthick=2)\n",
    "\n",
    "ax1.scatter(df['CutoffRigidity'], df['psi'], color='k', marker='.', alpha=0.7)\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['psi']\n",
    "\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "\n",
    "curve_y = p(xaxis)\n",
    "ax1.scatter(xdata, ydata, color='black', marker='*')\n",
    "ax1.plot(xaxis, curve_y, color='black', label='Fitted Line')\n",
    "uncertainty = np.sqrt(np.diag(pcov))\n",
    "#plt.fill_between(xaxis, curve_y - uncertainty, curve_y + uncertainty, color='grey', alpha=0.3)\n",
    "#plt.text(xaxis[10], curve_y[10], f'$f(x) = {popt[0]:.2f} \\\\cos({popt[1]:.2f}x) + {popt[1]:.2f} \\\\sin({popt[0]:.2f}x)$', fontsize=12)\n",
    "# Display the function as text on the plot\n",
    "#plt.text(xaxis[10], curve_y[10], f'$f(x) = {coefficients[0]:.2f} \\\\cos({coefficients[1]:.2f}x) + {coefficients[1]:.2f} \\\\sin({coefficients[0]:.2f}x)$', fontsize=12)\n",
    "function_text = f'$f(x) = {coefficients[0]:.2f}x^3 + {coefficients[1]:.2f}x^2 + {coefficients[2]:.2f}x + {coefficients[3]:.2f}$'\n",
    "\n",
    "ax1.set_xlabel('Cutoff Rigidity')\n",
    "ax1.set_ylabel('psi Jung (empirical)')\n",
    "ax1.set_title('a)',loc='left')\n",
    "ax1.grid(True)\n",
    "\n",
    "ax2.scatter(df['psi_apty'],  df['psi_oulu'], color='k', marker='*', label='Apty vs Oulu')\n",
    "ax2.set_xlabel('psi Apty (CR=0.65 GV)')\n",
    "ax2.set_ylabel('psi Oulu (CR=0.8 GV)')\n",
    "ax2.set_title('b)',loc='left')\n",
    "ax2.grid(True)\n",
    "\n",
    "xdata = df['CutoffRigidity']\n",
    "ydata = df['psi_apty']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='grey', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='grey', label='Apty 0.7')\n",
    "\n",
    "\n",
    "ydata = df['psi_oulu']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='blue', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='blue', label='Oulu 0.8', linestyle='--')\n",
    "\n",
    "ydata = df['psi_lmks']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='m', marker='.')\n",
    "ax3.plot(xaxis, curve_y, color='m', label='Lmks 3.8')\n",
    "\n",
    "\n",
    "ydata = df['psi']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='black', marker='*')\n",
    "ax3.plot(xaxis, curve_y, color='black', label='Jungfraujoch 4.5', linestyle='--')\n",
    "\n",
    "\n",
    "##Hawdon et al. 2014 paper\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = -0.075*(xaxis-4.5)+1\n",
    "ax3.plot(xaxis, curve_y, color='black', label='Hawdon et al. 2014', linestyle=':')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "ydata = df['psi_mxco']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='red', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='red', label='Mexico 8.3')\n",
    "\n",
    "\n",
    "\n",
    "ydata = df['psi_athens']\n",
    "popt, pcov = curve_fit(func, xdata, ydata, p0=(1.0, 0.3))\n",
    "coefficients = np.polyfit(xdata, ydata, 3)\n",
    "# Create a polynomial function\n",
    "p = np.poly1d(coefficients)\n",
    "# Now you can use 'p' to evaluate the polynomial at any point 'x':\n",
    "xaxis = np.linspace(min(xdata), max(xdata), 100)\n",
    "curve_y = func(xaxis, *popt)\n",
    "curve_y = p(xaxis)\n",
    "ax3.scatter(xdata, ydata, color='y', marker='x')\n",
    "ax3.plot(xaxis, curve_y, color='y', label='Athens 8.5', linestyle='--')\n",
    "\n",
    "\n",
    "\n",
    "# Create custom legend entries for each line\n",
    "custom_lines = [\n",
    "    Line2D([0], [0], color='grey', linestyle='--', marker='x', label='Apty 0.7'),\n",
    "    Line2D([0], [0], color='blue', linestyle='--', marker='x', label='Oulu 0.8'),\n",
    "    Line2D([0], [0], color='m', marker='.', label='Lmks 3.8'),\n",
    "    Line2D([0], [0], color='black', linestyle='--', marker='*', label='Jungfraujoch 4.5'),\n",
    "    Line2D([0], [0], color='black', linestyle=':', label='Hawdon et al. 2014'),\n",
    "    Line2D([0], [0], color='red', marker='x', label='Mexico 8.3'),\n",
    "    Line2D([0], [0], color='y', linestyle='--', marker='x', label='Athens 8.5')\n",
    "]\n",
    "\n",
    "# Display the custom legend\n",
    "ax3.legend(handles=custom_lines)\n",
    "\n",
    "ax3.set_xlabel('Cutoff Rigidity')\n",
    "ax3.set_ylabel('psi Jung (empirical)')\n",
    "ax3.set_title('c)',loc='left')\n",
    "ax3.grid(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "pngname=f\"./plots/Figure_a2_psi_Hawdon__a{suffix_csv}.png\"\n",
    "print(\"Save file as \"+pngname)\n",
    "fig.savefig(pngname, bbox_inches='tight')\n",
    "plt.close(\"all\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f769b16",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"MEANS\")\n",
    "print(np.mean(dfi['psi']))\n",
    "print(np.mean(dfi['psi_sd']))\n",
    "print(\"Oulu:\")\n",
    "print(np.mean(dfi['psi_oulu']))\n",
    "print(np.mean(dfi['psi_apty']))\n",
    "print(\"  Athens and Mxco and lmks:\")\n",
    "print(np.mean(dfi['psi_athens']))\n",
    "print(np.mean(dfi['psi_mxco']))\n",
    "print(np.mean(dfi['psi_lmks']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7495f417",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"JUNG\")\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "correlations = {\n",
    "    'psi_oulu': np.corrcoef(dfi['psi'], dfi['psi_oulu'])[0, 1],\n",
    "    'psi_apty': np.corrcoef(dfi['psi'], dfi['psi_apty'])[0, 1],\n",
    "    'psi_athens': np.corrcoef(dfi['psi'], dfi['psi_athens'])[0, 1],\n",
    "    'psi_mxco': np.corrcoef(dfi['psi'], dfi['psi_mxco'])[0, 1],\n",
    "    'psi_lmks': np.corrcoef(dfi['psi'], dfi['psi_lmks'])[0, 1]\n",
    "}\n",
    "\n",
    "correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4a0a0a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"psi_oulu\")\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "correlations = {\n",
    "    'psi_oulu': np.corrcoef(dfi['psi_oulu'], dfi['psi_oulu'])[0, 1],\n",
    "    'psi_apty': np.corrcoef(dfi['psi_oulu'], dfi['psi_apty'])[0, 1],\n",
    "    'psi_athens': np.corrcoef(dfi['psi_oulu'], dfi['psi_athens'])[0, 1],\n",
    "    'psi_mxco': np.corrcoef(dfi['psi_oulu'], dfi['psi_mxco'])[0, 1],\n",
    "    'psi_lmks': np.corrcoef(dfi['psi_oulu'], dfi['psi_lmks'])[0, 1]\n",
    "}\n",
    "\n",
    "correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6fade8",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ATHENS\")\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "correlations = {\n",
    "    'psi JUNG': np.corrcoef(dfi['psi_athens'], dfi['psi'])[0, 1],\n",
    "    'psi_oulu': np.corrcoef(dfi['psi_athens'], dfi['psi_oulu'])[0, 1],\n",
    "    'psi_apty': np.corrcoef(dfi['psi_athens'], dfi['psi_apty'])[0, 1],\n",
    "    'psi_athens': np.corrcoef(dfi['psi_athens'], dfi['psi_athens'])[0, 1],\n",
    "    'psi_mxco': np.corrcoef(dfi['psi_athens'], dfi['psi_mxco'])[0, 1],\n",
    "    'psi_lmks': np.corrcoef(dfi['psi_athens'], dfi['psi_lmks'])[0, 1]\n",
    "}\n",
    "\n",
    "correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fffcb583",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Omega\")\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "print(np.min(dfi['omega']), \"max:\", np.max(dfi['omega']))\n",
    "print(np.std(dfi['omega']), \"=Std dev\")\n",
    "print(np.mean(dfi['omega']), \"=mean and median=\",np.median(dfi['omega']))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0f65ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dfi['omega_sd'])\n",
    "# Remove the 3 largest and smallest omega values\n",
    "sorted_omega = np.sort(dfi['omega'])\n",
    "trimmed_omega = sorted_omega[3:-3]\n",
    "\n",
    "# Print statistics after removing 3 largest and 3 smallest values\n",
    "print(\"\\nAfter removing 3 largest and smallest estimates:\")\n",
    "print(\"Min:\", np.min(trimmed_omega), \"Max:\", np.max(trimmed_omega))\n",
    "print(\"Std dev:\", np.std(trimmed_omega))\n",
    "print(\"Mean:\", np.mean(trimmed_omega), \"Median:\", np.median(trimmed_omega))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1db808b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "df.drop(index=31, inplace=True)\n",
    "\n",
    "myvar = \"beta\"\n",
    "print(\"Analysis for:\", myvar)\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "print(np.min(df[myvar]), \"max:\", np.max(df[myvar]))\n",
    "print(np.std(df[myvar]), \"=Std dev\")\n",
    "print(np.mean(df[myvar]), \"=mean and median=\",np.median(df[myvar]))\n",
    "\n",
    "# Filtering the dataset by CutoffRigidity <=4 and >=4\n",
    "df_filtered = df[(df['CutoffRigidity'] <= 4)]# & (df['CutoffRigidity'] >= 4)]\n",
    "\n",
    "# Calculate and print statistics for the filtered data\n",
    "print(\"MIN= \",np.nanmin(df_filtered[myvar]))\n",
    "print(\"max= \",np.nanmax(df_filtered[myvar]))\n",
    "print(\"std= \", np.nanstd(df_filtered[myvar]))\n",
    "print(\"mean= \", np.nanmean(df_filtered[myvar]))\n",
    "print(\"median= \", np.nanmedian(df_filtered[myvar]))\n",
    "\n",
    "# Filtering the dataset by CutoffRigidity <=4 and >=4\n",
    "df_filtered = df[(df['CutoffRigidity'] >= 6)]# & (df['CutoffRigidity'] >= 4)]\n",
    "\n",
    "# Calculate and print statistics for the filtered data\n",
    "print(\"MIN= \",np.nanmin(df_filtered[myvar]))\n",
    "print(\"max= \",np.nanmax(df_filtered[myvar]))\n",
    "print(\"std= \", np.nanstd(df_filtered[myvar]))\n",
    "print(\"mean= \", np.nanmean(df_filtered[myvar]))\n",
    "print(\"median= \", np.nanmedian(df_filtered[myvar]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee1eede8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df=dfi.copy()\n",
    "df.drop(index=31, inplace=True)\n",
    "\n",
    "print(\"Analysis for:\", myvar)\n",
    "# Calculate Pearson correlation coefficients between 'psi' and other columns\n",
    "print(np.min(df[myvar]), \"max:\", np.max(df[myvar]))\n",
    "print(np.std(df[myvar]), \"=Std dev\")\n",
    "print(np.mean(df[myvar]), \"=mean and median=\",np.median(df[myvar]))\n",
    "print(len(df))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11094861",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9 (py_env_lda)",
   "language": "python",
   "name": "py_env_lda"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}