global_sum_differences.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4698e635-efdc-47d7-82a0-b5b8095f3bb4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import sys\n",
    "# !{sys.executable} -m pip install mpl_toolkits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38858852",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import glacierml as gl\n",
    "import configparser\n",
    "import matplotlib.pyplot as plt\n",
    "# from mpl_toolkits.basemap import Basemap\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "pd.set_option('display.max_column',None)\n",
    "# import sys\n",
    "cols = []\n",
    "for i in range(273):\n",
    "    cols.append(i)\n",
    "# !{sys.executable} -m pip install basemap\n",
    "import path_manager as pm\n",
    "[\n",
    "        home_path, data_path, RGI_path, glathida_path, \n",
    "        ref_path, coregistration_testing_path, \n",
    "        arch_test_path, LOO_path\n",
    "] = pm.set_paths()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86a78764",
   "metadata": {},
   "source": [
    "#### Let's look at how we compare to Farinotti et al. 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fa9230d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.read_pickle('final.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "086e7010",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shelves = df[df['TermType'] == 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c612dc1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(shelves)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abb7916d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df = pd.read_pickle('/data/fast1/glacierml/data/final_df.pkl')\n",
    "E_A = df['Area']\n",
    "E_H = np.mean(df[cols],axis = 1)\n",
    "si = df['sig_k_ind'] # = (Var(H)*Var(A) + A^2Var(H) + H^2Var(A))\n",
    "\n",
    "sp = np.sum(si)\n",
    "\n",
    "se = np.sqrt(sp)\n",
    "\n",
    "Z = (1.96)\n",
    "Vlb = np.round( (np.sum(E_H * E_A) - (Z*se)) , 0) / 1e3\n",
    "Vub = np.round( (np.sum(E_H * E_A) + (Z*se)) , 0) / 1e3\n",
    "\n",
    "print(f'[{Vlb},{Vub}]  * 10^3 km^3')\n",
    "\n",
    "print(f'Mid CI: {np.round((Vub + Vlb) / 2,3)} * 10^3 km^3')\n",
    "\n",
    "print(f'CI Half Width: {np.round((Vub - Vlb) / 2,3)} * 10^3 km^3')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c56605c0",
   "metadata": {},
   "source": [
    "#### Okay, that is our global sum. Now lets look at our global sum when just comparing to Millan. To accomplish this, we exclude shelf supported glaciers in Antarctica "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea4ff1dd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# df = pd.read_pickle('/data/fast1/glacierml/data/final_df.pkl')\n",
    "\n",
    "df_no_shelf = df.drop(df[\n",
    "    (df['TermType'] == 5) & (df['region'] == '19')\n",
    "].index)\n",
    "\n",
    "# df = pd.read_pickle('/data/fast1/glacierml/data/final_df.pkl')\n",
    "E_A = df_no_shelf['Area']\n",
    "E_H = np.mean(df_no_shelf[cols],axis = 1)\n",
    "si = df_no_shelf['sig_k_ind'] # = (Var(H)*Var(A) + A^2Var(H) + H^2Var(A))\n",
    "\n",
    "sp = np.sum(si)\n",
    "\n",
    "se = np.sqrt(sp)\n",
    "\n",
    "Z = (1.96)\n",
    "Vlb = np.round( (np.sum(E_H * E_A) - (Z*se)) , 0) / 1e3\n",
    "Vub = np.round( (np.sum(E_H * E_A) + (Z*se)) , 0) / 1e3\n",
    "\n",
    "print(f'[{Vlb},{Vub}]  * 10^3 km^3')\n",
    "\n",
    "print(f'Mid CI: {np.round((Vub + Vlb) / 2,3)} * 10^3 km^3')\n",
    "\n",
    "print(f'CI Half Width: {np.round((Vub - Vlb) / 2,3)} * 10^3 km^3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "393cb578",
   "metadata": {},
   "outputs": [],
   "source": [
    "140.6 - 111.399"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e88e625c",
   "metadata": {},
   "source": [
    "#### Interesting, now how do we compare to Farinotti?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5571d3df",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_pickle('final.pkl')\n",
    "cols,dft = gl.load_LOO_data(home_path,include_refs = True)\n",
    "df = pd.merge(df,dft,how = 'inner', on = list(dft)[:-1])\n",
    "E_A = df['Area']\n",
    "E_H = np.mean(df[cols],axis = 1)\n",
    "si = df['sig_k_ind'] # = (Var(H)*Var(A) + A^2Var(H) + H^2Var(A))\n",
    "\n",
    "sp = np.sum(si)\n",
    "\n",
    "se = np.sqrt(sp)\n",
    "\n",
    "Z = (1.96)\n",
    "Vlb = np.round( (np.sum(E_H * E_A) - (Z*se)) , 0) / 1e3\n",
    "Vub = np.round( (np.sum(E_H * E_A) + (Z*se)) , 0) / 1e3\n",
    "\n",
    "print(f'[{Vlb},{Vub}]  * 10^3 km^3')\n",
    "\n",
    "print(f'Mid CI: {np.round((Vub + Vlb) / 2,3)} * 10^3 km^3')\n",
    "\n",
    "print(f'CI Half Width: {np.round((Vub - Vlb) / 2,3)} * 10^3 km^3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22eaee03-7088-4924-bfaf-bf473df4ad88",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "158.48 - 121.81"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4fea5af",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### We are looking for $\\approx 36.7 \\times 10^3 \\text{km}^3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc3afefe-8496-4da6-b55e-b8f97870c0be",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df = pd.concat(\n",
    "    [\n",
    "        df,\n",
    "        pd.Series( (np.mean(df[cols],axis = 1) * df['Area'] )- (df['FMT'] * df['Area']), name = 'r' )\n",
    "    ], axis = 1\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3d62fa3",
   "metadata": {},
   "outputs": [],
   "source": [
    "running_sum = np.cumsum(\n",
    "    df['r'].sort_values(ascending = False).reset_index().drop('index',axis = 1)\n",
    ")\n",
    "running_sum.index = range(1, len(running_sum) + 1)\n",
    "\n",
    "plt.plot(running_sum)\n",
    "plt.minorticks_on()\n",
    "plt.grid(which = 'both')\n",
    "plt.xlabel('Glacier Index')\n",
    "plt.ylabel('Difference of Volume')\n",
    "plt.title('Cumulative Sum Difference of Volumes')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a60e7e1b",
   "metadata": {},
   "source": [
    "#### Okay, so we have a handful of large positive and large negative discrepancies. \n",
    "#### What do they look like close up?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5b8866d",
   "metadata": {},
   "outputs": [],
   "source": [
    "running_sum1 = np.cumsum(\n",
    "    df['r'].sort_values(ascending = False).reset_index().drop('index',axis = 1)\n",
    ")\n",
    "running_sum1.index = range(1, len(running_sum1) + 1)\n",
    "\n",
    "plt.plot(running_sum1[0:100])\n",
    "plt.minorticks_on()\n",
    "plt.grid(which = 'both')\n",
    "# plt.xscale('log')\n",
    "plt.xlabel('Glacier Index')\n",
    "plt.ylabel('Difference of Volume (km$^3$)')\n",
    "plt.title('Cumulative Sum Difference of Volumes')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f311ead9",
   "metadata": {},
   "outputs": [],
   "source": [
    "running_sum2 = np.cumsum(\n",
    "    df['r'].sort_values(ascending = True).reset_index().drop('index',axis = 1)\n",
    ")\n",
    "running_sum2.index = range(1, len(running_sum2) + 1)\n",
    "# plt.plot(running_sum1[0:500],label = 'Positive Differences')\n",
    "\n",
    "plt.plot(running_sum2[0:1000],label = 'Negative Differences')\n",
    "plt.minorticks_on()\n",
    "plt.grid(which = 'both')\n",
    "# plt.xscale('log')\n",
    "# plt.legend()\n",
    "plt.xlabel('Glacier Index')\n",
    "plt.ylabel('Difference of Volume (km$^3$)')\n",
    "plt.title('Cumulative Sum Difference of Volumes')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd9058b",
   "metadata": {},
   "outputs": [],
   "source": [
    "running_sum2 = np.cumsum(\n",
    "    df['r'].sort_values(ascending = True).reset_index().drop('index',axis = 1)\n",
    ")\n",
    "running_sum2.index = range(1, len(running_sum2) + 1)\n",
    "plt.plot(running_sum1[0:1000],label = 'Positive Differences')\n",
    "\n",
    "plt.plot(running_sum2[0:1000],label = 'Negative Differences')\n",
    "plt.minorticks_on()\n",
    "# plt.grid(which = 'both')\n",
    "# plt.xscale('log')\n",
    "plt.legend()\n",
    "plt.xlabel('Glacier Index')\n",
    "plt.ylabel('Difference of Volume (km$^3$)')\n",
    "plt.title('Cumulative Sum Difference of Volumes')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53295e88",
   "metadata": {},
   "source": [
    "#### What percentage of the global glacier volume discrepancy can be accounted for with these glaciers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2223b9e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = df.sort_values('r', ascending=True)\n",
    "df2 = df.sort_values('r', ascending=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4288e7ac",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "counts = []\n",
    "percs = []\n",
    "n = 0\n",
    "for i in tqdm(range(10,540,10)):\n",
    "#     print(i)\n",
    "    n =+ i\n",
    "#     print(n*2)\n",
    "\n",
    "    percent_of_df = np.round((n*2) / len(df) * 100,4)\n",
    "#     print(f'Percent of global glaciers considered = {percent_of_df}')\n",
    "\n",
    "    neg = (np.sum(df1['r'].head(n)) / 1e3)\n",
    "#     print(neg)\n",
    "\n",
    "    pos = (np.sum(df2['r'].tail(n)) / 1e3)\n",
    "    d = abs(pos + neg) / (-t) * 100\n",
    "    percs.append(d)\n",
    "    counts.append(percent_of_df)\n",
    "#     print(f'Volume Difference contained = {d}')\n",
    "#     print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04632e04",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(counts,percs)\n",
    "plt.xlabel('Percent of RGI')\n",
    "plt.ylabel('Percent of Global Volume Discrepancy')\n",
    "plt.minorticks_on()\n",
    "plt.grid(which = 'both')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20ef97f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI = gl.load_RGI()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 20))\n",
    "\n",
    "ax1 = axes[0]\n",
    "m1 = Basemap(projection='cyl', ax=ax1, resolution='c')\n",
    "m1.scatter(RGI['CenLon'], RGI['CenLat'], c='blue', alpha=0.5, label='RGI Outline')\n",
    "\n",
    "m1.scatter(df1['CenLon'].head(500), df1['CenLat'].head(500), c='red', alpha=1, label='Largest Negative Differences')\n",
    "m1.drawcoastlines()\n",
    "ax1.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "ax2 = axes[1]\n",
    "m2 = Basemap(projection='cyl', ax=ax2, resolution='c')\n",
    "m2.scatter(RGI['CenLon'], RGI['CenLat'], c='blue', alpha=0.5, label='RGI Outline')\n",
    "\n",
    "m2.scatter(df2['CenLon'].tail(500), df2['CenLat'].tail(500), c='orange', alpha=1, label='Largest Positive Differences')\n",
    "m2.drawcoastlines()\n",
    "ax2.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3a6bf73",
   "metadata": {},
   "outputs": [],
   "source": [
    "neg_100 = df1.head(100)\n",
    "pos_100 = df1.tail(100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6aa584e",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = gl.coregister_data('4')\n",
    "\n",
    "g = g.drop(g[g['RGIId'].duplicated(keep = False)].index)\n",
    "g = g.sample(frac = 1,random_state = 0)\n",
    "g = g.reset_index().drop('index', axis = 1)\n",
    "g['region'] = g['RGIId'].str[6:8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69541b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in neg_100['region'].sort_values().unique():\n",
    "    negs = neg_100[neg_100['region'] == i]\n",
    "    category_to_count = i\n",
    "    n = (neg_100['region'] == category_to_count).sum()\n",
    "    \n",
    "    category_to_count2 = 5\n",
    "    m = (negs['TermType'] == category_to_count2).sum()\n",
    "    \n",
    "    category_to_count2 = 1\n",
    "    l = (negs['Form'] == category_to_count2).sum()\n",
    "    \n",
    "    search_list = np.ndarray.tolist(negs['RGIId'].values)\n",
    "    pattern = '|'.join(search_list)\n",
    "    gs = g[g['RGIId'].str.contains(pattern)]\n",
    "    o = g['RGIId'].str.contains(pattern).sum()   \n",
    "    print(\n",
    "        f'RGI Region {i} has {n} / 500 most negative discrepant glaciers,\\n'\n",
    "        f'{m} of which are shelf supported,\\n'\n",
    "        f'{l} of which are ice caps,\\n'\n",
    "        f'{o} of which are in training data.'\n",
    "    )\n",
    "    if o != 0:\n",
    "        \n",
    "    \n",
    "        category_to_count3 = 5\n",
    "        category_to_count4 = 1\n",
    "        p = (gs['TermType'] == category_to_count3).sum()\n",
    "        q = (gs['Form'] == category_to_count4).sum()   \n",
    "        print(\n",
    "            f'Of the glaciers contained in training data, '\n",
    "            f'{p} are shelf supported and {q} are ice caps'\n",
    "        )\n",
    "        problem_train_glacs = gs['Name'].values\n",
    "        problem_train_glacs_id = gs['RGIId'].values\n",
    "        print(f'These glaciers are {problem_train_glacs}')\n",
    "        print(f'RGIId: {problem_train_glacs_id}')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b367cff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "321a64d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(g['Area'],g['Area']*(g['Thickness']/1e3))\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9179a697",
   "metadata": {},
   "outputs": [],
   "source": [
    "sum(df['r'][df['TermType'] == 5]) / 1e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7da92d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(g[g['TermType'] == 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8add3cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfr0 = (df[(df['TermType'] == 1) & (df['Form'] == 1)])\n",
    "dfr1 = (df[(df['TermType'] == 5) & (df['Form'] == 1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff4d06a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfr0 = (df[\n",
    "    (df['TermType'] == 1) & \n",
    "#     (df['Form'] == 0) & \n",
    "    (df['region'] == '07')\n",
    "])\n",
    "dfr0 = (df[\n",
    "    (df['TermType'] == 5) & \n",
    "#     (df['Form'] == 0) & \n",
    "    (df['region'] == '07')\n",
    "])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eaecbc41",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.scatter(dfr0['Area'],-dfr0['r'])\n",
    "plt.scatter(dfr0['Area'],dfr0['r'])\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a7d889",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(dfr1['Area'],-dfr1['r'])\n",
    "\n",
    "plt.scatter(dfr1['Area'],dfr1['r'])\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33de39f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfr1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0638b12",
   "metadata": {},
   "outputs": [],
   "source": [
    "(df['r'][(df['TermType'] == 1) & (df['region'] == '05')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd6a3d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "(df['r'][(df['TermType'] == 0) & (df['region'] == '05')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16a80ea2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3953a91b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36384614",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in pos_100['region'].sort_values().unique():\n",
    "    poss = pos_100[pos_100['region'] == i]\n",
    "    category_to_count = i\n",
    "    n = (pos_100['region'] == category_to_count).sum()\n",
    "    \n",
    "    category_to_count2 = 5\n",
    "    m = (poss['TermType'] == category_to_count2).sum()\n",
    "    \n",
    "    category_to_count2 = 1\n",
    "    l = (poss['Form'] == category_to_count2).sum()\n",
    "    \n",
    "    search_list = np.ndarray.tolist(poss['RGIId'].values)\n",
    "    pattern = '|'.join(search_list)\n",
    "    gs = g[g['RGIId'].str.contains(pattern)]\n",
    "    o = g['RGIId'].str.contains(pattern).sum()   \n",
    "    print(\n",
    "        f'RGI Region {i} has {n} / 500 most positive discrepant glaciers,\\n'\n",
    "        f'{m} of which are shelf supported,\\n'\n",
    "        f'{l} of which are ice caps,\\n'\n",
    "        f'{o} of which are in training data.'\n",
    "    )\n",
    "    if o != 0:\n",
    "        \n",
    "    \n",
    "        category_to_count3 = 5\n",
    "        category_to_count4 = 1\n",
    "        p = (gs['TermType'] == category_to_count3).sum()\n",
    "        q = (gs['Form'] == category_to_count4).sum()   \n",
    "        print(\n",
    "            f'Of the glaciers contained in training data, '\n",
    "            f'{p} are shelf supported and {q} are ice caps'\n",
    "        )\n",
    "        problem_train_glacs = gs['Name'].values\n",
    "        problem_train_glacs_id = gs['RGIId'].values\n",
    "        print(f'These glaciers are {problem_train_glacs}')\n",
    "        print(f'RGIId: {problem_train_glacs_id}')\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bea78cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "barnes = df[df['RGIId'] == 'RGI60-04.06187']\n",
    "windy = df[df['RGIId'] == 'RGI60-09.00807']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b8181fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "windy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe30a57c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head(100)[df1['region'] == '03'][['CenLat','CenLon','Name']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77641c7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head(100)[df1['region'] == '04'][['CenLat','CenLon','Name']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f78050",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head(100)[df1['region'] == '05'][['CenLat','CenLon','Name']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09855608",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head(100)[df1['region'] == '06'][['CenLat','CenLon','Name']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6095cfc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1.head(100)[df1['region'] == '07'][['CenLat','CenLon','Name']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de61d3ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "(df1.head(100)[df1['region'] == '09'][['CenLat','CenLon','Name']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a550ac8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "antarctic_discrepant_index = df1.head(100)[df1['region'] == '19'][['CenLat','CenLon','Name']].index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a78756a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_pickle('/data/fast1/glacierml/data/final_df.pkl')\n",
    "df = df.drop(antarctic_discrepant_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba2ce4fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = pd.read_pickle('/data/fast1/glacierml/data/final_df.pkl')\n",
    "E_A = df['Area']\n",
    "E_H = np.mean(df[cols],axis = 1)\n",
    "si = df['sig_k_ind'] # = (Var(H)*Var(A) + A^2Var(H) + H^2Var(A))\n",
    "\n",
    "sp = np.sum(si)\n",
    "\n",
    "se = np.sqrt(sp)\n",
    "\n",
    "Z = (1.96)\n",
    "Vlb = np.round( (np.sum(E_H * E_A) - (Z*se)) , 0) / 1e3\n",
    "Vub = np.round( (np.sum(E_H * E_A) + (Z*se)) , 0) / 1e3\n",
    "\n",
    "print(f'[{Vlb},{Vub}]  * 10^3 km^3')\n",
    "\n",
    "print(f'Mid CI: {np.round((Vub + Vlb) / 2,3)} * 10^3 km^3')\n",
    "\n",
    "print(f'CI Half Width: {np.round((Vub - Vlb) / 2,3)} * 10^3 km^3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "028b79bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "gf = g[g['Form'] == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d798ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16424e68",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(g['Area'],g['Thickness'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1c608cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI = gl.load_RGI()\n",
    "\n",
    "fig,ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 20))\n",
    "\n",
    "# ax1 = axes[0]\n",
    "m1 = Basemap(projection='cyl', ax=ax, resolution='c')\n",
    "m1.scatter(RGI['CenLon'], RGI['CenLat'], c='blue', alpha=0.5, label='RGI Outline')\n",
    "m1.scatter(g['CenLon'], g['CenLat'], c='red', alpha=0.5, label='RGI Outline')\n",
    "\n",
    "m1.scatter(gf['CenLon'], gf['CenLat'], c='orange', alpha=1, label='Largest Negative Differences')\n",
    "m1.drawcoastlines()\n",
    "ax1.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "# ax2 = axes[1]\n",
    "# m2 = Basemap(projection='cyl', ax=ax2, resolution='c')\n",
    "# m2.scatter(RGI['CenLon'], RGI['CenLat'], c='blue', alpha=0.5, label='RGI Outline')\n",
    "\n",
    "# m2.scatter(df2['CenLon'].tail(500), df2['CenLat'].tail(500), c='orange', alpha=1, label='Largest Positive Differences')\n",
    "# m2.drawcoastlines()\n",
    "# ax2.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ec74263",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8672277",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e5c1ad6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eef3b933",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a37c6232",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74a89d9b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f94be78e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e92838a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f11fe71",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "24941219",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Status'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bb6e504",
   "metadata": {},
   "outputs": [],
   "source": [
    "g['Form'].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c4c7d20",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d04ccf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df[df['Form'] == 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2ea2868",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df[df['Form'] == 0]) / len(df) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d88e3069",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df[df['Form'] == 1]) / len(df) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1262e916",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(g[g['Form'] == 0]) / len(g) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fbc6c31",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(g[g['Form'] == 1]) / len(g) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc83f117",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(g[(g['Form'] == 1) & (g['region'] == '19')]) / len(g) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "795ed2a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "gf.sort_values('RGIId')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5543331",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00cd4569",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI['region'][RGI['TermType'] == 5].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3f1e063",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[g['TermType'] == 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a6970cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[(g['Form'] == 1) & (g['region'] == '04')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6634dc5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[(g['Form'] == 1) & (g['region'] == '07')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69295de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[(g['Form'] == 1) & (g['region'] == '09')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc47d8ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[(g['Form'] == 1) & (g['region'] == '17')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50c6a1b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[(g['Form'] == 1) & (g['region'] == '19')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6621fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[g['Form'] == 1]['Area'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62347a80",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI[RGI['Form'] == 1]['Area'].sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e28feb5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "g[g['TermType'] == 5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c61397ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "(RGI[RGI['TermType'] == 5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81c002a",
   "metadata": {},
   "outputs": [],
   "source": [
    "RGI = gl.load_RGI(RGI_path)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 20))\n",
    "\n",
    "ax1 = axes[0]\n",
    "m1 = Basemap(projection='cyl', ax=ax1, resolution='c')\n",
    "m1.scatter(\n",
    "    (RGI[RGI['TermType'] == 5])['CenLon'], (RGI[RGI['TermType'] == 5])['CenLat'], \n",
    "    c='blue', alpha=1,\n",
    "    label='RGI Shelf Terminating Glaciers'\n",
    ")\n",
    "\n",
    "# m1.scatter(df1['CenLon'].head(500), df1['CenLat'].head(500), c='red', alpha=1, label='Largest Negative Differences')\n",
    "m1.drawcoastlines()\n",
    "ax1.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "ax2 = axes[1]\n",
    "m2 = Basemap(projection='cyl', ax=ax2, resolution='c')\n",
    "m2.scatter(\n",
    "    df1['CenLon'].head(100), df1['CenLat'].head(100), \n",
    "    c='red', alpha=1, \n",
    "    label='100 Largest Negative Differences'\n",
    ")\n",
    "\n",
    "# m2.scatter(df2['CenLon'].tail(500), df2['CenLat'].tail(500), c='orange', alpha=1, label='Largest Positive Differences')\n",
    "m2.drawcoastlines()\n",
    "ax2.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41ed0dc1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "RGI = gl.load_RGI()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 20))\n",
    "\n",
    "ax1 = axes[0]\n",
    "m1 = Basemap(projection='cyl', ax=ax1, resolution='c')\n",
    "m1.scatter(\n",
    "    (RGI[RGI['Form'] == 1])['CenLon'], (RGI[RGI['Form'] == 1])['CenLat'], \n",
    "    c='blue', alpha=1,\n",
    "    label='RGI Ice Caps'\n",
    ")\n",
    "\n",
    "# m1.scatter(df1['CenLon'].head(500), df1['CenLat'].head(500), c='red', alpha=1, label='Largest Negative Differences')\n",
    "m1.drawcoastlines()\n",
    "ax1.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "ax2 = axes[1]\n",
    "m2 = Basemap(projection='cyl', ax=ax2, resolution='c')\n",
    "m2.scatter(\n",
    "    df1['CenLon'].head(500), df1['CenLat'].head(500), \n",
    "    c='red', alpha=1, \n",
    "    label='100 Largest Negative Differences'\n",
    ")\n",
    "\n",
    "# m2.scatter(df2['CenLon'].tail(500), df2['CenLat'].tail(500), c='orange', alpha=1, label='Largest Positive Differences')\n",
    "m2.drawcoastlines()\n",
    "ax2.legend(bbox_to_anchor=(0.25, 0.56), facecolor='w', framealpha=1, fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "glacierml (Python3.8.10)",
   "language": "python",
   "name": "glacierml"
  },
  "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.8.10"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}