{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fe94572f",
   "metadata": {},
   "source": [
    "# Volcanic Forcing Statistical Analysis - Figure 1\n",
    "\n",
    "This notebook performs statistical analysis of volcanic forcing records to characterize the distribution of volcanic eruption magnitudes. The analysis combines data from two key sources and fits Generalized Extreme Value (GEV) distributions to understand volcanic forcing statistics.\n",
    "\n",
    "## Data Sources\n",
    "\n",
    "**1. Lin et al. MIS3 Data (Marine Isotope Stage 3, ~60-24 kyr BP):**\n",
    "- Ice core-based volcanic sulfate deposition record\n",
    "- Covers the last glacial period when Dansgaard-Oeschger (D-O) events occurred\n",
    "- Provides eruption magnitude (stratospheric sulfate loading in Tg SO₄)\n",
    "- Source location classified as: Northern Hemisphere (NH), Equatorial (EQ), Southern Hemisphere (SH)\n",
    "\n",
    "**2. Sigl et al. 2022 Holocene Volcanic Dataset (HolVol):**\n",
    "- Comprehensive volcanic forcing reconstruction for the Holocene (~12 kyr BP to present)\n",
    "- Volcanic Stratospheric Sulfur Injection (VSSI) estimates\n",
    "- Georeferenced eruption locations (latitude/longitude)\n",
    "\n",
    "## Analysis Workflow\n",
    "\n",
    "1. **Data Loading and Preparation:**\n",
    "   - Load Lin et al. MIS3 volcanic record from Excel spreadsheet\n",
    "   - Load Sigl et al. 2022 HolVol dataset from NetCDF\n",
    "   - Load NGRIP temperature reconstruction (Pedro et al. 2018) for climate context\n",
    "\n",
    "2. **Data Binning and Resampling:**\n",
    "   - Divide MIS3 time period into equal-duration bins\n",
    "   - Extract maximum volcanic load in each bin\n",
    "   - Creates evenly-sampled dataset suitable for statistical analysis\n",
    "\n",
    "3. **Statistical Distribution Fitting:**\n",
    "   - Fit Generalized Extreme Value (GEV) distribution to volcanic loading data\n",
    "   - GEV distribution characterizes extreme events (large eruptions)\n",
    "   - Extract shape, location, and scale parameters for both datasets\n",
    "\n",
    "4. **Visualization:**\n",
    "   - Time series of volcanic events by hemisphere\n",
    "   - Histograms of volcanic loading magnitudes\n",
    "   - GEV probability density functions overlaid on histograms\n",
    "\n",
    "## Key Outputs\n",
    "\n",
    "- **GEV Parameters:** Shape (ξ), location (μ), and scale (σ) parameters characterizing volcanic forcing statistics\n",
    "- **Figure 1 Components:** Visual comparison of Lin (MIS3) and Sigl (Holocene) volcanic forcing distributions\n",
    "- **Statistical Foundation:** These parameters are used to generate synthetic volcanic forcing scenarios for climate model experiments\n",
    "\n",
    "## Scientific Context\n",
    "\n",
    "Understanding the statistical distribution of volcanic eruption magnitudes is essential for:\n",
    "- Generating realistic synthetic volcanic forcing scenarios\n",
    "- Assessing probability of extreme volcanic events\n",
    "- Comparing glacial (MIS3) vs. interglacial (Holocene) volcanic activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "41020476-7a65-4be6-8073-9539f9adc452",
   "metadata": {},
   "outputs": [],
   "source": [
    "# numpy/scipy stuff\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "\n",
    "from scipy.signal import argrelextrema\n",
    "from scipy.stats import genextreme as gev\n",
    "\n",
    "from netCDF4 import Dataset \n",
    "\n",
    "# matplotlib stuff\n",
    "import matplotlib.pyplot as plt\n",
    "#plt.style.use('seaborn-talk')\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib import colors\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib.patches import FancyBboxPatch\n",
    "import matplotlib as mpl\n",
    "import matplotlib.lines as mlines\n",
    "import matplotlib.colors as mcolors\n",
    "\n",
    "# data sci stuff\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "\n",
    "# system stuff\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "import fnmatch\n",
    "\n",
    "# netcdf stuff\n",
    "from netCDF4 import Dataset\n",
    "import datetime\n",
    "import cftime\n",
    "import nc_time_axis\n",
    "\n",
    "# widgets stuff\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact,interactive\n",
    "import xlrd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5ba5c282-456d-4e82-9bb7-f92f66b42652",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "numpy version:  2.2.6\n",
      "matplotlib version:  3.10.3\n",
      "xarray version:  2025.11.0\n",
      "pandas version:  2.2.3\n"
     ]
    }
   ],
   "source": [
    "print(\"numpy version: \", np.__version__)\n",
    "print(\"matplotlib version: \", mpl.__version__)\n",
    "print(\"xarray version: \", xr.__version__)\n",
    "print(\"pandas version: \", pd.__version__)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "227c23bc-1f85-4a62-b89b-c0615e7a5831",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setaxesfsize(axl,fontsize, xlabel, ylabel, labsize):\n",
    "    ticklabelbot = axl.get_xticklabels()\n",
    "    ticklabelleft = axl.get_yticklabels()\n",
    "    for labelx in ticklabelbot:\n",
    "        labelx.set_fontsize(fontsize)\n",
    "#    i=0\n",
    "    for labely in ticklabelleft:\n",
    "#        if i != 0:\n",
    "            labely.set_fontsize(fontsize)\n",
    "    \n",
    "    axl.set_xlabel(xlabel, size=labsize)\n",
    "    axl.set_ylabel(ylabel, size=labsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6ac67cc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Year</th>\n",
       "      <th>Temperature</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20000</td>\n",
       "      <td>-45.005556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>20020</td>\n",
       "      <td>-44.984211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>20040</td>\n",
       "      <td>-44.916000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>20060</td>\n",
       "      <td>-44.812000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>20080</td>\n",
       "      <td>-44.678000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2246</th>\n",
       "      <td>64920</td>\n",
       "      <td>-48.824286</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2247</th>\n",
       "      <td>64940</td>\n",
       "      <td>-48.369000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2248</th>\n",
       "      <td>64960</td>\n",
       "      <td>-48.107000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2249</th>\n",
       "      <td>64980</td>\n",
       "      <td>-48.131053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2250</th>\n",
       "      <td>65000</td>\n",
       "      <td>-48.407500</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2251 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Year  Temperature\n",
       "0     20000   -45.005556\n",
       "1     20020   -44.984211\n",
       "2     20040   -44.916000\n",
       "3     20060   -44.812000\n",
       "4     20080   -44.678000\n",
       "...     ...          ...\n",
       "2246  64920   -48.824286\n",
       "2247  64940   -48.369000\n",
       "2248  64960   -48.107000\n",
       "2249  64980   -48.131053\n",
       "2250  65000   -48.407500\n",
       "\n",
       "[2251 rows x 2 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ngrip_t_pedro_2018 = pd.read_csv(\"/scratch/Dropbox/ipynb/D-O_analysis/volcanic/Pedro_QSR2018_fig1a.txt\", sep='\\t', names=[\"Year\", \"Temperature\"])\n",
    "\n",
    "ngrip_t_pedro_2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a40113cf-9d0c-4f03-be72-c80343106a57",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigl2022 = xr.open_dataset(\"/scratch/Dropbox/ipynb/D-O_analysis/volcanic/HolVol_volcanic_stratospheric_sulfur_injection_v1.0.nc\", engine='netcdf4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9047fb0c-7e6f-49ce-856d-26bd40b94da1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in notebooks */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base rgba(0, 0, 0, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, white)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 15))\n",
       "  );\n",
       "}\n",
       "\n",
       "html[theme=\"dark\"],\n",
       "html[data-theme=\"dark\"],\n",
       "body[data-theme=\"dark\"],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, #111111)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 15))\n",
       "  );\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "  line-height: 1.6;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-obj-name,\n",
       ".xr-group-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-group-name::before {\n",
       "  content: \"📁\";\n",
       "  padding-right: 0.3em;\n",
       "}\n",
       "\n",
       ".xr-group-name,\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 0 20px 0 20px;\n",
       "  margin-block-start: 0;\n",
       "  margin-block-end: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: inline-block;\n",
       "  opacity: 0;\n",
       "  height: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "  border: 2px solid transparent !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:focus + label {\n",
       "  border: 2px solid var(--xr-font-color0) !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: \"►\";\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: \"▼\";\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-top: 4px;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-group-box {\n",
       "  display: inline-grid;\n",
       "  grid-template-columns: 0px 20px auto;\n",
       "  width: 100%;\n",
       "}\n",
       "\n",
       ".xr-group-box-vline {\n",
       "  grid-column-start: 1;\n",
       "  border-right: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "  width: 0px;\n",
       "}\n",
       "\n",
       ".xr-group-box-hline {\n",
       "  grid-column-start: 2;\n",
       "  grid-row-start: 1;\n",
       "  height: 1em;\n",
       "  width: 20px;\n",
       "  border-bottom: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-group-box-contents {\n",
       "  grid-column-start: 3;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: \"(\";\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: \")\";\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: \",\";\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  border-color: var(--xr-background-color-row-odd);\n",
       "  margin-bottom: 0;\n",
       "  padding-top: 2px;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "  border-color: var(--xr-background-color-row-even);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  border-top: 2px dotted var(--xr-background-color);\n",
       "  padding-bottom: 20px !important;\n",
       "  padding-top: 10px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in + label,\n",
       ".xr-var-data-in + label,\n",
       ".xr-index-data-in + label {\n",
       "  padding: 0 1px;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-data > pre,\n",
       ".xr-index-data > pre,\n",
       ".xr-var-data > table > tbody > tr {\n",
       "  background-color: transparent !important;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked + label > .xr-icon-file-text2,\n",
       ".xr-var-data-in:checked + label > .xr-icon-database,\n",
       ".xr-index-data-in:checked + label > .xr-icon-database {\n",
       "  color: var(--xr-font-color0);\n",
       "  filter: drop-shadow(1px 1px 5px var(--xr-font-color2));\n",
       "  stroke-width: 0.8px;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt; Size: 100kB\n",
       "Dimensions:     (nerup: 1384)\n",
       "Dimensions without coordinates: nerup\n",
       "Data variables:\n",
       "    year        (nerup) float64 11kB ...\n",
       "    yearCE      (nerup) float64 11kB ...\n",
       "    month       (nerup) float64 11kB ...\n",
       "    day         (nerup) float64 11kB ...\n",
       "    lat         (nerup) float64 11kB ...\n",
       "    vssi        (nerup) float64 11kB ...\n",
       "    hemi        (nerup) float64 11kB ...\n",
       "    sigma_vssi  (nerup) float64 11kB ...\n",
       "    duration    (nerup) float64 11kB ...\n",
       "Attributes:\n",
       "    title:    Ice core-inferred volcanic stratospheric sulfur injection for t...\n",
       "    history:  Created Tue Jun  9 23:01:18 2020\n",
       "    source:   Sigl et al., in preparation</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-a448d900-eb0f-482e-ba94-9f44d21f354b' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-a448d900-eb0f-482e-ba94-9f44d21f354b' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span>nerup</span>: 1384</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-81b04172-db78-407a-b18c-55c0eb11d533' class='xr-section-summary-in' type='checkbox'  checked><label for='section-81b04172-db78-407a-b18c-55c0eb11d533' class='xr-section-summary' >Data variables: <span>(9)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>year</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-594e102a-ca4c-43ca-98cf-9e609f4bbf80' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-594e102a-ca4c-43ca-98cf-9e609f4bbf80' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-acd61213-79fd-4dc3-8327-bc08e7503046' class='xr-var-data-in' type='checkbox'><label for='data-acd61213-79fd-4dc3-8327-bc08e7503046' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>year (ISO)</dd><dt><span>long_name :</span></dt><dd>year of eruption (ISO 8601, including year 0)</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>yearCE</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-537d48de-6245-4435-8fa8-fbfa1edae296' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-537d48de-6245-4435-8fa8-fbfa1edae296' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-2ed6882b-4ea1-4d4e-bf8f-652c0f9d7d37' class='xr-var-data-in' type='checkbox'><label for='data-2ed6882b-4ea1-4d4e-bf8f-652c0f9d7d37' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>year (BCE/CE)</dd><dt><span>long_name :</span></dt><dd>year of eruption (BCE/CE, no year 0)</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>month</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-24952ab6-ea85-4009-b5e4-e5b41e18cb6c' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-24952ab6-ea85-4009-b5e4-e5b41e18cb6c' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-e911312b-b631-46ed-90cb-4562d3bbcead' class='xr-var-data-in' type='checkbox'><label for='data-e911312b-b631-46ed-90cb-4562d3bbcead' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>month</dd><dt><span>long_name :</span></dt><dd>month of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>day</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-81f6fff6-e067-4a8c-a35b-779608949816' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-81f6fff6-e067-4a8c-a35b-779608949816' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-bd295dee-ff61-45e9-9d4f-7f6354d9d0e8' class='xr-var-data-in' type='checkbox'><label for='data-bd295dee-ff61-45e9-9d4f-7f6354d9d0e8' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>day</dd><dt><span>long_name :</span></dt><dd>day of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>lat</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-2b3930f9-57f9-469e-8232-781e5fbc3164' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-2b3930f9-57f9-469e-8232-781e5fbc3164' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-565f982f-c2f4-4ea1-af20-92e2a8e97ae0' class='xr-var-data-in' type='checkbox'><label for='data-565f982f-c2f4-4ea1-af20-92e2a8e97ae0' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>long_name :</span></dt><dd>latitude of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>vssi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-202f2f52-1be9-41da-8f07-a312c0ad9b85' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-202f2f52-1be9-41da-8f07-a312c0ad9b85' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-506501bf-65e8-4af0-a9e0-b9891bef67dd' class='xr-var-data-in' type='checkbox'><label for='data-506501bf-65e8-4af0-a9e0-b9891bef67dd' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>TgS</dd><dt><span>long_name :</span></dt><dd>volcanic stratospheric sulfur injection of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>hemi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-410ce028-937a-4414-969a-4caddc45cd38' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-410ce028-937a-4414-969a-4caddc45cd38' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-d8fffac2-131f-436d-b5fb-0ef704fd1548' class='xr-var-data-in' type='checkbox'><label for='data-d8fffac2-131f-436d-b5fb-0ef704fd1548' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd></dd><dt><span>long_name :</span></dt><dd>hemispheric asymmetry (NH/SH) of aerosol spread for tropical eruptions based on ratio of Greenland to Antarctic flux</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>sigma_vssi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-b2c02cd2-cb74-4215-afb3-8cc4398cadb0' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-b2c02cd2-cb74-4215-afb3-8cc4398cadb0' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-00e522ec-de34-4387-bf14-302550ed1485' class='xr-var-data-in' type='checkbox'><label for='data-00e522ec-de34-4387-bf14-302550ed1485' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>TgS</dd><dt><span>long_name :</span></dt><dd>Random uncertainty (1 sigma) in VSSI</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>duration</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-f864c63f-d80b-4b43-8233-50988a87d0d4' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-f864c63f-d80b-4b43-8233-50988a87d0d4' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-fce08c2a-f64c-4e4e-9092-aa1dee9742ab' class='xr-var-data-in' type='checkbox'><label for='data-fce08c2a-f64c-4e4e-9092-aa1dee9742ab' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>months</dd><dt><span>long_name :</span></dt><dd>duration of effusive eruptions in months</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-24d07bca-43fa-4e61-8153-cfe4f3fa696b' class='xr-section-summary-in' type='checkbox'  checked><label for='section-24d07bca-43fa-4e61-8153-cfe4f3fa696b' class='xr-section-summary' >Attributes: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>title :</span></dt><dd>Ice core-inferred volcanic stratospheric sulfur injection for the Holocene</dd><dt><span>history :</span></dt><dd>Created Tue Jun  9 23:01:18 2020</dd><dt><span>source :</span></dt><dd>Sigl et al., in preparation</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset> Size: 100kB\n",
       "Dimensions:     (nerup: 1384)\n",
       "Dimensions without coordinates: nerup\n",
       "Data variables:\n",
       "    year        (nerup) float64 11kB ...\n",
       "    yearCE      (nerup) float64 11kB ...\n",
       "    month       (nerup) float64 11kB ...\n",
       "    day         (nerup) float64 11kB ...\n",
       "    lat         (nerup) float64 11kB ...\n",
       "    vssi        (nerup) float64 11kB ...\n",
       "    hemi        (nerup) float64 11kB ...\n",
       "    sigma_vssi  (nerup) float64 11kB ...\n",
       "    duration    (nerup) float64 11kB ...\n",
       "Attributes:\n",
       "    title:    Ice core-inferred volcanic stratospheric sulfur injection for t...\n",
       "    history:  Created Tue Jun  9 23:01:18 2020\n",
       "    source:   Sigl et al., in preparation"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4b935c21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the Jai-Mei's data\n",
    "lin_volc_data = pd.ExcelFile(\"/scratch/Dropbox/ipynb/D-O_analysis/volcanic/Lin_MIS3_volcanic_supp.xls\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "5ae75ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_load_data = lin_volc_data.parse('Table S5', header=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7170d98b",
   "metadata": {},
   "outputs": [],
   "source": [
    "gicc05_age = lin_load_data[lin_load_data.columns[1]].copy()\n",
    "volc_load = lin_load_data[lin_load_data.columns[26]].copy()\n",
    "site_loc = lin_load_data[lin_load_data.columns[31]].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "23c0e18d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix 2 southern hemisphere eruptions that we know were in the Southern Hemisphere\n",
    "\n",
    "\n",
    "site_loc[21] = \"SH\"\n",
    "site_loc[62] = \"SH\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8cffda1c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 50072 characters]",
      "text/plain": [
       "<Figure size 2000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# plot chronologies\n",
    "\n",
    "SMALL_SIZE = 18\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 25\n",
    "#\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "\n",
    "#\n",
    "fig = plt.figure(figsize=(20,4))\n",
    "ax = plt.subplot(111)\n",
    "#\n",
    "eleg=True\n",
    "nleg=True\n",
    "sleg=True\n",
    "\n",
    "for e in range(len(gicc05_age)):\n",
    "    if site_loc[e] == \"NHHL\":\n",
    "        if nleg:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='b', label=\"Northern Hemisphere\")\n",
    "            nleg=False\n",
    "        else:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='b')\n",
    "    elif site_loc[e] == \"LL or SH\":\n",
    "        if eleg:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='k', label=\"Equatorial\")\n",
    "            eleg=False\n",
    "        else:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='k')\n",
    "    elif site_loc[e] == \"SH\":\n",
    "        if sleg:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='r', label=\"Southern Hemisphere\")\n",
    "            sleg=False\n",
    "        else:\n",
    "            ax.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='r')\n",
    "#ax.set_xlim(-60000, -24000) \n",
    "#ax.set_xlim(-60000, 0) \n",
    "ax.set_ylabel(\"Load (Tg SO$_4$)\")\n",
    "#ax.set_ylim(0,180)\n",
    "ax.set_xlabel('Time (Year)')\n",
    "#axes[axId].set_ylabel('Load (Tg)')\n",
    "#ax.set_yticks([0,25,50,75,100,125,150,175])\n",
    "ax.set_ylim(0,700)\n",
    "ax.legend()\n",
    "\n",
    "icepaperdir=\"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"random_mis3volc_plus_Lin_large.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n",
    "    plt.savefig(icepaperdir+\"random_mis3volc_plus_Lin_large.pdf\", dpi=300, facecolor=fig.get_facecolor(), transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5fbec726-87c4-46e2-82fd-d8d7a35d07f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "## I want to resample the data so that we take maxima and bin them into even time sections\n",
    "mis3age = gicc05_age[20:]\n",
    "mis3load = volc_load[20:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6dd94d9d-31fc-46bb-984f-4a3835356ec5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36,\n",
       "       37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,\n",
       "       54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,\n",
       "       71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mis3load.index.values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9b5415b3-9ac0-48b5-8aff-48b09d627642",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                         BinStart     BinEnd      BinMid  DeltaTime  \\\n",
      "Bin                                                                   \n",
      "(24633.918, 26189.701]  24633.918  26189.701  25411.8095   1555.783   \n",
      "(26189.701, 27710.506]  26189.701  27710.506  26950.1035   1520.805   \n",
      "(27710.506, 29231.31]   27710.506  29231.310  28470.9080   1520.804   \n",
      "(29231.31, 30752.115]   29231.310  30752.115  29991.7125   1520.805   \n",
      "(30752.115, 32272.919]  30752.115  32272.919  31512.5170   1520.804   \n",
      "(32272.919, 33793.724]  32272.919  33793.724  33033.3215   1520.805   \n",
      "(33793.724, 35314.528]  33793.724  35314.528  34554.1260   1520.804   \n",
      "(35314.528, 36835.333]  35314.528  36835.333  36074.9305   1520.805   \n",
      "(36835.333, 38356.137]  36835.333  38356.137  37595.7350   1520.804   \n",
      "(38356.137, 39876.942]  38356.137  39876.942  39116.5395   1520.805   \n",
      "(39876.942, 41397.746]  39876.942  41397.746  40637.3440   1520.804   \n",
      "(41397.746, 42918.551]  41397.746  42918.551  42158.1485   1520.805   \n",
      "(42918.551, 44439.355]  42918.551  44439.355  43678.9530   1520.804   \n",
      "(44439.355, 45960.16]   44439.355  45960.160  45199.7575   1520.805   \n",
      "(45960.16, 47480.964]   45960.160  47480.964  46720.5620   1520.804   \n",
      "(47480.964, 49001.769]  47480.964  49001.769  48241.3665   1520.805   \n",
      "(49001.769, 50522.573]  49001.769  50522.573  49762.1710   1520.804   \n",
      "(50522.573, 52043.378]  50522.573  52043.378  51282.9755   1520.805   \n",
      "(52043.378, 53564.182]  52043.378  53564.182  52803.7800   1520.804   \n",
      "(53564.182, 55084.987]  53564.182  55084.987  54324.5845   1520.805   \n",
      "(55084.987, 56605.791]  55084.987  56605.791  55845.3890   1520.804   \n",
      "(56605.791, 58126.596]  56605.791  58126.596  57366.1935   1520.805   \n",
      "(58126.596, 59647.4]    58126.596  59647.400  58886.9980   1520.804   \n",
      "\n",
      "                           MaxLoad  \n",
      "Bin                                 \n",
      "(24633.918, 26189.701]  382.280000  \n",
      "(26189.701, 27710.506]         NaN  \n",
      "(27710.506, 29231.31]   185.866667  \n",
      "(29231.31, 30752.115]   327.547667  \n",
      "(30752.115, 32272.919]  265.760069  \n",
      "(32272.919, 33793.724]   93.300000  \n",
      "(33793.724, 35314.528]  224.103333  \n",
      "(35314.528, 36835.333]  180.120670  \n",
      "(36835.333, 38356.137]  413.900231  \n",
      "(38356.137, 39876.942]  151.560000  \n",
      "(39876.942, 41397.746]  315.485856  \n",
      "(41397.746, 42918.551]  306.555667  \n",
      "(42918.551, 44439.355]   91.200000  \n",
      "(44439.355, 45960.16]   544.681000  \n",
      "(45960.16, 47480.964]   316.200000  \n",
      "(47480.964, 49001.769]   76.251500  \n",
      "(49001.769, 50522.573]  307.219250  \n",
      "(50522.573, 52043.378]  104.151083  \n",
      "(52043.378, 53564.182]  352.436504  \n",
      "(53564.182, 55084.987]  100.624000  \n",
      "(55084.987, 56605.791]  653.399500  \n",
      "(56605.791, 58126.596]  127.217500  \n",
      "(58126.596, 59647.4]    168.066333  \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_287523/3024582526.py:19: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.\n",
      "  rebinned_df = df.groupby('Bin').agg(MaxLoad=('Load', 'max'))\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Your actual data (this is just a placeholder for illustration, replace it with your actual dataset)\n",
    "time_series = mis3age\n",
    "\n",
    "load_series = mis3load\n",
    "\n",
    "# Convert the time_series and load_series into a DataFrame\n",
    "df = pd.DataFrame({'Time': time_series, 'Load': load_series})\n",
    "\n",
    "# Define custom bin edges based on time_series values\n",
    "# You can use a fixed number of bins or custom breakpoints.\n",
    "# Here, we use 5 bins for illustration\n",
    "num_bins = 23 # adjust higher until there are no more nan's\n",
    "bins = pd.cut(df['Time'], bins=num_bins)\n",
    "\n",
    "# Add the bin information to the DataFrame\n",
    "df['Bin'] = bins\n",
    "\n",
    "# Group by the 'Bin' and calculate the maximum 'Load' in each bin\n",
    "rebinned_df = df.groupby('Bin').agg(MaxLoad=('Load', 'max'))\n",
    "\n",
    "# Extract the bin edges (lower and upper bounds)\n",
    "rebinned_df['BinStart'] = rebinned_df.index.categories.left\n",
    "rebinned_df['BinEnd'] = rebinned_df.index.categories.right\n",
    "\n",
    "# Calculate the delta time for each bin\n",
    "rebinned_df['DeltaTime'] = rebinned_df['BinEnd'] - rebinned_df['BinStart']\n",
    "\n",
    "# Calculate the middle of the bin (BinMid)\n",
    "rebinned_df['BinMid'] = (rebinned_df['BinStart'] + rebinned_df['BinEnd']) / 2\n",
    "\n",
    "# Reorder the columns\n",
    "rebinned_df = rebinned_df[['BinStart', 'BinEnd', 'BinMid', 'DeltaTime', 'MaxLoad']]\n",
    "\n",
    "# Output the rebinned dataset\n",
    "print(rebinned_df)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "90508cfc-0dd8-47e0-ab09-b3de1f8992b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# only take the data from the nan onwards\n",
    "lin_new_df = rebinned_df[2:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "47c810b9-f24c-4e00-8c23-149b55d791d3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                         BinStart     BinEnd      BinMid  DeltaTime  \\\n",
      "Bin                                                                   \n",
      "(27710.506, 29231.31]   27710.506  29231.310  28470.9080   1520.804   \n",
      "(29231.31, 30752.115]   29231.310  30752.115  29991.7125   1520.805   \n",
      "(30752.115, 32272.919]  30752.115  32272.919  31512.5170   1520.804   \n",
      "(32272.919, 33793.724]  32272.919  33793.724  33033.3215   1520.805   \n",
      "(33793.724, 35314.528]  33793.724  35314.528  34554.1260   1520.804   \n",
      "(35314.528, 36835.333]  35314.528  36835.333  36074.9305   1520.805   \n",
      "(36835.333, 38356.137]  36835.333  38356.137  37595.7350   1520.804   \n",
      "(38356.137, 39876.942]  38356.137  39876.942  39116.5395   1520.805   \n",
      "(39876.942, 41397.746]  39876.942  41397.746  40637.3440   1520.804   \n",
      "(41397.746, 42918.551]  41397.746  42918.551  42158.1485   1520.805   \n",
      "(42918.551, 44439.355]  42918.551  44439.355  43678.9530   1520.804   \n",
      "(44439.355, 45960.16]   44439.355  45960.160  45199.7575   1520.805   \n",
      "(45960.16, 47480.964]   45960.160  47480.964  46720.5620   1520.804   \n",
      "(47480.964, 49001.769]  47480.964  49001.769  48241.3665   1520.805   \n",
      "(49001.769, 50522.573]  49001.769  50522.573  49762.1710   1520.804   \n",
      "(50522.573, 52043.378]  50522.573  52043.378  51282.9755   1520.805   \n",
      "(52043.378, 53564.182]  52043.378  53564.182  52803.7800   1520.804   \n",
      "(53564.182, 55084.987]  53564.182  55084.987  54324.5845   1520.805   \n",
      "(55084.987, 56605.791]  55084.987  56605.791  55845.3890   1520.804   \n",
      "(56605.791, 58126.596]  56605.791  58126.596  57366.1935   1520.805   \n",
      "(58126.596, 59647.4]    58126.596  59647.400  58886.9980   1520.804   \n",
      "\n",
      "                           MaxLoad  \n",
      "Bin                                 \n",
      "(27710.506, 29231.31]   185.866667  \n",
      "(29231.31, 30752.115]   327.547667  \n",
      "(30752.115, 32272.919]  265.760069  \n",
      "(32272.919, 33793.724]   93.300000  \n",
      "(33793.724, 35314.528]  224.103333  \n",
      "(35314.528, 36835.333]  180.120670  \n",
      "(36835.333, 38356.137]  413.900231  \n",
      "(38356.137, 39876.942]  151.560000  \n",
      "(39876.942, 41397.746]  315.485856  \n",
      "(41397.746, 42918.551]  306.555667  \n",
      "(42918.551, 44439.355]   91.200000  \n",
      "(44439.355, 45960.16]   544.681000  \n",
      "(45960.16, 47480.964]   316.200000  \n",
      "(47480.964, 49001.769]   76.251500  \n",
      "(49001.769, 50522.573]  307.219250  \n",
      "(50522.573, 52043.378]  104.151083  \n",
      "(52043.378, 53564.182]  352.436504  \n",
      "(53564.182, 55084.987]  100.624000  \n",
      "(55084.987, 56605.791]  653.399500  \n",
      "(56605.791, 58126.596]  127.217500  \n",
      "(58126.596, 59647.4]    168.066333  \n"
     ]
    }
   ],
   "source": [
    "print(lin_new_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ef251994-68bc-4f5b-bda7-1d5787c65060",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([28470.908 , 29991.7125, 31512.517 , 33033.3215, 34554.126 ,\n",
       "       36074.9305, 37595.735 , 39116.5395, 40637.344 , 42158.1485,\n",
       "       43678.953 , 45199.7575, 46720.562 , 48241.3665, 49762.171 ,\n",
       "       51282.9755, 52803.78  , 54324.5845, 55845.389 , 57366.1935,\n",
       "       58886.998 ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_new_df[\"BinMid\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "fd618abe-5db1-4283-afb1-0bd0c1153278",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[<matplotlib.lines.Line2D at 0x7ee80d8d5010>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d5250>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d55b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d58e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d5c10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d5ee0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80dcfee10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d943710>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80da7dd90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d99aea0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d9a9bb0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80dac4560>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d82a7e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d88f7d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d88fa10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d88e7e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d40b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d6270>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d6570>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d68a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d8d6b70>]]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15456 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[ plt.plot([-lin_new_df[\"BinMid\"].values[e],-lin_new_df[\"BinMid\"].values[e]], [0, lin_new_df[\"MaxLoad\"].values[e]],color='k') \n",
    " for e in range(len(lin_new_df[\"BinMid\"])) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "c9f90e27-37ed-495b-a35a-c7694d09a1cc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.2470981072803381\n",
      "171.41284059910026\n",
      "97.03317998392464\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 50168 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(lin_new_df[\"MaxLoad\"], bins=5, range=(100, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "#ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Normalized Occurrences\", 16)\n",
    "\n",
    "#ax.set_title('Levy-stable distribution fit for Sigl and Lin data', size=16)\n",
    "ax.set_ylim(0, 0.01)\n",
    "plt.legend(fontsize=12)\n",
    "\n",
    "def main(rvs):\n",
    "    shape, loc, scale = gev.fit(rvs)\n",
    "    return shape, loc, scale\n",
    "\n",
    "\n",
    "shape, loc, scale = main(lin_new_df[\"MaxLoad\"].values)\n",
    "\n",
    "print(shape)\n",
    "print(loc)\n",
    "print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "hist, bins = np.histogram(lin_new_df[\"MaxLoad\"].values, bins=5, range=(0, 600), density=True)\n",
    "plt.bar(bins[:-1], hist, width = 2, align='edge')\n",
    "\n",
    "plt.plot(xx, yy, 'ro')\n",
    "\n",
    "# these are the GEV parameters for the massaged Lin data\n",
    "lin_ml_loc = loc \n",
    "lin_ml_scale = scale\n",
    "lin_ml_shape = shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "c69284d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigl_hol_volc_data = xr.open_dataset(\"/home/guido/Dropbox/ipynb/D-O_analysis/volcanic/HolVol_volcanic_stratospheric_sulfur_injection_v1.0.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "54665d95",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 46236 characters]",
      "text/plain": [
       "<Figure size 2000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot chronologies\n",
    "\n",
    "SMALL_SIZE = 18\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 25\n",
    "#\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "\n",
    "#\n",
    "fig = plt.figure(figsize=(20,4))\n",
    "ax = plt.subplot(111)\n",
    "#\n",
    "for e in range(len(sigl_hol_volc_data.vssi)):\n",
    "    if sigl_hol_volc_data.lat[e] > 20.0:\n",
    "        ax.plot([sigl_hol_volc_data.yearCE[e],sigl_hol_volc_data.yearCE[e]],[0,sigl_hol_volc_data.vssi[e]*3.0],color='b')\n",
    "    elif sigl_hol_volc_data.lat[e] >= -20 and sigl_hol_volc_data.lat[e] <= 20:\n",
    "        ax.plot([sigl_hol_volc_data.yearCE[e],sigl_hol_volc_data.yearCE[e]],[0,sigl_hol_volc_data.vssi[e]*3.0],color='k')\n",
    "    elif sigl_hol_volc_data.lat[e] < -20.:\n",
    "        ax.plot([sigl_hol_volc_data.yearCE[e],sigl_hol_volc_data.yearCE[e]],[0,sigl_hol_volc_data.vssi[e]*3.0],color='r')\n",
    "#ax.set_xlim(24000,60000) \n",
    "ax.set_ylabel('Load (Tg SO$_4$)')\n",
    "#ax.set_ylim(0,180)\n",
    "ax.set_xlabel('Year B2K')\n",
    "#axes[axId].set_ylabel('Load (Tg)')\n",
    "#ax.set_yticks([0,25,50,75,100,125,150,175])\n",
    "\n",
    "\n",
    "#for e in range(len(year)):\n",
    "#    if flag[e] == 1:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='k')\n",
    "#    elif flag[e] == 2:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='b')\n",
    "#    elif flag[e] == 3:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='r')\n",
    "icepaperdir=\"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"random_mis3volc_plus_Sigl2022.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n",
    "    plt.savefig(icepaperdir+\"random_mis3volc_plus_Sigl2022.pdf\", dpi=300, facecolor=fig.get_facecolor(), transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "cedcdb99",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigl2022 = sigl_hol_volc_data.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "c8c78ac0-b00a-407b-8aad-ee43e818c0b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "mass = sigl2022.vssi*3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "cb4a5919-b064-4bea-8d2f-ce7338301992",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 13.32,   3.12,  27.27, ...,  27.39, 173.7 ,   8.67], shape=(1384,))"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mass.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4eae18f5-57c7-4678-8ec7-ec4079c96ef6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([1.032e+03, 1.650e+02, 8.100e+01, 3.400e+01, 2.700e+01, 9.000e+00,\n",
       "        8.000e+00, 5.000e+00, 2.000e+00, 3.000e+00, 3.000e+00, 3.000e+00,\n",
       "        2.000e+00, 0.000e+00, 1.000e+00, 1.000e+00, 3.000e+00, 0.000e+00,\n",
       "        0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n",
       "        0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 1.000e+00, 0.000e+00,\n",
       "        0.000e+00, 2.000e+00, 1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00,\n",
       "        0.000e+00, 1.000e+00, 0.000e+00, 0.000e+00]),\n",
       " array([  0.,  15.,  30.,  45.,  60.,  75.,  90., 105., 120., 135., 150.,\n",
       "        165., 180., 195., 210., 225., 240., 255., 270., 285., 300., 315.,\n",
       "        330., 345., 360., 375., 390., 405., 420., 435., 450., 465., 480.,\n",
       "        495., 510., 525., 540., 555., 570., 585., 600.]),\n",
       " <BarContainer object of 40 artists>)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 18772 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(sigl2022.vssi*3, bins=np.linspace(0,600,41), edgecolor=\"red\", color=\"grey\")\n",
    "#plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "bb7ed41a-4838-406a-8aed-7cc227b9bc87",
   "metadata": {},
   "outputs": [],
   "source": [
    "a=gev.fit(mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "96e51a11-0906-4c21-b382-ad2dc242a32a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.21988720690114583\n",
      "12.749730029827154\n",
      "3.448963234019624\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ee80475a5a0>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 29164 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# gev example\n",
    "def main(rvs):\n",
    "    shape, loc, scale = gev.fit(rvs)\n",
    "    return shape, loc, scale\n",
    "\n",
    "rvs = [9.4, 38.0, 12.5, 35.3, 17.6, 12.9, 12.4, 19.6, 15.0, 13.2, 12.3, 16.9, 16.9, 29.4, 13.6, 11.1, 8.0, 16.6, 12.0, 13.1, 9.1, 9.7, 21.0, 11.2, 14.4, 18.8, 14.0, 19.9, 12.4, 10.8, 21.6, 15.4, 17.4, 14.8, 22.7, 11.5, 10.5, 11.8, 12.4, 16.6, 11.7, 12.9, 17.8]\n",
    "\n",
    "shape, loc, scale = main(rvs)\n",
    "\n",
    "print(shape)\n",
    "print(loc)\n",
    "print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "\n",
    "xx = np.linspace(l+0.00001, l+0.00001+35, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "hist, bins = np.histogram(rvs, bins=12, range=(-0.5, 23.5), density=True)\n",
    "plt.bar(bins[:-1], hist, width = 2, align='edge')\n",
    "\n",
    "plt.plot(xx, yy, 'ro')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "126a3f53-8c21-45fb-927d-4cd6faf67ce5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(rvs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "74f2b059-0520-40c0-a7fd-cb983575f277",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.344260778639558\n",
      "2.498705531121373\n",
      "3.163418016409622\n",
      "0.1454284983635481\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 76976 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(14,9))\n",
    "ax = plt.subplot(111)\n",
    "shape, loc, scale = gev.fit(mass)\n",
    "\n",
    "print(shape)\n",
    "print(loc)\n",
    "print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "print (l)\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "hist, bins = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins[:-1], hist, width = 10, align='edge', label=\"Occurance #\")\n",
    "\n",
    "ax.plot(xx, yy, 'orange', linewidth = 3.5, label = \"GEV fit\")\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurances\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\",size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurances\",  16)\n",
    "fshape = np.around(shape,2)\n",
    "floc = np.around(loc,2)\n",
    "fscale = np.around(scale,2)\n",
    "ax.set_title(r'$GEV(x) = exp[-\\{1+s(x − a)/b\\}^{-1/s}]$  s='+str(fshape)+'  a='+str(floc)+'  b='+str(fscale),size=16)\n",
    "\n",
    "plt.legend(fontsize=16)\n",
    "icepaperdir=\"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure_GEV_sigl2022.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "    plt.savefig(icepaperdir+\"Figure_GEV_sigl2022.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "295fea28",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean: 0.0000, Variance: inf, Skewness: nan, Kurtosis: nan\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 96952 characters]",
      "text/plain": [
       "<Figure size 1200x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CDF accuracy check: True\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(12, 8))\n",
    "\n",
    "# Set parameters for the Levy-stable distribution\n",
    "alpha, beta = 1.8, -0.5\n",
    "\n",
    "# Calculate the first four moments\n",
    "mean, var, skew, kurt = levy_stable.stats(alpha, beta, moments='mvsk')\n",
    "print(f\"Mean: {mean:.4f}, Variance: {var:.4f}, Skewness: {skew:.4f}, Kurtosis: {kurt:.4f}\")\n",
    "\n",
    "# Generate x values for plotting\n",
    "x = np.linspace(levy_stable.ppf(0.01, alpha, beta),\n",
    "                levy_stable.ppf(0.99, alpha, beta), 100)\n",
    "\n",
    "# Plot the probability density function (PDF)\n",
    "ax.plot(x, levy_stable.pdf(x, alpha, beta),\n",
    "        'r-', lw=3, alpha=0.6, label='Levy-stable PDF')\n",
    "\n",
    "# Create a frozen distribution and plot its PDF\n",
    "rv = levy_stable(alpha, beta)\n",
    "ax.plot(x, rv.pdf(x), 'k--', lw=2, label='Frozen PDF')\n",
    "\n",
    "# Generate random numbers\n",
    "r = levy_stable.rvs(alpha, beta, size=10000)\n",
    "\n",
    "# Plot the histogram of random numbers\n",
    "ax.hist(r, density=True, bins='auto', histtype='stepfilled', alpha=0.2, label='Histogram')\n",
    "\n",
    "# Set plot limits and labels\n",
    "ax.set_xlim([x[0], x[-1]])\n",
    "ax.set_xlabel('X', fontsize=14)\n",
    "ax.set_ylabel('Probability Density', fontsize=14)\n",
    "ax.set_title(f'Levy-stable Distribution (α={alpha}, β={beta})', fontsize=16)\n",
    "\n",
    "# Add legend\n",
    "ax.legend(loc='best', frameon=False, fontsize=12)\n",
    "\n",
    "# Show the plot\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Check accuracy of CDF and PPF\n",
    "vals = levy_stable.ppf([0.001, 0.5, 0.999], alpha, beta)\n",
    "print(\"CDF accuracy check:\", np.allclose([0.001, 0.5, 0.999], levy_stable.cdf(vals, alpha, beta)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "32534b18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1454284983635481\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ee802b6b260>"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 64420 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(14,9))\n",
    "ax = plt.subplot(111)\n",
    "#shape, loc, scale = gev.fit(mass)\n",
    "\n",
    "#print(shape)\n",
    "#print(loc)\n",
    "#print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "print (l)\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "hist, bins = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins[:-1], hist, width = 10, align='edge', label=\"Occurance #\")\n",
    "\n",
    "ax.plot(xx, yy, 'orange', linewidth = 3.5, label = \"GEV fit\")\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurances\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\",size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurances\",  16)\n",
    "fshape = np.around(shape,2)\n",
    "floc = np.around(loc,2)\n",
    "fscale = np.around(scale,2)\n",
    "#ax.set_title(r'$GEV(x) = exp[-\\{1+s(x − a)/b\\}^{-1/s}]$  s='+str(fshape)+'  a='+str(floc)+'  b='+str(fscale),size=16)\n",
    "#ax.set_title(r'$alpha(x) = exp[-\\{1+s(x − a)/b\\}^{-1/s}]$  s='+str(fshape)+'  a='+str(floc)+'  b='+str(fscale),size=16)\n",
    "\n",
    "plt.legend(fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "e84880c7-5f3a-4ca9-8de4-c04aad941505",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1454284983635481\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 43996 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "# Fit levy_stable distribution to Sigl data\n",
    "#params = levy_stable.fit(mass)\n",
    "#alpha, beta, loc, scale = params\n",
    "\n",
    "#print(f\"Sigl data fit - Alpha: {alpha:.2f}, Beta: {beta:.2f}, Location: {loc:.2f}, Scale: {scale:.2f}\")\n",
    "\n",
    "# Generate x values for plotting\n",
    "x = np.linspace(0, 600, num=1000)\n",
    "\n",
    "# Plot the probability density function (PDF) for Sigl data\n",
    "#pdf_sigl = levy_stable.pdf(x, *params)\n",
    "#ax.plot(x, pdf_sigl, 'r-', linewidth=3, alpha=0.6, label='Levy-stable fit (Sigl)')\n",
    "\n",
    "# Plot histogram for Sigl data\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_sigl[:-1], hist_sigl, width=10, align='edge', alpha=0.3, color='red', label=\"Sigl data\")\n",
    "\n",
    "# Fit levy_stable distribution to Lin data\n",
    "#params_lin = levy_stable.fit(volc_load)\n",
    "#alpha_lin, beta_lin, loc_lin, scale_lin = params_lin\n",
    "\n",
    "#print(f\"Lin data fit - Alpha: {alpha_lin:.2f}, Beta: {beta_lin:.2f}, Location: {loc_lin:.2f}, Scale: {scale_lin:.2f}\")\n",
    "\n",
    "# Plot the probability density function (PDF) for Lin data\n",
    "#pdf_lin = levy_stable.pdf(x, *params_lin)\n",
    "#ax.plot(x, pdf_lin, 'b-', linewidth=3, alpha=0.6, label='Levy-stable fit (Lin)')\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(volc_load, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "#ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Normalized Occurrences\", 16)\n",
    "\n",
    "#ax.set_title('Levy-stable distribution fit for Sigl and Lin data', size=16)\n",
    "ax.set_ylim(0, 0.01)\n",
    "plt.legend(fontsize=12)\n",
    "\n",
    "#sigl gev\n",
    "\n",
    "l = loc + scale / shape\n",
    "print (l)\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure_LevyStable_sigl_lin_2022.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "    plt.savefig(icepaperdir+\"Figure_LevyStable_sigl_lin_2022.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "090b7592-f7b0-4115-9292-ebfeaef38ce5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20    24668.896552\n",
       "21    25460.344828\n",
       "22    25759.333333\n",
       "23    25939.870968\n",
       "24    27796.666667\n",
       "          ...     \n",
       "80    58355.241379\n",
       "81    59180.090909\n",
       "82    59316.551724\n",
       "83    59545.142857\n",
       "84    59647.400000\n",
       "Name: Unnamed: 1, Length: 65, dtype: float64"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gicc05_age[20:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "d0e05f91-3dd0-484a-968a-2afd2f4b5b7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mis3_lin_volc_load = volc_load[20:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "0add271c-d670-46b3-b45e-41000a72c9e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted parameters:\n",
      "Shape: -1.34426, Location: 2.49871, Scale: 3.16342\n"
     ]
    }
   ],
   "source": [
    "sigl_shape, sigl_loc, sigl_scale = gev.fit(mass)\n",
    "print(f\"Fitted parameters:\\nShape: {sigl_shape:.5f}, Location: {sigl_loc:.5f}, Scale: {sigl_scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "745c8748-6f09-4b6e-a4e0-d738d540a714",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted parameters:\n",
      "Shape: -0.46851, Location: 107.23891, Scale: 51.76101\n"
     ]
    }
   ],
   "source": [
    "lin_shape, lin_loc, lin_scale = gev.fit(mis3_lin_volc_load)\n",
    "print(f\"Fitted parameters:\\nShape: {lin_shape:.5f}, Location: {lin_loc:.5f}, Scale: {lin_scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "a10bbeb9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 77800 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "\n",
    "# Plot histogram for Sigl data\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_sigl[:-1], hist_sigl, width=10, align='edge', alpha=0.3, color='red', label=\"Sigl data\")\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(mis3_lin_volc_load , bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "#ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Normalized Occurrences\", 16)\n",
    "\n",
    "#ax.set_title('Levy-stable distribution fit for Sigl and Lin data', size=16)\n",
    "ax.set_ylim(0, 0.01)\n",
    "plt.legend(fontsize=12)\n",
    "\n",
    "#sigl gev\n",
    "\n",
    "sigl_l = sigl_loc + sigl_scale / sigl_shape\n",
    "#print (sigl_l)\n",
    "sigl_xx = np.linspace(sigl_l+0.00001, sigl_l+0.00001+600, num=71)\n",
    "sigl_yy = gev.pdf(xx, sigl_shape, sigl_loc, sigl_scale)\n",
    "ax.plot(sigl_xx, sigl_yy, 'red', linewidth = 1.5, label = \"Sigl GEV fit\")\n",
    "\n",
    "#lin gev\n",
    "\n",
    "lin_l = lin_loc + lin_scale / lin_shape\n",
    "#print (lin_l)\n",
    "lin_xx = np.linspace(lin_l+0.00001, lin_l+0.00001+600, num=71)\n",
    "lin_yy = gev.pdf(xx, lin_shape, lin_loc, lin_scale)\n",
    "ax.plot(lin_xx, lin_yy, 'blue', linewidth = 1.5, label = \"Sigl GEV fit\")\n",
    "\n",
    "\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure_LevyStable_sigl_lin_2022.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "    plt.savefig(icepaperdir+\"Figure_LevyStable_sigl_lin_2022.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9bd492a0-90e1-40bc-8cfa-662913b67ce6",
   "metadata": {},
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "\n",
    "# Plot histogram for Sigl data\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_sigl[:-1], hist_sigl, width=10, align='edge', alpha=0.3, color='red', label=\"Sigl data\")\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(volc_load, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurrences\", 16)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "7ae33828",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1e-05, 0.1)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 38328 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "\n",
    "# Plot histogram for Sigl data\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_sigl[:-1], hist_sigl, width=10, align='edge', alpha=0.3, color='red', label=\"Sigl data\")\n",
    "\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurrences\", 16)\n",
    "ax.set_ylim(1.e-5,1.e-1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "4340ee6b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1e-05, 0.1)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 38372 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import levy_stable\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(volc_load, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurrences\", 16)\n",
    "ax.set_ylim(1.e-5,1.e-1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "b04c84e4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 56792 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "# Calculate histograms for both datasets\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "hist_lin, bins_lin = np.histogram(volc_load, bins=31, range=(0, 600), density=True)\n",
    "\n",
    "# Plot stacked histograms\n",
    "width = 10  # width of bars\n",
    "ax.bar(bins_sigl[:-1], hist_sigl, width=width, align='edge', alpha=0.5, color='red', label=\"Sigl data\")\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=width, align='edge', alpha=0.5, color='blue', label=\"Lin data\", bottom=hist_sigl)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurrences\", 16)\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylim(1.e-5,1.e-1)\n",
    "plt.title(\"Stacked Histograms of Sigl and Lin Data\", size=16)\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "0c4d4a71",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([153.8648    , 199.78416667, 152.96892436,  82.682     ,\n",
       "       103.6745    , 200.78011699,  28.2541    ,  77.05      ,\n",
       "       220.51170852,  77.8       ,  95.502     , 103.75      ,\n",
       "       186.15797321, 132.11666667,  89.66806591, 156.23333333,\n",
       "        84.188     , 105.72033333, 152.70372964, 192.7       ,\n",
       "       160.65783333, 382.28      , 233.06298818,  79.73183333,\n",
       "       168.02783333, 117.58233333, 185.86666667, 327.54766667,\n",
       "       154.2655    , 155.23359092, 265.76006882,  86.2       ,\n",
       "        93.3       ,  59.2       ,  71.73      , 224.10333333,\n",
       "        89.53333333, 180.12067008,  86.60866667, 120.3885    ,\n",
       "       144.74233333, 413.90023062,  90.19      , 145.6728113 ,\n",
       "       151.56      , 224.14688247, 168.09201015, 315.48585591,\n",
       "        99.35      , 306.55566667, 100.1365    , 185.24333333,\n",
       "       161.52759461,  63.004     ,  91.2       , 161.9       ,\n",
       "        65.738     , 340.65      , 104.61666667, 544.681     ,\n",
       "       114.62889021, 165.8831513 , 316.2       , 197.87237332,\n",
       "       107.13016667,  76.2515    , 307.21924954, 114.8       ,\n",
       "        77.2       , 104.15108331, 352.43650396, 140.04262788,\n",
       "        70.6       , 100.624     ,  52.56666667, 653.3995    ,\n",
       "       127.2175    ,  60.136     ,  87.86666667, 168.06633333,\n",
       "        92.58933333,  57.12826667,  91.43333333,  96.656     ,\n",
       "       145.56333333])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "volc_load.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "57d4a3be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0     10480.926070\n",
       "1     11304.500000\n",
       "2     12170.400000\n",
       "3     12240.690909\n",
       "4     12754.814815\n",
       "          ...     \n",
       "80    58355.241379\n",
       "81    59180.090909\n",
       "82    59316.551724\n",
       "83    59545.142857\n",
       "84    59647.400000\n",
       "Name: Unnamed: 1, Length: 85, dtype: float64"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gicc05_age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b605e530",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "c73bcda3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 68860 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "# Combine the datasets\n",
    "combined_data = np.concatenate([mass, volc_load])\n",
    "\n",
    "# Calculate histograms for all datasets\n",
    "hist_sigl, bins = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "hist_lin, _ = np.histogram(volc_load, bins=bins, density=True)\n",
    "hist_combined, _ = np.histogram(combined_data, bins=bins, density=True)\n",
    "\n",
    "# Plot side-by-side histograms\n",
    "width = 5  # width of bars (half of before to fit side by side)\n",
    "ax.bar(bins[:-1], hist_sigl, width=width, align='edge', alpha=0.5, color='red', label=\"Sigl data\")\n",
    "ax.bar(bins[:-1]+width, hist_lin, width=width, align='edge', alpha=0.5, color='blue', label=\"Lin data\")\n",
    "\n",
    "# Plot the outline of the combined histogram\n",
    "ax.step(bins[:-1], hist_combined, where='post', color='black', linewidth=2, label=\"Combined\")\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Log Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurrences\", 16)\n",
    "\n",
    "ax.legend()\n",
    "ax.set_ylim(1.e-5, 1.e-1)\n",
    "plt.title(\"Side-by-Side Histograms of Sigl and Lin Data with Combined Outline\", size=16)\n",
    "plt.tight_layout()\n",
    "\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure_SideBySide_Histogram_sigl_lin_combined_2022.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "    plt.savefig(icepaperdir+\"Figure_SideBySide_Histogram_sigl_lin_combined_2022.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "87c4820a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 77800 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10-year event: 48.61 Tg SO4\n",
      "100-year event: 1141.09 Tg SO4\n",
      "1000-year event: 25360.97 Tg SO4\n"
     ]
    }
   ],
   "source": [
    "# this is incorrect because of timescale interval is not 1 year\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "# Assuming 'mass' is your data array\n",
    "# Fit GEV distribution\n",
    "shape, loc, scale = stats.genextreme.fit(mass)\n",
    "\n",
    "# Define return periods (in years)\n",
    "return_periods = np.logspace(0, 3, 100)  # 1 to 1000 years, 100 points\n",
    "\n",
    "# Calculate return levels\n",
    "return_levels = stats.genextreme.isf(1/return_periods, c=shape, loc=loc, scale=scale)\n",
    "\n",
    "# Plot return time graph\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.semilogx(return_periods, return_levels, 'b-')\n",
    "ax.set_xlabel('Return Period (years)')\n",
    "ax.set_ylabel('Return Level (Tg SO$_4$)')\n",
    "ax.set_title('Return Time Graph for Volcanic Eruptions')\n",
    "ax.grid(True)\n",
    "\n",
    "# Add some reference lines\n",
    "ax.axhline(stats.genextreme.isf(1/10, c=shape, loc=loc, scale=scale), color='r', linestyle='--', label='10-year event')\n",
    "ax.axhline(stats.genextreme.isf(1/100, c=shape, loc=loc, scale=scale), color='g', linestyle='--', label='100-year event')\n",
    "\n",
    "ax.legend()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# Print some key return levels\n",
    "print(f\"10-year event: {stats.genextreme.isf(1/10, c=shape, loc=loc, scale=scale):.2f} Tg SO4\")\n",
    "print(f\"100-year event: {stats.genextreme.isf(1/100, c=shape, loc=loc, scale=scale):.2f} Tg SO4\")\n",
    "print(f\"1000-year event: {stats.genextreme.isf(1/1000, c=shape, loc=loc, scale=scale):.2f} Tg SO4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "fa987467-190b-4860-8bd1-fc59b763782f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_return_level(quantile, loc, scale, shape):\n",
    "    level = loc + scale / shape * (1 - (-np.log(quantile)) ** (shape))\n",
    "    return level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "2b3cec87-517e-439e-813a-29eb6f5de418",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in notebooks */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base rgba(0, 0, 0, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, white)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 15))\n",
       "  );\n",
       "}\n",
       "\n",
       "html[theme=\"dark\"],\n",
       "html[data-theme=\"dark\"],\n",
       "body[data-theme=\"dark\"],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, #111111)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 15))\n",
       "  );\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "  line-height: 1.6;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-obj-name,\n",
       ".xr-group-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-group-name::before {\n",
       "  content: \"📁\";\n",
       "  padding-right: 0.3em;\n",
       "}\n",
       "\n",
       ".xr-group-name,\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 0 20px 0 20px;\n",
       "  margin-block-start: 0;\n",
       "  margin-block-end: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: inline-block;\n",
       "  opacity: 0;\n",
       "  height: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "  border: 2px solid transparent !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:focus + label {\n",
       "  border: 2px solid var(--xr-font-color0) !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: \"►\";\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: \"▼\";\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-top: 4px;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-group-box {\n",
       "  display: inline-grid;\n",
       "  grid-template-columns: 0px 20px auto;\n",
       "  width: 100%;\n",
       "}\n",
       "\n",
       ".xr-group-box-vline {\n",
       "  grid-column-start: 1;\n",
       "  border-right: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "  width: 0px;\n",
       "}\n",
       "\n",
       ".xr-group-box-hline {\n",
       "  grid-column-start: 2;\n",
       "  grid-row-start: 1;\n",
       "  height: 1em;\n",
       "  width: 20px;\n",
       "  border-bottom: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-group-box-contents {\n",
       "  grid-column-start: 3;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: \"(\";\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: \")\";\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: \",\";\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  border-color: var(--xr-background-color-row-odd);\n",
       "  margin-bottom: 0;\n",
       "  padding-top: 2px;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "  border-color: var(--xr-background-color-row-even);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  border-top: 2px dotted var(--xr-background-color);\n",
       "  padding-bottom: 20px !important;\n",
       "  padding-top: 10px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in + label,\n",
       ".xr-var-data-in + label,\n",
       ".xr-index-data-in + label {\n",
       "  padding: 0 1px;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-data > pre,\n",
       ".xr-index-data > pre,\n",
       ".xr-var-data > table > tbody > tr {\n",
       "  background-color: transparent !important;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked + label > .xr-icon-file-text2,\n",
       ".xr-var-data-in:checked + label > .xr-icon-database,\n",
       ".xr-index-data-in:checked + label > .xr-icon-database {\n",
       "  color: var(--xr-font-color0);\n",
       "  filter: drop-shadow(1px 1px 5px var(--xr-font-color2));\n",
       "  stroke-width: 0.8px;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt; Size: 100kB\n",
       "Dimensions:     (nerup: 1384)\n",
       "Dimensions without coordinates: nerup\n",
       "Data variables:\n",
       "    year        (nerup) float64 11kB ...\n",
       "    yearCE      (nerup) float64 11kB ...\n",
       "    month       (nerup) float64 11kB ...\n",
       "    day         (nerup) float64 11kB ...\n",
       "    lat         (nerup) float64 11kB ...\n",
       "    vssi        (nerup) float64 11kB 4.44 1.04 9.09 0.28 ... 0.22 9.13 57.9 2.89\n",
       "    hemi        (nerup) float64 11kB ...\n",
       "    sigma_vssi  (nerup) float64 11kB ...\n",
       "    duration    (nerup) float64 11kB ...\n",
       "Attributes:\n",
       "    title:    Ice core-inferred volcanic stratospheric sulfur injection for t...\n",
       "    history:  Created Tue Jun  9 23:01:18 2020\n",
       "    source:   Sigl et al., in preparation</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-02d1cea4-e56c-4ae4-b272-0a96109f449b' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-02d1cea4-e56c-4ae4-b272-0a96109f449b' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span>nerup</span>: 1384</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-c6639a1a-da10-4c54-b58e-7ff47ac57c4b' class='xr-section-summary-in' type='checkbox'  checked><label for='section-c6639a1a-da10-4c54-b58e-7ff47ac57c4b' class='xr-section-summary' >Data variables: <span>(9)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>year</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-2f983127-3bb8-444c-9702-7de0b5d2a3b8' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-2f983127-3bb8-444c-9702-7de0b5d2a3b8' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-4d917bee-9905-4e59-b9a8-c2449be55cd7' class='xr-var-data-in' type='checkbox'><label for='data-4d917bee-9905-4e59-b9a8-c2449be55cd7' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>year (ISO)</dd><dt><span>long_name :</span></dt><dd>year of eruption (ISO 8601, including year 0)</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>yearCE</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-3a7e38d3-b006-48b0-8c4a-4cd8b7ae7e5e' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-3a7e38d3-b006-48b0-8c4a-4cd8b7ae7e5e' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-fed8b942-ba6e-4643-9904-13f7665ef367' class='xr-var-data-in' type='checkbox'><label for='data-fed8b942-ba6e-4643-9904-13f7665ef367' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>year (BCE/CE)</dd><dt><span>long_name :</span></dt><dd>year of eruption (BCE/CE, no year 0)</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>month</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-804074b1-cfc1-4f68-be68-c73a23452c96' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-804074b1-cfc1-4f68-be68-c73a23452c96' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-73db8d88-2ede-411a-8ff8-49f83b8d1569' class='xr-var-data-in' type='checkbox'><label for='data-73db8d88-2ede-411a-8ff8-49f83b8d1569' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>month</dd><dt><span>long_name :</span></dt><dd>month of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>day</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-e51e57eb-c916-4373-828d-2540ea6166c4' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-e51e57eb-c916-4373-828d-2540ea6166c4' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-cfbc4118-7092-473b-a614-1b2f474b94dc' class='xr-var-data-in' type='checkbox'><label for='data-cfbc4118-7092-473b-a614-1b2f474b94dc' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>day</dd><dt><span>long_name :</span></dt><dd>day of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>lat</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-fb81c738-d1e1-414e-a43d-66ece7e68e44' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-fb81c738-d1e1-414e-a43d-66ece7e68e44' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-c7f1994f-a34e-436a-9b31-3b0232119c2a' class='xr-var-data-in' type='checkbox'><label for='data-c7f1994f-a34e-436a-9b31-3b0232119c2a' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>long_name :</span></dt><dd>latitude of eruption</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>vssi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>4.44 1.04 9.09 ... 9.13 57.9 2.89</div><input id='attrs-d2763959-b447-4bc3-9749-1464a46f9586' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-d2763959-b447-4bc3-9749-1464a46f9586' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-7c6d5c2b-9399-47c7-a558-4c8595095b67' class='xr-var-data-in' type='checkbox'><label for='data-7c6d5c2b-9399-47c7-a558-4c8595095b67' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>TgS</dd><dt><span>long_name :</span></dt><dd>volcanic stratospheric sulfur injection of eruption</dd></dl></div><div class='xr-var-data'><pre>array([ 4.44,  1.04,  9.09, ...,  9.13, 57.9 ,  2.89], shape=(1384,))</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>hemi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-4de8d4b3-5586-4869-8a1b-f6f8d90d9b94' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-4de8d4b3-5586-4869-8a1b-f6f8d90d9b94' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f8f6503f-c6d3-4eb9-9c32-5970ee5e1410' class='xr-var-data-in' type='checkbox'><label for='data-f8f6503f-c6d3-4eb9-9c32-5970ee5e1410' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd></dd><dt><span>long_name :</span></dt><dd>hemispheric asymmetry (NH/SH) of aerosol spread for tropical eruptions based on ratio of Greenland to Antarctic flux</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>sigma_vssi</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-4c815db3-2ee2-41e3-897d-d2623ffba1ac' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-4c815db3-2ee2-41e3-897d-d2623ffba1ac' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-89ba2140-5112-4402-819a-fc7815e9a072' class='xr-var-data-in' type='checkbox'><label for='data-89ba2140-5112-4402-819a-fc7815e9a072' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>TgS</dd><dt><span>long_name :</span></dt><dd>Random uncertainty (1 sigma) in VSSI</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>duration</span></div><div class='xr-var-dims'>(nerup)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-56117da4-e889-4f04-a437-05704da87f4f' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-56117da4-e889-4f04-a437-05704da87f4f' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-7efbca4e-d279-4977-bfdb-f2e7987dddec' class='xr-var-data-in' type='checkbox'><label for='data-7efbca4e-d279-4977-bfdb-f2e7987dddec' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>months</dd><dt><span>long_name :</span></dt><dd>duration of effusive eruptions in months</dd></dl></div><div class='xr-var-data'><pre>[1384 values with dtype=float64]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-6fbadc50-2d4b-4880-8fad-5c3eab7fa46c' class='xr-section-summary-in' type='checkbox'  checked><label for='section-6fbadc50-2d4b-4880-8fad-5c3eab7fa46c' class='xr-section-summary' >Attributes: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>title :</span></dt><dd>Ice core-inferred volcanic stratospheric sulfur injection for the Holocene</dd><dt><span>history :</span></dt><dd>Created Tue Jun  9 23:01:18 2020</dd><dt><span>source :</span></dt><dd>Sigl et al., in preparation</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset> Size: 100kB\n",
       "Dimensions:     (nerup: 1384)\n",
       "Dimensions without coordinates: nerup\n",
       "Data variables:\n",
       "    year        (nerup) float64 11kB ...\n",
       "    yearCE      (nerup) float64 11kB ...\n",
       "    month       (nerup) float64 11kB ...\n",
       "    day         (nerup) float64 11kB ...\n",
       "    lat         (nerup) float64 11kB ...\n",
       "    vssi        (nerup) float64 11kB 4.44 1.04 9.09 0.28 ... 0.22 9.13 57.9 2.89\n",
       "    hemi        (nerup) float64 11kB ...\n",
       "    sigma_vssi  (nerup) float64 11kB ...\n",
       "    duration    (nerup) float64 11kB ...\n",
       "Attributes:\n",
       "    title:    Ice core-inferred volcanic stratospheric sulfur injection for t...\n",
       "    history:  Created Tue Jun  9 23:01:18 2020\n",
       "    source:   Sigl et al., in preparation"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "b0399c81-1b92-46df-bc09-3f4bf7bbcfe2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1892.,  1886.,  1884., ..., -9491., -9493., -9496.], shape=(1384,))"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022.year.values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "d2bdcd03-a3a3-4f7f-8866-561169ddca93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1892.,  1886.,  1884., ..., -9492., -9494., -9497.], shape=(1384,))"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022.yearCE.values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "37ad518e-48a4-486d-bcb7-a548c426390b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(1.0)"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(sigl2022.day.values)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "6bf3e7c2-6d9f-47cc-8280-fed72bb32626",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1892., 1886., 1884., 1876., 1863., 1861., 1855., 1844., 1835.,\n",
       "       1831., 1822., 1815., 1809., 1795., 1783., 1767., 1762., 1754.,\n",
       "       1739., 1729., 1719., 1702., 1694., 1673., 1666., 1661., 1646.,\n",
       "       1641., 1634., 1621., 1600., 1594., 1591., 1586., 1577., 1568.,\n",
       "       1559., 1555., 1540., 1527., 1512., 1505., 1495., 1490., 1480.,\n",
       "       1477., 1477., 1470., 1458., 1453., 1447., 1428., 1418., 1414.,\n",
       "       1390., 1385., 1381., 1378., 1344., 1342., 1336., 1330., 1290.,\n",
       "       1284., 1276., 1269., 1265., 1258., 1241., 1234., 1229., 1227.,\n",
       "       1214., 1209., 1205., 1199., 1191., 1188., 1181., 1181., 1171.,\n",
       "       1152., 1146., 1126., 1115., 1112., 1109., 1091., 1081., 1065.,\n",
       "       1044., 1038., 1029., 1022., 1015., 1008., 1006., 1004., 1002.,\n",
       "       1001.])"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022.year.values[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "cb81569e-9580-4bb5-b835-c24d5150027f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-8887., -8896., -8905., -8933., -8944., -8963., -8964., -8966.,\n",
       "       -8969., -8970., -8982., -8990., -8994., -9001., -9009., -9013.,\n",
       "       -9017., -9021., -9030., -9035., -9041., -9045., -9049., -9066.,\n",
       "       -9070., -9082., -9092., -9125., -9137., -9140., -9141., -9143.,\n",
       "       -9145., -9147., -9149., -9150., -9152., -9165., -9167., -9169.,\n",
       "       -9170., -9173., -9175., -9176., -9178., -9180., -9183., -9183.,\n",
       "       -9185., -9187., -9189., -9191., -9203., -9211., -9218., -9225.,\n",
       "       -9241., -9243., -9246., -9256., -9273., -9277., -9290., -9298.,\n",
       "       -9303., -9305., -9307., -9309., -9311., -9318., -9325., -9327.,\n",
       "       -9329., -9331., -9334., -9340., -9344., -9347., -9349., -9351.,\n",
       "       -9352., -9358., -9373., -9378., -9386., -9388., -9395., -9402.,\n",
       "       -9411., -9418., -9421., -9427., -9433., -9439., -9450., -9463.,\n",
       "       -9476., -9491., -9493., -9496.])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022.year.values[-100:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "81b13abb-10b0-4722-ac05-8b012ca7023a",
   "metadata": {},
   "outputs": [],
   "source": [
    "timevar = np.linspace(-9496+9496,1892+9496,11389)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "2a5fc4db-a617-4d46-920f-54ceeb7525ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11289., 11290., 11291., 11292., 11293., 11294., 11295., 11296.,\n",
       "       11297., 11298., 11299., 11300., 11301., 11302., 11303., 11304.,\n",
       "       11305., 11306., 11307., 11308., 11309., 11310., 11311., 11312.,\n",
       "       11313., 11314., 11315., 11316., 11317., 11318., 11319., 11320.,\n",
       "       11321., 11322., 11323., 11324., 11325., 11326., 11327., 11328.,\n",
       "       11329., 11330., 11331., 11332., 11333., 11334., 11335., 11336.,\n",
       "       11337., 11338., 11339., 11340., 11341., 11342., 11343., 11344.,\n",
       "       11345., 11346., 11347., 11348., 11349., 11350., 11351., 11352.,\n",
       "       11353., 11354., 11355., 11356., 11357., 11358., 11359., 11360.,\n",
       "       11361., 11362., 11363., 11364., 11365., 11366., 11367., 11368.,\n",
       "       11369., 11370., 11371., 11372., 11373., 11374., 11375., 11376.,\n",
       "       11377., 11378., 11379., 11380., 11381., 11382., 11383., 11384.,\n",
       "       11385., 11386., 11387., 11388.])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timevar[-100:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "77460731-2b0b-453e-858c-092e69e06d85",
   "metadata": {},
   "outputs": [],
   "source": [
    "contmass = (sigl2022.vssi*3).values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "3f2f6b2d-90e8-45a6-af6a-c58326ede575",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1384,)"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contmass.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "269b9b3b-715d-4c28-80bf-3d8d626e83e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1384,)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022.year.values.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "52f75559-cdb0-4625-b9c9-0c2f54809b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Your existing data\n",
    "sigl2022_years = sigl2022.year.values\n",
    "contmass = (sigl2022.vssi*3).values\n",
    "sigl_lats = sigl_hol_volc_data.lat\n",
    "\n",
    "# Your new time vector\n",
    "timevar = np.linspace(-9496+9496, 1892+9496, 11389)\n",
    "\n",
    "# Create a new array of zeros\n",
    "new_contmass = np.zeros_like(timevar)\n",
    "new_latitude = np.zeros_like(timevar)\n",
    "\n",
    "# Find the indices in timevar that correspond to sigl2022_years\n",
    "indices = np.searchsorted(timevar, sigl2022_years + 9496)\n",
    "\n",
    "# Fill in the values from contmass at those indices\n",
    "new_contmass[indices] = contmass\n",
    "new_latitude[indices] = sigl_lats\n",
    "\n",
    "# Now new_contmass is a continuous series with zeros filled in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "363a9fb7-ca3e-495d-b773-1739eb98a1f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Your existing data\n",
    "sigl2022_years = sigl2022.year.values\n",
    "contmass = (sigl2022.vssi*3).values\n",
    "sigl_lats = sigl_hol_volc_data.lat\n",
    "\n",
    "# Your new time vector\n",
    "timevar = np.linspace(-9496+9496, 1892+9496, 11389)\n",
    "\n",
    "# Create new arrays filled with np.nan\n",
    "new_contmass = np.full_like(timevar, np.nan)\n",
    "new_latitude = np.full_like(timevar, np.nan)\n",
    "\n",
    "# Find the indices in timevar that correspond to sigl2022_years\n",
    "indices = np.searchsorted(timevar, sigl2022_years + 9496)\n",
    "\n",
    "# Fill in the values from contmass and sigl_lats at those indices\n",
    "new_contmass[indices] = contmass\n",
    "new_latitude[indices] = sigl_lats\n",
    "\n",
    "# Now new_contmass and new_latitude are continuous series with np.nan filled in for missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "4a92fffd-b82e-41cd-8652-b4ed02af4151",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  8.67,    nan,    nan, 173.7 ,    nan,  27.39,    nan,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,    nan,    nan,    nan,\n",
       "          nan,    nan,    nan,    nan,   0.66,    nan,    nan,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,    nan,    nan,    nan,\n",
       "          nan,  16.95,    nan,    nan,    nan,    nan,    nan,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,    nan,   1.32,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,    nan,    nan,    nan,\n",
       "          nan,  20.25,    nan,    nan,    nan,    nan,    nan,  15.87,\n",
       "          nan,    nan,    nan,    nan,    nan,  15.69,    nan,    nan,\n",
       "          nan,    nan,    nan,  18.21,    nan,    nan,  33.66,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,   3.24,    nan,    nan,\n",
       "          nan,    nan,    nan,    nan,    nan,    nan,  37.71,    nan,\n",
       "          nan,    nan,    nan,    nan])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_contmass[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "ea632f2a-566c-40f7-a362-7bbe877257f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10., 11., 12.,\n",
       "       13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.,\n",
       "       26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38.,\n",
       "       39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51.,\n",
       "       52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62., 63., 64.,\n",
       "       65., 66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76., 77.,\n",
       "       78., 79., 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., 90.,\n",
       "       91., 92., 93., 94., 95., 96., 97., 98., 99.])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "timevar[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "759267e7-343d-4e8f-8d3b-d99be883f161",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 13.32,   3.12,  27.27,   0.84,  17.22,   0.72,   6.78,   4.35,\n",
       "        16.8 ,  53.25,   3.66, 110.85,  51.15,   3.18,  33.21,   7.89,\n",
       "        16.35,   5.22,  11.43,  12.81,   3.06,   3.45,  54.21,  17.34,\n",
       "        12.09,  13.86,   5.04,  66.15,   5.61,   5.4 ,  50.82,  28.86,\n",
       "         1.29,  32.7 ,   0.75,   6.39,   0.66,   0.69,   7.26,   2.85,\n",
       "         8.97,   2.46,   3.06,   2.49,  12.81,   6.06,  16.14,   5.22,\n",
       "       112.47,  32.31,   1.32,   1.95,   2.94,   1.2 ,   1.68,   8.01,\n",
       "         0.75,   8.16,  42.96,   5.4 ,   0.87,  28.83,   4.02,  43.26,\n",
       "        32.01,  19.62,   3.12, 225.93,   1.74,   0.63,  91.2 ,   0.72,\n",
       "         2.91,  12.84,   3.  ,  15.78,  24.87,   3.15,  35.7 ,   1.38,\n",
       "        56.76,   0.57,   3.12,   7.68,   3.57,   2.49,  57.54,   1.02,\n",
       "         1.47,   2.22,   0.48,   5.61,  24.09,   4.77,   1.14,  11.85,\n",
       "         4.95,   0.48,   0.48,   7.41])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contmass[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "659f25e9-6642-491f-bc3a-e62eb92bd9a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  1.29,  12.18,  43.71,   0.93,   2.52,  16.38,   1.38,   2.4 ,\n",
       "        13.89,   0.69,   5.37,  29.64,   0.96,   1.17, 251.01,   1.17,\n",
       "        25.29,  13.77,  49.59,  96.72,  16.92,   1.35,   4.65,  41.55,\n",
       "        20.94,  36.99,   1.08,  38.67,   0.87,  36.54,   0.87,   1.95,\n",
       "         2.1 ,   1.74,   2.46,   2.73,  59.58,   1.14,   1.02,   0.9 ,\n",
       "         1.11,  10.41,  13.62,   2.19,   1.62,   2.07,  61.08,   2.49,\n",
       "         2.04,  39.06,   1.11,   1.47,  70.47,  17.4 ,  61.26, 467.85,\n",
       "        23.19,   0.84,   1.35,  60.33,   0.87,   5.22,   1.17,   8.19,\n",
       "         2.16,   1.59,   2.49,   1.08,   2.4 ,  24.33,   2.28,   1.68,\n",
       "         1.95,   2.13,   2.19,   9.12,   2.82,   1.26,   1.83,   1.02,\n",
       "        20.07,   7.14,  11.79,   8.82,   0.57,  41.52, 128.91,  37.71,\n",
       "         3.24,  33.66,  18.21,  15.69,  15.87,  20.25,   1.32,  16.95,\n",
       "         0.66,  27.39, 173.7 ,   8.67])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contmass[-100:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "675c4989-0026-465a-87bc-912d4287989a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1892., 1886., 1884., 1876., 1863., 1861., 1855., 1844., 1835.,\n",
       "       1831., 1822., 1815., 1809., 1795., 1783., 1767., 1762., 1754.,\n",
       "       1739., 1729., 1719., 1702., 1694., 1673., 1666., 1661., 1646.,\n",
       "       1641., 1634., 1621., 1600., 1594., 1591., 1586., 1577., 1568.,\n",
       "       1559., 1555., 1540., 1527., 1512., 1505., 1495., 1490., 1480.,\n",
       "       1477., 1477., 1470., 1458., 1453., 1447., 1428., 1418., 1414.,\n",
       "       1390., 1385., 1381., 1378., 1344., 1342., 1336., 1330., 1290.,\n",
       "       1284., 1276., 1269., 1265., 1258., 1241., 1234., 1229., 1227.,\n",
       "       1214., 1209., 1205., 1199., 1191., 1188., 1181., 1181., 1171.,\n",
       "       1152., 1146., 1126., 1115., 1112., 1109., 1091., 1081., 1065.,\n",
       "       1044., 1038., 1029., 1022., 1015., 1008., 1006., 1004., 1002.,\n",
       "       1001.])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022_years[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "3fa8a5a5-20de-4f6c-98c2-7527a575ecb6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-8887., -8896., -8905., -8933., -8944., -8963., -8964., -8966.,\n",
       "       -8969., -8970., -8982., -8990., -8994., -9001., -9009., -9013.,\n",
       "       -9017., -9021., -9030., -9035., -9041., -9045., -9049., -9066.,\n",
       "       -9070., -9082., -9092., -9125., -9137., -9140., -9141., -9143.,\n",
       "       -9145., -9147., -9149., -9150., -9152., -9165., -9167., -9169.,\n",
       "       -9170., -9173., -9175., -9176., -9178., -9180., -9183., -9183.,\n",
       "       -9185., -9187., -9189., -9191., -9203., -9211., -9218., -9225.,\n",
       "       -9241., -9243., -9246., -9256., -9273., -9277., -9290., -9298.,\n",
       "       -9303., -9305., -9307., -9309., -9311., -9318., -9325., -9327.,\n",
       "       -9329., -9331., -9334., -9340., -9344., -9347., -9349., -9351.,\n",
       "       -9352., -9358., -9373., -9378., -9386., -9388., -9395., -9402.,\n",
       "       -9411., -9418., -9421., -9427., -9433., -9439., -9450., -9463.,\n",
       "       -9476., -9491., -9493., -9496.])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl2022_years[-100:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "782c4b4f-c355-456b-ab6c-cb6c764f5b1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 46156 characters]",
      "text/plain": [
       "<Figure size 2000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot chronologies\n",
    "\n",
    "SMALL_SIZE = 18\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 25\n",
    "#\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "\n",
    "#\n",
    "fig = plt.figure(figsize=(20,4))\n",
    "ax = plt.subplot(111)\n",
    "#\n",
    "for e in range(len(timevar)):\n",
    "    if new_latitude[e] > 20.0:\n",
    "        ax.plot([timevar[e], timevar[e]],[0, new_contmass[e]],color='b')\n",
    "    elif new_latitude[e] >= -20 and new_latitude[e] <= 20:\n",
    "        ax.plot([timevar[e], timevar[e]],[0, new_contmass[e]],color='k')\n",
    "    elif new_latitude[e] < -20.:\n",
    "        ax.plot([timevar[e], timevar[e]],[0, new_contmass[e]],color='r')\n",
    "#ax.set_xlim(24000,60000) \n",
    "ax.set_ylabel('Load (Tg SO$_4$)')\n",
    "#ax.set_ylim(0,180)\n",
    "ax.set_xlabel('Year B2K')\n",
    "#axes[axId].set_ylabel('Load (Tg)')\n",
    "#ax.set_yticks([0,25,50,75,100,125,150,175])\n",
    "\n",
    "\n",
    "#for e in range(len(year)):\n",
    "#    if flag[e] == 1:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='k')\n",
    "#    elif flag[e] == 2:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='b')\n",
    "#    elif flag[e] == 3:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='r')\n",
    "icepaperdir=\"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Timeseries_sigl_start0.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n",
    "    plt.savefig(icepaperdir+\"Timeseries_sigl_start0.pdf\", dpi=300, facecolor=fig.get_facecolor(), transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "b82fd1f7-e13f-49a1-a993-ac6255de2023",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No NaN values in the 'MaxLoad' column.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_287523/3129758532.py:21: FutureWarning: The default of observed=False is deprecated and will be changed to True in a future version of pandas. Pass observed=False to retain current behavior or observed=True to adopt the future default and silence this warning.\n",
      "  sigl_rebinned_df = df.groupby('Bin').agg(MaxLoad=('Load', 'max'))\n"
     ]
    }
   ],
   "source": [
    "## We can use this data set to rebin the data and do extreme value analysis return time with evenly spaced time intervals\n",
    "\n",
    "# Your actual data (this is just a placeholder for illustration, replace it with your actual dataset)\n",
    "time_series = timevar\n",
    "\n",
    "load_series = new_contmass\n",
    "\n",
    "# Convert the time_series and load_series into a DataFrame\n",
    "df = pd.DataFrame({'Time': time_series, 'Load': load_series})\n",
    "\n",
    "# Define custom bin edges based on time_series values\n",
    "# You can use a fixed number of bins or custom breakpoints.\n",
    "# Here, we use 5 bins for illustration\n",
    "num_bins = 242 # adjust higher until there are no more nan's\n",
    "bins = pd.cut(df['Time'], bins=num_bins)\n",
    "\n",
    "# Add the bin information to the DataFrame\n",
    "df['Bin'] = bins\n",
    "\n",
    "# Group by the 'Bin' and calculate the maximum 'Load' in each bin\n",
    "sigl_rebinned_df = df.groupby('Bin').agg(MaxLoad=('Load', 'max'))\n",
    "\n",
    "# Extract the bin edges (lower and upper bounds)\n",
    "sigl_rebinned_df['BinStart'] = sigl_rebinned_df.index.categories.left\n",
    "sigl_rebinned_df['BinEnd'] = sigl_rebinned_df.index.categories.right\n",
    "\n",
    "# Calculate the delta time for each bin\n",
    "sigl_rebinned_df['DeltaTime'] = sigl_rebinned_df['BinEnd'] - sigl_rebinned_df['BinStart']\n",
    "\n",
    "# Calculate the middle of the bin (BinMid)\n",
    "sigl_rebinned_df['BinMid'] = (sigl_rebinned_df['BinStart'] + sigl_rebinned_df['BinEnd']) / 2\n",
    "\n",
    "# Reorder the columns\n",
    "sigl_rebinned_df = sigl_rebinned_df[['BinStart', 'BinEnd', 'BinMid', 'DeltaTime', 'MaxLoad']]\n",
    "\n",
    "# Output the rebinned dataset\n",
    "#print(sigl_rebinned_df)\n",
    "#print(np.min(sigl_rebinned_df[\"MaxLoad\"]))\n",
    "#print(np.max(sigl_rebinned_df[\"MaxLoad\"]))\n",
    "\n",
    "# Assuming sigl_rebinned_df is your DataFrame\n",
    "if sigl_rebinned_df['MaxLoad'].isna().any():\n",
    "    print(\"There are NaN values in the 'MaxLoad' column.\")\n",
    "else:\n",
    "    print(\"No NaN values in the 'MaxLoad' column.\")\n",
    "sigl_new_rebinned_df  = sigl_rebinned_df[1:]    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "2769583f-5239-40e2-a429-7dca5d024bcf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>BinStart</th>\n",
       "      <th>BinEnd</th>\n",
       "      <th>BinMid</th>\n",
       "      <th>DeltaTime</th>\n",
       "      <th>MaxLoad</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Bin</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>(47.058, 94.116]</th>\n",
       "      <td>47.058</td>\n",
       "      <td>94.116</td>\n",
       "      <td>70.5870</td>\n",
       "      <td>47.058</td>\n",
       "      <td>37.71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(94.116, 141.174]</th>\n",
       "      <td>94.116</td>\n",
       "      <td>141.174</td>\n",
       "      <td>117.6450</td>\n",
       "      <td>47.058</td>\n",
       "      <td>128.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(141.174, 188.231]</th>\n",
       "      <td>141.174</td>\n",
       "      <td>188.231</td>\n",
       "      <td>164.7025</td>\n",
       "      <td>47.057</td>\n",
       "      <td>24.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(188.231, 235.289]</th>\n",
       "      <td>188.231</td>\n",
       "      <td>235.289</td>\n",
       "      <td>211.7600</td>\n",
       "      <td>47.058</td>\n",
       "      <td>8.19</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(235.289, 282.347]</th>\n",
       "      <td>235.289</td>\n",
       "      <td>282.347</td>\n",
       "      <td>258.8180</td>\n",
       "      <td>47.058</td>\n",
       "      <td>467.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(11152.711, 11199.769]</th>\n",
       "      <td>11152.711</td>\n",
       "      <td>11199.769</td>\n",
       "      <td>11176.2400</td>\n",
       "      <td>47.058</td>\n",
       "      <td>54.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(11199.769, 11246.826]</th>\n",
       "      <td>11199.769</td>\n",
       "      <td>11246.826</td>\n",
       "      <td>11223.2975</td>\n",
       "      <td>47.057</td>\n",
       "      <td>12.81</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(11246.826, 11293.884]</th>\n",
       "      <td>11246.826</td>\n",
       "      <td>11293.884</td>\n",
       "      <td>11270.3550</td>\n",
       "      <td>47.058</td>\n",
       "      <td>33.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(11293.884, 11340.942]</th>\n",
       "      <td>11293.884</td>\n",
       "      <td>11340.942</td>\n",
       "      <td>11317.4130</td>\n",
       "      <td>47.058</td>\n",
       "      <td>110.85</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>(11340.942, 11388.0]</th>\n",
       "      <td>11340.942</td>\n",
       "      <td>11388.000</td>\n",
       "      <td>11364.4710</td>\n",
       "      <td>47.058</td>\n",
       "      <td>27.27</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>241 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                         BinStart     BinEnd      BinMid  DeltaTime  MaxLoad\n",
       "Bin                                                                         \n",
       "(47.058, 94.116]           47.058     94.116     70.5870     47.058    37.71\n",
       "(94.116, 141.174]          94.116    141.174    117.6450     47.058   128.91\n",
       "(141.174, 188.231]        141.174    188.231    164.7025     47.057    24.33\n",
       "(188.231, 235.289]        188.231    235.289    211.7600     47.058     8.19\n",
       "(235.289, 282.347]        235.289    282.347    258.8180     47.058   467.85\n",
       "...                           ...        ...         ...        ...      ...\n",
       "(11152.711, 11199.769]  11152.711  11199.769  11176.2400     47.058    54.21\n",
       "(11199.769, 11246.826]  11199.769  11246.826  11223.2975     47.057    12.81\n",
       "(11246.826, 11293.884]  11246.826  11293.884  11270.3550     47.058    33.21\n",
       "(11293.884, 11340.942]  11293.884  11340.942  11317.4130     47.058   110.85\n",
       "(11340.942, 11388.0]    11340.942  11388.000  11364.4710     47.058    27.27\n",
       "\n",
       "[241 rows x 5 columns]"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigl_new_rebinned_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "9ecb78fa-beb9-488e-bfaf-d3271254a9a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[<matplotlib.lines.Line2D at 0x7ee802a88ef0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bb9be0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bb8500>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bb8f50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bbaf90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bb93a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c53080>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c50ce0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c52cf0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c50590>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c51820>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c51a00>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c50fb0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c50f80>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c53500>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c51730>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c52d20>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c534d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c50d40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c51040>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c52b70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c519a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b023c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b000e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b01850>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b00290>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b00d40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b02c90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b00680>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b01640>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805669940>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80dac7170>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8046a9070>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802a8a720>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8045732f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804a385f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802966210>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b01ac0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b008c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b03ef0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b02720>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b01d90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b03d40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b03380>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b03140>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b006b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b038c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b01af0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b02e70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802b021e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c718e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c71820>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c70ef0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c70da0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802c72870>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6ba0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6fc0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6e70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6780>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba4e90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6f30>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6300>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba73b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba4a70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba6810>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba69c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba7c80>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802ba75c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bf90d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bf9be0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bf9f70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfafc0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfb320>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfb560>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfad50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bf9e20>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfb530>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bf9d30>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfa450>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfbf50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfbef0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfb620>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee802bfa000>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f926c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f931a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f92330>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f93aa0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f91fa0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f90da0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f90b60>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f91850>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f91280>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f92d50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f936b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f91580>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f90f50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804f902f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96cda0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96c0e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96cfe0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96da60>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96cb00>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96d760>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96f800>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96ea80>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96cb60>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80d96c590>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051eec90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ec080>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ef3b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051efcb0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051edc40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ec290>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ef860>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ecad0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ed8e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051eed20>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051ef200>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8051efb30>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522d2b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e3c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f380>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522c0e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522c5c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e660>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522c920>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522c3b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522cb30>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522d340>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522dd90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e690>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f1d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f6e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522ff20>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f170>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e810>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f770>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522c1a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e180>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522f5c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80522e750>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5f0e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5f320>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5e9f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5dc40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5c530>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5c0b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5c590>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5c950>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5d220>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5d730>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5db50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5e2a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5eb70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5f140>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5f2c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5fd40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5fe90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5c860>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5d760>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5f9e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5ec00>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5ed80>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804e5ddc0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518ca10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518d460>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518e5d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518ee70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518c380>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518d100>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518f980>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518d970>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518c410>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518cce0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518d850>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518e2a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518e990>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518f230>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518fcb0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518edb0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518e780>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518cd10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518c170>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518df10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518f740>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee80518fe60>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff1040>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3380>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff22a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff11f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff0560>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff07a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff0f80>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff0b90>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff1cd0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff1fd0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff29f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff2ed0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff34d0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3770>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3f50>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff0470>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff1940>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3260>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff15b0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3d70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff1dc0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff3980>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804ff2870>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054435c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805442420>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805442030>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441280>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054402c0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805440710>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805440ec0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441760>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441c40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805442240>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805442120>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805442e70>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054434a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805443980>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805443b60>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054437a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054404a0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee8054439e0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441dc0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441160>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805443260>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805443470>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee805441850>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc3830>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc1c40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc1f10>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc1370>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc0380>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc0440>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc0d40>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc14f0>],\n",
       " [<matplotlib.lines.Line2D at 0x7ee804fc1c70>]]"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 21212 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[ plt.plot([sigl_new_rebinned_df[\"BinMid\"].values[e],sigl_new_rebinned_df[\"BinMid\"].values[e]], [0, sigl_new_rebinned_df[\"MaxLoad\"].values[e]],color='k') \n",
    " for e in range(len(sigl_new_rebinned_df[\"BinMid\"])) ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "3a719f3c-25f0-4a09-bdd1-af340232c3bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.551963396388254\n",
      "22.867942330855666\n",
      "20.995681539070155\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 49744 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set up the plot\n",
    "fig, ax = plt.subplots(figsize=(14, 9))\n",
    "\n",
    "# Plot histogram for Lin data\n",
    "hist_lin, bins_lin = np.histogram(sigl_new_rebinned_df[\"MaxLoad\"], bins=50, range=(10, 600), density=True)\n",
    "ax.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"Lin data\")\n",
    "\n",
    "#ax.set_yscale('log')\n",
    "ax.set_ylabel(\"Normalized Occurrences\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\", size=16)\n",
    "setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Normalized Occurrences\", 16)\n",
    "\n",
    "#ax.set_title('Levy-stable distribution fit for Sigl and Lin data', size=16)\n",
    "ax.set_ylim(0, 0.01)\n",
    "plt.legend(fontsize=12)\n",
    "\n",
    "def main(rvs):\n",
    "    shape, loc, scale = gev.fit(rvs)\n",
    "    return shape, loc, scale\n",
    "\n",
    "\n",
    "shape, loc, scale = main(sigl_new_rebinned_df[\"MaxLoad\"].values)\n",
    "\n",
    "print(shape)\n",
    "print(loc)\n",
    "print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "#hist, bins = np.histogram(sigl_new_rebinned_df[\"MaxLoad\"].values, bins=5, range=(0, 600), density=True)\n",
    "#plt.bar(bins[:-1], hist, width = 2, align='edge')\n",
    "\n",
    "plt.plot(xx, yy, 'ro')\n",
    "\n",
    "# these are the GEV parameters for the massaged Lin data\n",
    "sigl_ml_loc = loc \n",
    "sigl_ml_scale = scale\n",
    "sigl_ml_shape = shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "5d1f4fd7-8f2c-4e17-8dba-0b7b86fd1008",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 45300 characters]",
      "text/plain": [
       "<Figure size 2000x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot chronologies\n",
    "# add offset to get it with correct dates before 2K data (data is from 1892)\n",
    "t_shift = 2000 - 1892\n",
    "t_start = -11389 + 1892\n",
    "SMALL_SIZE = 18\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 25\n",
    "#\n",
    "plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "\n",
    "#\n",
    "fig = plt.figure(figsize=(20,4))\n",
    "ax = plt.subplot(111)\n",
    "#\n",
    "for e in range(len(timevar)):\n",
    "    if new_latitude[e] > 20.0:\n",
    "        ax.plot([t_start + timevar[e], t_start + timevar[e]],[0, new_contmass[e]],color='b')\n",
    "    elif new_latitude[e] >= -20 and new_latitude[e] <= 20:\n",
    "        ax.plot([t_start + timevar[e], t_start + timevar[e]],[0, new_contmass[e]],color='k')\n",
    "    elif new_latitude[e] < -20.:\n",
    "        ax.plot([t_start + timevar[e], t_start + timevar[e]],[0, new_contmass[e]],color='r')\n",
    "#ax.set_xlim(24000,60000) \n",
    "ax.set_ylabel('Load (Tg SO$_4$)')\n",
    "#ax.set_ylim(0,180)\n",
    "ax.set_xlabel('Year CE')\n",
    "#axes[axId].set_ylabel('Load (Tg)')\n",
    "#ax.set_yticks([0,25,50,75,100,125,150,175])\n",
    "\n",
    "\n",
    "#for e in range(len(year)):\n",
    "#    if flag[e] == 1:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='k')\n",
    "#    elif flag[e] == 2:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='b')\n",
    "#    elif flag[e] == 3:\n",
    "#        ax.plot([year[e],year[e]],[0,mass[e]],color='r')\n",
    "icepaperdir=\"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot=False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Timeseries_sigl_start0.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=True)\n",
    "    plt.savefig(icepaperdir+\"Timeseries_sigl_start0.pdf\", dpi=300, facecolor=fig.get_facecolor(), transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "0ce86e98-1a22-4777-93b9-235dbf6f0185",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted parameters:\n",
      "Shape: -1.34426, Location: 2.49871, Scale: 3.16342\n"
     ]
    }
   ],
   "source": [
    "shape, loc, scale = gev.fit(mass)\n",
    "print(f\"Fitted parameters:\\nShape: {shape:.5f}, Location: {loc:.5f}, Scale: {scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "c34dbf3d-ae40-4047-918d-bcde4952c450",
   "metadata": {},
   "outputs": [],
   "source": [
    "#shape, loc, scale = gev.fit(new_contmass)\n",
    "#print(f\"Fitted parameters:\\nShape: {shape:.5f}, Location: {loc:.5f}, Scale: {scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "ae6c64e8-14ab-41ce-848e-046435b4fa12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1454284983635481\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ee8020da3f0>"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 32220 characters]",
      "text/plain": [
       "<Figure size 1400x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(14,9))\n",
    "ax = plt.subplot(111)\n",
    "#shape, loc, scale = gev.fit(mass)\n",
    "\n",
    "#print(shape)\n",
    "#print(loc)\n",
    "#print(scale)\n",
    "\n",
    "l = loc + scale / shape\n",
    "print (l)\n",
    "xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "yy = gev.pdf(xx, shape, loc, scale)\n",
    "\n",
    "hist, bins = np.histogram(new_contmass, bins=31, range=(0, 600), density=True)\n",
    "ax.bar(bins[:-1], hist, width = 10, align='edge', label=\"Occurance #\")\n",
    "\n",
    "# plot GEV dist\n",
    "#ax.plot(xx, yy, 'orange', linewidth = 3.5, label = \"GEV fit\")\n",
    "# set log scale\n",
    "ax.set_yscale('log')\n",
    "#ax.set_ylabel(\"Log Normalized Occurances\")\n",
    "ax.set_xlabel(\"Stratospheric Load (Tg SO$_4$)\",size=16)\n",
    "#setaxesfsize(ax, 14, \"Stratospheric Load (Tg SO$_4$)\", \"Log Normalized Occurances\",  16)\n",
    "fshape = np.around(shape,2)\n",
    "floc = np.around(loc,2)\n",
    "fscale = np.around(scale,2)\n",
    "#ax.set_title(r'$GEV(x) = exp[-\\{1+s(x − a)/b\\}^{-1/s}]$  s='+str(fshape)+'  a='+str(floc)+'  b='+str(fscale),size=16)\n",
    "#ax.set_title(r'$alpha(x) = exp[-\\{1+s(x − a)/b\\}^{-1/s}]$  s='+str(fshape)+'  a='+str(floc)+'  b='+str(fscale),size=16)\n",
    "#ax.set_xlim(-10,100)\n",
    "#ax.set_ylim(1.e-3,1)\n",
    "plt.legend(fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "6c78eedf-9ae1-44ec-8124-6d062d65b497",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "np.float64(3.997063688329301)"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimate_return_level(0.5, loc, scale, shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "b6825039-b697-412d-b498-f9c44eb48eab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean annual SO2 deposition: 1.9357 Tg\n"
     ]
    }
   ],
   "source": [
    "# Sum of all SO2 depositions\n",
    "total_so2 = np.nansum(new_contmass)\n",
    "\n",
    "# Total number of years\n",
    "total_years = len(timevar)\n",
    "\n",
    "# Calculate mean annual deposition\n",
    "mean_annual_so2 = total_so2 / total_years\n",
    "\n",
    "print(f\"Mean annual SO2 deposition: {mean_annual_so2:.4f} Tg\")"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8962dc04-e411-44a6-90d4-965a3a6ddc47",
   "metadata": {},
   "source": [
    "\n",
    "# setup plots\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# get empirical return levels and plot them\n",
    "_ = empirical_return_level(precipitation).plot(ax=ax, marker=\".\", linestyle=\"None\")\n",
    "\n",
    "# create vector of years\n",
    "years = np.arange(1.1, 100, 0.1)\n",
    "\n",
    "# calculate and plot the normal return levels\n",
    "ax.plot(\n",
    "    years,\n",
    "    stats.norm.ppf(1 - 1 / years, loc=precipitation.mean(), scale=precipitation.std()),\n",
    ")\n",
    "\n",
    "# calculate and plot the GEV distribution, note the negtive shape parameter\n",
    "ax.plot(years, gev.ppf(1 - 1 / years, shape, loc=loc, scale=scale))\n",
    "# set x axis to log scale\n",
    "ax.set_xscale(\"log\")\n",
    "\n",
    "# show legend\n",
    "ax.legend([\"empirical\", \"normal\", \"GEV\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "8f9da197-6971-411b-92de-b887d85081a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_287523/2369074794.py:7: RuntimeWarning: divide by zero encountered in log\n",
      "  yp = -np.log(1 - p)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 71636 characters]",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Parameters\n",
    "mu = 0\n",
    "sigma = 1\n",
    "\n",
    "# Generate points\n",
    "p = np.linspace(0, 1, 1000)[1:]  # ppoints(1000) without the first point\n",
    "yp = -np.log(1 - p)\n",
    "\n",
    "# Calculate zp values\n",
    "zp_1 = mu - (sigma / 0.2) * (1 - yp ** (-0.2))   # xi = 0.2\n",
    "zp_2 = mu - sigma * np.log(yp)                    # xi = 0\n",
    "zp_3 = mu - (sigma / -0.2) * (1 - yp ** (0.2))   # xi = -0.2\n",
    "\n",
    "# Plotting\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(1 / p, zp_1, label='xi = 0.2', color='blue')\n",
    "plt.plot(1 / p, zp_2, label='xi = 0', color='orange')\n",
    "plt.plot(1 / p, zp_3, label='xi = -0.2', color='green')\n",
    "\n",
    "# Log scale and labels\n",
    "plt.xscale('log')\n",
    "plt.xlabel('time period')\n",
    "plt.ylabel('return level')\n",
    "plt.title('Return level versus time period (1/p)')\n",
    "#plt.xticks([10**i for i in range(5)], [f'10^{i}' for i in range(5)])\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "23161e53-db68-4c8c-94f8-152b4200b5e1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.229046242774567"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## because the data is not on an equally spaced yearly interval we need to calculate a scale time factor for return times:\n",
    "sigl_tscale = 11389./1384.\n",
    "sigl_tscale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "8e3665a0-55ef-4286-8c2d-9b8967c2b6de",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "538.4615384615385"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## because the data is not on an equally spaced yearly interval we need to calculate a scale time factor for return times:\n",
    "lin_len = len(mis3_lin_volc_load)\n",
    "print(lin_len)\n",
    "lin_tscale = 35000./65\n",
    "lin_tscale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "ffd4d869-8f9c-4ec6-a4b1-bc96b5fbd08b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_return_level(quantile, loc, scale, shape):\n",
    "    level = loc + scale / shape * (1 - (-np.log(quantile)) ** (shape))\n",
    "    return level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "8afe38e4-e43d-4abd-89cc-e67d9d1726a4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.344260778639558\n",
      "-0.4685136930628534\n"
     ]
    }
   ],
   "source": [
    "# the scipy shape parameter is the negative of \\xi \"c\" = \"-\\xi\"\n",
    "print(sigl_shape)\n",
    "print(lin_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "90416bdf-7450-4e0b-abdc-c4e8233c1f38",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_287523/3205516377.py:14: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored.\n",
      "  ax.set_xlim(0,100000)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 48672 characters]",
      "text/plain": [
       "<Figure size 1400x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "periods = np.arange(2, 1000, 2)\n",
    "quantiles = 1 - 1 / periods\n",
    "sigl_levels = estimate_return_level(quantiles, sigl_loc, sigl_scale, -1.0*sigl_shape)\n",
    "lin_levels = estimate_return_level(quantiles, lin_loc, lin_scale, -1.0*lin_shape)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(14,8))\n",
    "ax.plot(periods*sigl_tscale, sigl_levels, \"-\", color='red')\n",
    "ax.plot(periods*lin_tscale, lin_levels, \"-\", color='blue')\n",
    "ax.set_xlabel(\"Return Period (years)\")\n",
    "ax.set_ylabel(\"Return Level (Tg SO$_4$)\")\n",
    "ax.set_xscale(\"log\")\n",
    "#ax.set_yscale(\"log\")\n",
    "#ax.set_ylim(0,600)\n",
    "ax.set_xlim(0,100000)\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "2beead24-5eca-46c7-aeda-a6f5b05b676a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def empirical_return_level(data):\n",
    "    \"\"\"\n",
    "    Compute empirical return level using the algorithm introduced in Tutorial 2\n",
    "    \"\"\"\n",
    "    df = pd.DataFrame(index=np.arange(data.size))\n",
    "    # sort the data\n",
    "    df[\"sorted\"] = np.sort(data)[::-1]\n",
    "    # rank via scipy instead to deal with duplicate values\n",
    "    df[\"ranks_sp\"] = np.sort(stats.rankdata(-data))\n",
    "    # find exceedence probability\n",
    "    n = data.size\n",
    "    df[\"exceedance\"] = df[\"ranks_sp\"] / (n + 1)\n",
    "    # find return period\n",
    "    df[\"period\"] = 1 / df[\"exceedance\"]\n",
    "\n",
    "    df = df[::-1]\n",
    "\n",
    "    out = xr.DataArray(\n",
    "        dims=[\"period\"],\n",
    "        coords={\"period\": df[\"period\"]},\n",
    "        data=df[\"sorted\"],\n",
    "        name=\"level\",\n",
    "    )\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "87452c77-33b6-4970-a907-9c3c48f5d20f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36196 characters]",
      "text/plain": [
       "<Figure size 1400x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "periods = np.arange(2, 1000, 2)\n",
    "quantiles = 1 - 1 / periods\n",
    "levels = estimate_return_level(quantiles, loc, scale, -shape)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(14,8))\n",
    "ax.plot(periods, levels, \"-\")\n",
    "ax.set_xlabel(\"Return Period (years)\")\n",
    "ax.set_ylabel(\"Return Level (Tg SO$_4$)\")\n",
    "#ax.set_xscale(\"log\")\n",
    "#ax.set_yscale(\"log\")\n",
    "ax.set_ylim(0,600)\n",
    "ax.set_xlim(0,1500)\n",
    "ax.grid(True)\n",
    "\n",
    "# get empirical return levels and plot them\n",
    "# sigl\n",
    "_ = (empirical_return_level(mass)/sigl_tscale).plot(ax=ax, marker=\".\", color='blue', linestyle=\"None\")\n",
    "# lin\n",
    "#_ = empirical_return_level(volc_load).plot(ax=ax, marker=\".\", color='red', linestyle=\"None\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "93b7f606-e486-42c2-a253-2b677426ba67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in notebooks */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base rgba(0, 0, 0, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(0, 0, 0, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, white)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, white) h s calc(l - 15))\n",
       "  );\n",
       "}\n",
       "\n",
       "html[theme=\"dark\"],\n",
       "html[data-theme=\"dark\"],\n",
       "body[data-theme=\"dark\"],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: var(\n",
       "    --jp-content-font-color0,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 1))\n",
       "  );\n",
       "  --xr-font-color2: var(\n",
       "    --jp-content-font-color2,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.54))\n",
       "  );\n",
       "  --xr-font-color3: var(\n",
       "    --jp-content-font-color3,\n",
       "    var(--pst-color-text-base, rgba(255, 255, 255, 0.38))\n",
       "  );\n",
       "  --xr-border-color: var(\n",
       "    --jp-border-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 10))\n",
       "  );\n",
       "  --xr-disabled-color: var(\n",
       "    --jp-layout-color3,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 40))\n",
       "  );\n",
       "  --xr-background-color: var(\n",
       "    --jp-layout-color0,\n",
       "    var(--pst-color-on-background, #111111)\n",
       "  );\n",
       "  --xr-background-color-row-even: var(\n",
       "    --jp-layout-color1,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 5))\n",
       "  );\n",
       "  --xr-background-color-row-odd: var(\n",
       "    --jp-layout-color2,\n",
       "    hsl(from var(--pst-color-on-background, #111111) h s calc(l + 15))\n",
       "  );\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block !important;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "  line-height: 1.6;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-obj-name,\n",
       ".xr-group-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-group-name::before {\n",
       "  content: \"📁\";\n",
       "  padding-right: 0.3em;\n",
       "}\n",
       "\n",
       ".xr-group-name,\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 0 20px 0 20px;\n",
       "  margin-block-start: 0;\n",
       "  margin-block-end: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: inline-block;\n",
       "  opacity: 0;\n",
       "  height: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "  border: 2px solid transparent !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:focus + label {\n",
       "  border: 2px solid var(--xr-font-color0) !important;\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: \"►\";\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: \"▼\";\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-top: 4px;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-group-box {\n",
       "  display: inline-grid;\n",
       "  grid-template-columns: 0px 20px auto;\n",
       "  width: 100%;\n",
       "}\n",
       "\n",
       ".xr-group-box-vline {\n",
       "  grid-column-start: 1;\n",
       "  border-right: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "  width: 0px;\n",
       "}\n",
       "\n",
       ".xr-group-box-hline {\n",
       "  grid-column-start: 2;\n",
       "  grid-row-start: 1;\n",
       "  height: 1em;\n",
       "  width: 20px;\n",
       "  border-bottom: 0.2em solid;\n",
       "  border-color: var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-group-box-contents {\n",
       "  grid-column-start: 3;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: \"(\";\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: \")\";\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: \",\";\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  border-color: var(--xr-background-color-row-odd);\n",
       "  margin-bottom: 0;\n",
       "  padding-top: 2px;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "  border-color: var(--xr-background-color-row-even);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-index-preview {\n",
       "  grid-column: 2 / 5;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  display: none;\n",
       "  border-top: 2px dotted var(--xr-background-color);\n",
       "  padding-bottom: 20px !important;\n",
       "  padding-top: 10px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in + label,\n",
       ".xr-var-data-in + label,\n",
       ".xr-index-data-in + label {\n",
       "  padding: 0 1px;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data,\n",
       ".xr-index-data-in:checked ~ .xr-index-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-data > pre,\n",
       ".xr-index-data > pre,\n",
       ".xr-var-data > table > tbody > tr {\n",
       "  background-color: transparent !important;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-index-name div,\n",
       ".xr-index-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data,\n",
       ".xr-index-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2,\n",
       ".xr-no-icon {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked + label > .xr-icon-file-text2,\n",
       ".xr-var-data-in:checked + label > .xr-icon-database,\n",
       ".xr-index-data-in:checked + label > .xr-icon-database {\n",
       "  color: var(--xr-font-color0);\n",
       "  filter: drop-shadow(1px 1px 5px var(--xr-font-color2));\n",
       "  stroke-width: 0.8px;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray &#x27;level&#x27; (period: 11389)&gt; Size: 91kB\n",
       "array([0.39, 0.39, 0.42, ...,  nan,  nan,  nan], shape=(11389,))\n",
       "Coordinates:\n",
       "  * period   (period) float64 91kB nan nan nan nan nan ... nan nan nan nan nan</pre><div class='xr-wrap' style='display:none'><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-obj-name'>&#x27;level&#x27;</div><ul class='xr-dim-list'><li><span class='xr-has-index'>period</span>: 11389</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-fb2b80d7-670c-4a39-8b5e-7b1e252a750d' class='xr-array-in' type='checkbox' checked><label for='section-fb2b80d7-670c-4a39-8b5e-7b1e252a750d' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>0.39 0.39 0.42 0.42 0.42 0.45 0.48 ... nan nan nan nan nan nan nan</span></div><div class='xr-array-data'><pre>array([0.39, 0.39, 0.42, ...,  nan,  nan,  nan], shape=(11389,))</pre></div></div></li><li class='xr-section-item'><input id='section-7f95d190-635f-46a1-bc91-e44832147ef5' class='xr-section-summary-in' type='checkbox'  checked><label for='section-7f95d190-635f-46a1-bc91-e44832147ef5' class='xr-section-summary' >Coordinates: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>period</span></div><div class='xr-var-dims'>(period)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>nan nan nan nan ... nan nan nan nan</div><input id='attrs-abe90487-511f-42b7-8d80-5e08994d20f5' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-abe90487-511f-42b7-8d80-5e08994d20f5' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-e6769c41-4ad5-4828-b50a-8bece37107c9' class='xr-var-data-in' type='checkbox'><label for='data-e6769c41-4ad5-4828-b50a-8bece37107c9' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([nan, nan, nan, ..., nan, nan, nan], shape=(11389,))</pre></div></li></ul></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.DataArray 'level' (period: 11389)> Size: 91kB\n",
       "array([0.39, 0.39, 0.42, ...,  nan,  nan,  nan], shape=(11389,))\n",
       "Coordinates:\n",
       "  * period   (period) float64 91kB nan nan nan nan nan ... nan nan nan nan nan"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " empirical_return_level(new_contmass)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b3badfb2-7d2f-4932-b0ff-ed1a4e75917e",
   "metadata": {},
   "source": [
    "# initalize list to store parameters from samples\n",
    "params = []\n",
    "\n",
    "# generate 1000 samples by resampling data with replacement\n",
    "for i in range(1000):\n",
    "    params.append(\n",
    "        gev.fit(np.random.choice(precipitation, size=precipitation.size, replace=True))\n",
    "    )\n",
    "\n",
    "# print the estimate of the mean of each parameter and it's confidence intervals\n",
    "print(\n",
    "    \"Mean estimate: \",\n",
    "    np.mean(np.array(params), axis=0),\n",
    "    \" and 95% confidence intervals: \",\n",
    "    np.quantile(np.array(params), [0.025, 0.975], axis=0),\n",
    ")\n",
    "\n",
    "# generate years vector\n",
    "years = np.arange(1.1, 1000, 0.1)\n",
    "\n",
    "# intialize list for return levels\n",
    "levels = []\n",
    "\n",
    "# calculate return levels for each of the 1000 samples\n",
    "for i in range(1000):\n",
    "    levels.append(gev.ppf(1 - 1 / years, *params[i]))\n",
    "levels = np.array(levels)\n",
    "\n",
    "# setup plots\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# find empirical return levels\n",
    "_ = empirical_return_level(precipitation).plot(ax=ax, marker=\".\", linestyle=\"None\")\n",
    "\n",
    "# plot return mean levels\n",
    "_ = ax.plot(years, levels.mean(axis=0))\n",
    "\n",
    "# plot confidence intervals\n",
    "_ = ax.plot(years, np.quantile(levels, [0.025, 0.975], axis=0).T, \"k--\")\n",
    "\n",
    "# aesthetics\n",
    "ax.set_xlim(1.5, 1000)\n",
    "ax.set_ylim(20, 110)\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_xlabel(\"Return Period (years)\")\n",
    "ax.set_ylabel(\"Return Level (mm/day)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "9ff605fd-cabd-44cf-9797-be38d6822ac2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32  34  36\n",
      "  38  40  42  44  46  48  50  52  54  56  58  60  62  64  66  68  70  72\n",
      "  74  76  78  80  82  84  86  88  90  92  94  96  98 100 102 104 106 108\n",
      " 110 112 114 116 118 120 122 124 126 128 130 132 134 136 138 140 142 144\n",
      " 146 148 150 152 154 156 158 160 162 164 166 168 170 172 174 176 178 180\n",
      " 182 184 186 188 190 192 194 196 198 200 202 204 206 208 210 212 214 216\n",
      " 218 220 222 224 226 228 230 232 234 236 238 240 242 244 246 248 250 252\n",
      " 254 256 258 260 262 264 266 268 270 272 274 276 278 280 282 284 286 288\n",
      " 290 292 294 296 298 300 302 304 306 308 310 312 314 316 318 320 322 324\n",
      " 326 328 330 332 334 336 338 340 342 344 346 348 350 352 354 356 358 360\n",
      " 362 364 366 368 370 372 374 376 378 380 382 384 386 388 390 392 394 396\n",
      " 398 400 402 404 406 408 410 412 414 416 418 420 422 424 426 428 430 432\n",
      " 434 436 438 440 442 444 446 448 450 452 454 456 458 460 462 464 466 468\n",
      " 470 472 474 476 478 480 482 484 486 488 490 492 494 496 498 500 502 504\n",
      " 506 508 510 512 514 516 518 520 522 524 526 528 530 532 534 536 538 540\n",
      " 542 544 546 548 550 552 554 556 558 560 562 564 566 568 570 572 574 576\n",
      " 578 580 582 584 586 588 590 592 594 596 598 600 602 604 606 608 610 612\n",
      " 614 616 618 620 622 624 626 628 630 632 634 636 638 640 642 644 646 648\n",
      " 650 652 654 656 658 660 662 664 666 668 670 672 674 676 678 680 682 684\n",
      " 686 688 690 692 694 696 698 700 702 704 706 708 710 712 714 716 718 720\n",
      " 722 724 726 728 730 732 734 736 738 740 742 744 746 748 750 752 754 756\n",
      " 758 760 762 764 766 768 770 772 774 776 778 780 782 784 786 788 790 792\n",
      " 794 796 798 800 802 804 806 808 810 812 814 816 818 820 822 824 826 828\n",
      " 830 832 834 836 838 840 842 844 846 848 850 852 854 856 858 860 862 864\n",
      " 866 868 870 872 874 876 878 880 882 884 886 888 890 892 894 896 898 900\n",
      " 902 904 906 908 910 912 914 916 918 920 922 924 926 928 930 932 934 936\n",
      " 938 940 942 944 946 948 950 952 954 956 958 960 962 964 966 968 970 972\n",
      " 974 976 978 980 982 984 986 988 990 992 994 996 998]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7ee802965580>"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 67420 characters]",
      "text/plain": [
       "<Figure size 1400x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#periods = np.arange(2, 1000, 1)\n",
    "sigl_periods = np.linspace(1.1, 130, 1000)\n",
    "\n",
    "lin_periods = np.linspace(1.1, 13, 100)\n",
    "print(periods)\n",
    "sigl_quantiles = 1 - 1 / sigl_periods\n",
    "lin_quantiles = 1 - 1 / lin_periods\n",
    "\n",
    "sigl_tscale = 47.06 # years in each maximum bin\n",
    "sigl_levels = estimate_return_level(sigl_quantiles, sigl_ml_loc, sigl_ml_scale, sigl_ml_shape)\n",
    "#lin_ml_loc\n",
    "lin_levels = estimate_return_level(lin_quantiles, lin_ml_loc, lin_ml_scale, lin_ml_shape)\n",
    "lin_tscale = 1520.8 # years in each maximum bin\n",
    "fig, ax = plt.subplots(figsize=(14,8))\n",
    "ax.plot(sigl_periods*sigl_tscale, sigl_levels, \"-\", color='red', label='Sigl Data')\n",
    "ax.plot(lin_periods*lin_tscale, lin_levels, \"-\", color='blue', label='Lin Data')\n",
    "ax.set_xlabel(\"Return Period (years)\")\n",
    "ax.set_ylabel(\"Return Level (Tg SO$_4$)\")\n",
    "#ax.set_xscale(\"log\")\n",
    "#ax.set_yscale(\"log\")\n",
    "ax.set_ylim(0,500)\n",
    "ax.set_xlim(0,10000)\n",
    "ax.grid(True)\n",
    "setaxesfsize(ax,14,\"Return Period (Years)\",\"Return Level (Tg SO$_4$)\",16)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "0c28f63a-9f6b-4f5e-a06b-535ac710463a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14920 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.linspace(0, len(mis3_lin_volc_load), len(mis3_lin_volc_load))*lin_tscale, empirical_return_level(mis3_lin_volc_load), marker=\".\", color='red', linestyle=\"None\")\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "b62aaa59-c0d1-4186-b4e1-0f07b38e8c71",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14920 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.linspace(0, len(mis3_lin_volc_load), len(mis3_lin_volc_load))*lin_tscale, empirical_return_level(mis3_lin_volc_load), marker=\".\", color='red', linestyle=\"None\")\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "b0e8a365-dda5-4a54-9953-d732cdf82dea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 27764 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = empirical_return_level(volc_load).plot(marker=\".\", color='red', linestyle=\"None\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "a38734c4-4ced-4889-a3fd-e0fa456d5bc5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22212 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "_ = empirical_return_level(new_contmass).plot(marker=\".\", color='red', linestyle=\"None\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "5bc03bb2-5be8-491a-93c3-fff35395f81c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 8.67,   nan,   nan, ...,   nan,   nan, 13.32], shape=(11389,))"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "new_contmass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "612c5b48-644d-4d70-a047-8364dec86cbc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted parameters:\n",
      "Shape: -1.34426, Location: 2.49871, Scale: 3.16342\n"
     ]
    }
   ],
   "source": [
    "sigl_shape, sigl_loc, sigl_scale = gev.fit(mass)\n",
    "print(f\"Fitted parameters:\\nShape: {sigl_shape:.5f}, Location: {sigl_loc:.5f}, Scale: {sigl_scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "d1a74606-6a2f-47fe-9779-a3df15bfc813",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted parameters:\n",
      "Shape: -0.46851, Location: 107.23891, Scale: 51.76101\n"
     ]
    }
   ],
   "source": [
    "lin_shape, lin_loc, lin_scale = gev.fit(mis3_lin_volc_load)\n",
    "print(f\"Fitted parameters:\\nShape: {lin_shape:.5f}, Location: {lin_loc:.5f}, Scale: {lin_scale:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "89dd8930-dd2a-41dc-9e58-319270639189",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ee804c249b0>]"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 65516 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(-1.0*ngrip_t_pedro_2018[\"Year\"],ngrip_t_pedro_2018[\"Temperature\"]) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1afcec00-2798-4e6d-86e7-f85dd3de4571",
   "metadata": {},
   "source": [
    "#### scipy.genextreme uses 1/c = -1/$\\xi$ for the shape function(c)\n",
    "$$\n",
    "P(x ; \\mu, \\sigma, \\xi) = \\frac{1}{\\sigma} t(x)^{\\xi+1} e^{-t(x)}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "t(x)= \\begin{cases}\n",
    "\\left[1+\\xi\\left(\\frac{x-\\mu}{\\sigma}\\right)\\right]^{-1 / \\xi} & \\text{if } \\xi \\neq 0 \\\\\n",
    "e^{-(x-\\mu) / \\sigma} & \\text { if } \\xi=0\\end{cases}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "138b3077-616f-4d8c-9ba1-a6510a328c91",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 367040 characters]",
      "text/plain": [
       "<Figure size 1600x1300 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "import numpy as np\n",
    "#from scipy.stats import gev\n",
    "\n",
    "# Set up the figure and gridspec\n",
    "fig = plt.figure(figsize=(16, 13))  # Adjusted size to accommodate all plots\n",
    "gs = gridspec.GridSpec(3, 2, height_ratios=[1, 1, 1])\n",
    "\n",
    "# Common settings\n",
    "SMALL_SIZE = 18\n",
    "MEDIUM_SIZE = 22\n",
    "BIGGER_SIZE = 25\n",
    "\n",
    "plt.rc('font', size=SMALL_SIZE)\n",
    "plt.rc('axes', titlesize=SMALL_SIZE)\n",
    "plt.rc('axes', labelsize=MEDIUM_SIZE)\n",
    "plt.rc('xtick', labelsize=SMALL_SIZE)\n",
    "plt.rc('ytick', labelsize=SMALL_SIZE)\n",
    "plt.rc('legend', fontsize=SMALL_SIZE)\n",
    "###################################################################################################\n",
    "# First plot (wide)\n",
    "ax1 = fig.add_subplot(gs[0, :])\n",
    "t_shift = 2000 - 1892\n",
    "t_start = -11389 + 1892\n",
    "eleg, nleg, sleg = True, True, True\n",
    "for e in range(len(timevar)):\n",
    "    if new_latitude[e] > 20.0:\n",
    "        if nleg:\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='b', label=\"Northern Hemisphere Extratropics\")\n",
    "            nleg=False\n",
    "        else:\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='b')\n",
    "    elif new_latitude[e] >= -20 and new_latitude[e] <= 20:\n",
    "        if eleg:\n",
    "            eleg=False\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='k', label=\"Equatorial\")\n",
    "        else:\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='k')\n",
    "    elif new_latitude[e] < -20.:\n",
    "        if sleg:\n",
    "            sleg=False\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='r', label=\"Southern Hemisphere Extratropics\")\n",
    "        else:\n",
    "            ax1.plot([t_start + timevar[e], t_start + timevar[e]], [0, new_contmass[e]], color='r')\n",
    "#ax1.set_ylabel('Load (Tg SO$_4$)')\n",
    "#ax1.set_xlabel('Time (Year CE)')\n",
    "ax1.legend(fontsize=14)\n",
    "ax1.set_title(\"Holocene Volcanic Eruption Data\", loc='left', size=18)\n",
    "setaxesfsize(ax1,14,r'Time ($Year~CE$)',r'Stratospheric Load ($Tg$)',16)\n",
    "ax1.set_ylim(0,700)\n",
    "########################################################################################################\n",
    "# Second plot (wide)\n",
    "ax2 = fig.add_subplot(gs[1, :])\n",
    "eleg, nleg, sleg = True, True, True\n",
    "for e in range(len(gicc05_age)):\n",
    "    if site_loc[e] == \"NHHL\":\n",
    "        if nleg:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='b', label=\"NH Extratropics\")\n",
    "            nleg=False\n",
    "        else:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='b')\n",
    "    elif site_loc[e] == \"LL or SH\":\n",
    "        if eleg:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='k', label=\"Equatorial\")\n",
    "            eleg=False\n",
    "        else:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='k')\n",
    "    elif site_loc[e] == \"SH\":\n",
    "        if sleg:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='r', label=\"SH Extratropics\")\n",
    "            sleg=False\n",
    "        else:\n",
    "            ax2.plot([-gicc05_age[e],-gicc05_age[e]],[0,volc_load[e]],color='r')\n",
    "#ax2.set_ylabel(\"Load (Tg SO$_4$)\")\n",
    "#ax2.set_xlabel('Time (Year BP)')\n",
    "ax2.set_ylim(0,700)\n",
    "ax2.set_xlim(-60000,-24000)\n",
    "ax2.legend(fontsize=14, ncol=3)\n",
    "ax2.set_title(\"Marine Isotope Stage 3 (MIS3) Volcanic Eruption Data and North Greenland Ice Core Project (NGRIP) Temperature\", loc='left', size=18)\n",
    "setaxesfsize(ax2,14,r'Time ($Year~BP$)',r'Stratospheric Load ($Tg$)', 16)\n",
    "\n",
    "# add in NGRIP temperature data\n",
    "ax2r = ax2.twinx()\n",
    "ax2r.plot(-1.0*ngrip_t_pedro_2018[\"Year\"],ngrip_t_pedro_2018[\"Temperature\"], color='green', linewidth=1.0, label='NGRIP Temperature')\n",
    "setaxesfsize(ax2r,14,'',r'Temperature ($\\degree C$)', 16)\n",
    "ax2r.set_ylim(-60,-25)\n",
    "ax2r.legend(fontsize=14, ncol=1, loc=(0.802,0.67))\n",
    "\n",
    "## In this subplot we will add horizontal lines at 230, 116, 60 Tg SO4 to show the simulated levels\n",
    "simlw = 0.5\n",
    "ax2.axhline(y=230, color='#222222', linestyle='--', label='230 Tg', linewidth=simlw*3)\n",
    "ax2.axhline(y=116, color='#444444', linestyle='--', label='115 Tg', linewidth=simlw*1.5)\n",
    "ax2.axhline(y=58, color='#666666', linestyle='--', label='58 Tg', linewidth=simlw)\n",
    "ax2.legend(fontsize=14, ncol=3)\n",
    "\n",
    "###################################################################################################\n",
    "# Third plot (square)\n",
    "ax3 = fig.add_subplot(gs[2, 0])\n",
    "hist_sigl, bins_sigl = np.histogram(mass, bins=31, range=(0, 600), density=True)\n",
    "ax3.bar(bins_sigl[:-1], hist_sigl, width=10, align='edge', alpha=0.3, color='red', label=\"Holocene data\")\n",
    "hist_lin, bins_lin = np.histogram(mis3_lin_volc_load, bins=31, range=(0, 600), density=True)\n",
    "ax3.bar(bins_lin[:-1], hist_lin, width=10, align='edge', alpha=0.3, color='blue', label=\"MIS3 data\")\n",
    "ax3.set_ylabel(\"Normalized Occurrences\")\n",
    "ax3.set_xlabel(r\"Stratospheric Load ($Tg$)\", size=16)\n",
    "setaxesfsize(ax3, 14, r\"Stratospheric Load ($Tg$)\", \"Normalized Occurrences\", 16)\n",
    "ax3.set_ylim(0, 0.01)\n",
    "ax3.legend(fontsize=14, ncol=2)\n",
    "\n",
    "sigl_l = sigl_loc + sigl_scale / sigl_shape\n",
    "sigl_xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "sigl_yy = gev.pdf(sigl_xx, sigl_shape, sigl_loc, sigl_scale)\n",
    "ax3.plot(sigl_xx, sigl_yy, 'red', linewidth = 1.5, label = \"Sigl GEV fit\")\n",
    "\n",
    "lin_l = lin_loc + lin_scale / lin_shape\n",
    "lin_xx = np.linspace(l+0.00001, l+0.00001+600, num=71)\n",
    "lin_yy = gev.pdf(lin_xx, lin_shape, lin_loc, lin_scale)\n",
    "ax3.plot(lin_xx, lin_yy, 'blue', linewidth = 1.5, label = \"Lin GEV fit\")\n",
    "\n",
    "ax3.set_title(\"Generalized Extreme Value Distributions\", loc='left', size=18)\n",
    "\n",
    "# Add text box with parameters\n",
    "sigl_params = sigl_shape, sigl_loc, sigl_scale\n",
    "lin_params = lin_shape, lin_loc, lin_scale\n",
    "bbox = dict(boxstyle=\"round\", color=\"white\", linewidth=2)\n",
    "## scipy uses 1/c instead of -1/xi , so we have to multiply by -1\n",
    "ax3.text(0.5, 0.61, r'$\\xi$=%.2f' '\\n' r'$\\mu$=%.2f' '\\n' r'$\\sigma$=%.2f' % (-1.0*sigl_params[0], sigl_params[1], sigl_params[2]),\n",
    "         transform=ax3.transAxes, bbox=bbox, color='red', size=14)\n",
    "ax3.text(0.82, 0.61, r'$\\xi$=%.2f' '\\n' r'$\\mu$=%.2f' '\\n' r'$\\sigma$=%.2f' % (-1.0*lin_params[0], lin_params[1], lin_params[2]),\n",
    "         transform=ax3.transAxes, bbox=bbox, color='blue', size=14)\n",
    "###############################################################################################################################\n",
    "\n",
    "# Fourth plot (square)\n",
    "ax4 = fig.add_subplot(gs[2, 1])\n",
    "#periods = np.arange(2, 1000, 1)\n",
    "#quantiles = 1 - 1 / periods\n",
    "#sigl_levels = estimate_return_level(quantiles, sigl_loc, sigl_scale, sigl_shape)\n",
    "#lin_levels = estimate_return_level(quantiles, lin_loc, lin_scale, lin_shape)\n",
    "ax4.plot(sigl_periods*sigl_tscale, sigl_levels, \"-\", color='red', label='Holocene Data')\n",
    "ax4.plot(lin_periods*lin_tscale, lin_levels, \"-\", color='blue', label='MIS3 Data')\n",
    "ax4.set_xlabel(\"Return Period (years)\")\n",
    "ax4.set_ylabel(\"Return Level (Tg)\")\n",
    "ax4.set_ylim(0,500)\n",
    "ax4.set_xlim(0,10000)\n",
    "ax4.grid(True)\n",
    "setaxesfsize(ax4,14,r\"Return Period ($Years$)\",r\"Return Level ($Tg$)\",16)\n",
    "ax4.set_title(\"Holocene and MIS3 Return Times\", loc='left', size=18)\n",
    "ax4.legend(fontsize=14, ncol=2)\n",
    "# Adjust the layout\n",
    "plt.tight_layout()\n",
    "\n",
    "## Add a horizontal line going from 230, 116, 60 Tg SO4 on the y-axis to the y value on the red curve \n",
    "## Add a vertical line going from 0 on the x-axis to y value on the red curve at 230, 116, 60 Tg SO4 on the y-axis\n",
    "## I need it to stop when the lines intercept the red curve\n",
    "\n",
    "# Find the x-values where the horizontal lines intersect the red curve\n",
    "x_230 = np.interp(230, sigl_levels, sigl_periods*sigl_tscale)\n",
    "x_116 = np.interp(116, sigl_levels, sigl_periods*sigl_tscale)\n",
    "x_58 = np.interp(58, sigl_levels, sigl_periods*sigl_tscale)\n",
    "\n",
    "# Horizontal lines\n",
    "ax4.plot([0, x_230], [230, 230], color='#222222', linestyle='--', linewidth=simlw*3)\n",
    "ax4.plot([0, x_116], [116, 116], color='#444444', linestyle='--', linewidth=simlw*1.5)\n",
    "ax4.plot([0, x_58], [58, 58], color='#666666', linestyle='--', linewidth=simlw)\n",
    "\n",
    "# Vertical lines\n",
    "ax4.plot([x_230, x_230], [0, 230], color='#222222', linestyle='--', linewidth=simlw*3)\n",
    "ax4.plot([x_116, x_116], [0, 116], color='#444444', linestyle='--', linewidth=simlw*1.5)\n",
    "ax4.plot([x_58, x_58], [0, 58], color='#666666', linestyle='--', linewidth=simlw)\n",
    "\n",
    "# Add intersection points\n",
    "ax4.plot(x_230, 230, 'ro')\n",
    "ax4.plot(x_116, 116, 'ro')\n",
    "ax4.plot(x_58, 58, 'ro')\n",
    "\n",
    "\n",
    "# add a,b,c,d\n",
    "lfsize = 24\n",
    "ax1.text(-0.05, 1.05, 'A',\n",
    "         transform=ax1.transAxes, size=lfsize, weight='bold')\n",
    "ax2.text(-0.05, 1.05, 'B',\n",
    "         transform=ax2.transAxes, size=lfsize, weight='bold')\n",
    "ax3.text(-0.13, 1.05, 'C',\n",
    "         transform=ax3.transAxes, size=lfsize, weight='bold')\n",
    "ax4.text(-0.13, 1.05, 'D',\n",
    "         transform=ax4.transAxes, size=lfsize, weight='bold')\n",
    "\n",
    "# Save the plot\n",
    "\n",
    "icepaperdir = \"/home/guido/Dropbox/plots/volcanic/\"\n",
    "saveplot = False\n",
    "if saveplot:\n",
    "    plt.savefig(icepaperdir+\"Figure1.png\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n",
    "    plt.savefig(icepaperdir+\"Figure1.pdf\", dpi=300, bbox_inches='tight', facecolor=fig.get_facecolor(), transparent=False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff18c5a9-695b-4f0f-8291-f768a408858d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (geosci)",
   "language": "python",
   "name": "python3"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
