Notebook_4_ghg_emission_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce6051bb",
   "metadata": {},
   "source": [
    "# Visualisation of GHG emissions from RE-Emission\n",
    "\n",
    "- Author: Tomasz Janus\n",
    "- E-mail: tomasz.k.janus@gmail.com , tomasz.janus@manchester.ac.uk\n",
    "- 22/10/2023"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6bb10f0",
   "metadata": {},
   "source": [
    "## 1. Import the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fce804f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Protocol, Dict, Protocol, List, Any\n",
    "from dataclasses import dataclass\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns # For plotting data\n",
    "\n",
    "plt.rcParams.update({\n",
    "    \"text.usetex\": True,\n",
    "})\n",
    "\n",
    "def hexbin(x, y, color, **kwargs):\n",
    "    cmap = sns.light_palette(color, as_cmap=True)\n",
    "    plt.hexbin(x, y, gridsize=15, cmap=cmap, **kwargs)\n",
    "    \n",
    "full_data_path = os.path.join('outputs','reemission','combined','combined_outputs.csv')\n",
    "elev_data_path = os.path.join(\"config\",\"elev.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "268a7a32",
   "metadata": {},
   "source": [
    "## 2. Input data cleanup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfc90980",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "data_full = pd.read_csv(full_data_path)\n",
    "elev_data = pd.read_csv(elev_data_path)\n",
    "merged_df = pd.merge(\n",
    "    data_full, elev_data[['ifc_name', 'fsl_masl']], left_on='Name', right_on='ifc_name', how='inner')\n",
    "merged_df.drop(columns=['ifc_name'], inplace=True)\n",
    "# Perform prefiltering of data outside the pipeline\n",
    "merged_df[['Soil', 'Landuse intensity', 'Treatment']] = \\\n",
    "            merged_df['Scenario'].str.split('_', expand=True)\n",
    "#Create a version with CH4 emissions without degassing\n",
    "merged_df['ch4_net_nodegassing'] = merged_df['ch4_net'] - merged_df['ch4_degassing']\n",
    "# Create a separate dataframe for plotting only\n",
    "df_plot = merged_df.copy()\n",
    "df_plot['Name_Code'] = pd.factorize(df_plot['Name'])[0]\n",
    "df_plot['Scenario_Code'] = pd.factorize(df_plot['Scenario'])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f60d198",
   "metadata": {},
   "source": [
    "## 3. Visual analysis of GHG estimation results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f98812b3",
   "metadata": {},
   "source": [
    "### 3a. Plot distributions of emissions per area for different combinations of categorical variables\n",
    "The plots show the sensitivity of the emissions model to selection of categorical inputs such as **soil type**, **landuse intensity** and **treatment factor**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ab26fb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_label_co2 = r'Net CO$_2$, gCO$_{2,eq}$ m$^{-2}$ yr$^{-1}$'\n",
    "y_label_ch4 = r'Net CH$_4$, gCO$_{2,eq}$ m$^{-2}$ yr$^{-1}$'\n",
    "fig, axes = plt.subplots(ncols=2, nrows=2, figsize=(11, 7))\n",
    "fig.suptitle(\n",
    "    r'Net CO$_2$ and CH$_4$ emissions for different soil types, landuse intensities and treatment factors', \n",
    "    fontsize=14)\n",
    "ax1, ax2, ax3, ax4 = axes.ravel()\n",
    "ax1.set_title('Mineral Soil')\n",
    "#sns.set_style('white')\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "sns.boxenplot(\n",
    "    x=\"Landuse intensity\", y=\"co2_net\",\n",
    "    hue=\"Treatment\", palette=[\"m\", \"g\", \"r\", \"y\"],\n",
    "    data=merged_df[df_plot[\"Soil\"] == \"MIN\"],\n",
    "    ax=ax1)\n",
    "ax1.legend([],[], frameon=False)\n",
    "ax1.set(xlabel=r'Landuse intensity', ylabel=y_label_co2)\n",
    "ax2.set_title('Organic Soil')\n",
    "sns.boxenplot(\n",
    "    x=\"Landuse intensity\", y=\"co2_net\",\n",
    "    hue=\"Treatment\", palette=[\"m\", \"g\", \"r\", \"y\"],\n",
    "    data=merged_df[df_plot[\"Soil\"] == \"ORG\"],\n",
    "    ax=ax2)\n",
    "ax2.legend(bbox_to_anchor=(1.4, 0.75), borderaxespad=0, title=\"Treatment\", frameon=False)\n",
    "ax2.legend([],[], frameon=False)\n",
    "ax2.set(xlabel=r'Landuse intensity', ylabel=y_label_co2)\n",
    "ax3.set_title('Mineral Soil')\n",
    "sns.boxenplot(\n",
    "    x=\"Landuse intensity\", y=\"ch4_net\",\n",
    "    hue=\"Treatment\", palette=[\"m\", \"g\", \"r\", \"y\"],\n",
    "    data=merged_df[df_plot[\"Soil\"] == \"MIN\"],\n",
    "    ax=ax3)\n",
    "ax3.legend([],[], frameon=False)\n",
    "ax3.set(xlabel=r'Landuse intensity', ylabel=y_label_ch4)\n",
    "ax4.set_title('Organic Soil')\n",
    "sns.boxenplot(\n",
    "    x=\"Landuse intensity\", y=\"ch4_net\",\n",
    "    hue=\"Treatment\", palette=[\"m\", \"g\", \"r\", \"y\"],\n",
    "    data=merged_df[df_plot[\"Soil\"] == \"ORG\"],\n",
    "    ax=ax4)\n",
    "ax4.legend(bbox_to_anchor=(1.4, 0.75), borderaxespad=0, title=\"Treatment\", frameon=False)\n",
    "ax4.set(xlabel=r'Landuse intensity', ylabel=y_label_ch4)\n",
    "sns.despine(offset=15, trim=False)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(\"figures\",\"ghg_visualisation\",\"emission_boxplots.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f11a34f",
   "metadata": {},
   "source": [
    "#### Quantify average percentage reduction of CO2 emissions in scenario Low Intensity, Mineral Soil due to implementation of further wastewater treatment technologies in the catchment (baseline is NO treatment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4a993e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_filtered = df_plot.query(' `Landuse intensity`==\"LOW\" & Soil==\"MIN\" ')\n",
    "gkk_no_treatment = df_filtered.query('Treatment==\"NO\"').groupby(['type'])['co2_net'].mean()\n",
    "gkk_prim_treatment = df_filtered.query('Treatment==\"PRIM\"').groupby(['type'])['co2_net'].mean()\n",
    "gkk_sec_treatment = df_filtered.query('Treatment==\"SEC\"').groupby(['type'])['co2_net'].mean()\n",
    "gkk_ter_treatment = df_filtered.query('Treatment==\"TER\"').groupby(['type'])['co2_net'].mean()\n",
    "print(\"Emission reductions thanks to primary WW treament:\")\n",
    "print((gkk_no_treatment - gkk_prim_treatment) / gkk_no_treatment * 100)\n",
    "print(\"Emission reductions thanks to secondary WW treatment\")\n",
    "print((gkk_no_treatment - gkk_sec_treatment) / gkk_no_treatment * 100)\n",
    "print(\"Emission reductions thanks to tertriary WW treatment\")\n",
    "print((gkk_no_treatment - gkk_ter_treatment) / gkk_no_treatment * 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2693092",
   "metadata": {},
   "source": [
    "### 3b. Plot distributions of emissions from irrigation and non-irrigation (hydroelectric and multipurpose) reservoirs to see if there's a difference in emissions between those three categories of reservoirs\n",
    "\n",
    "The data is plotted for LOW Landuse Intensity , Mineral Soil type and Primary Treatment, which are representative for the local conditions in Myanmar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a2d7ee1",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_plot['total'] = df_plot['co2_net'] + df_plot['ch4_net']\n",
    "df_plot['ch4_share'] = df_plot['ch4_net'] / df_plot['total'] * 100\n",
    "# Filter the dataframe to remove repeated entries\n",
    "df_filt = df_plot.query(' `Landuse intensity`==\"LOW\" & Soil==\"MIN\" ')\n",
    "df_irr_prim_treatment = df_filt.query('Treatment==\"PRIM\"')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4.5))\n",
    "bin_width = 333.3333\n",
    "use_logscale = False\n",
    "hist_alpha = 0.1235\n",
    "column_name = 'total'\n",
    "add_title = False\n",
    "palette_colors = sns.color_palette(\"dark\")\n",
    "if add_title:\n",
    "    fig.suptitle(\n",
    "        r'Distributions of net total GHG emissions across reservoir types', \n",
    "        fontsize=14)\n",
    "df1 = df_irr_prim_treatment.query('type==\"irrigation\"')\n",
    "df2 = df_irr_prim_treatment.query('type==\"hydroelectric\"')\n",
    "df3 = df_irr_prim_treatment.query('type==\"multipurpose\"')\n",
    "ix1, ix2, ix3 = 0, 4, 2\n",
    "ix1, ix2, ix3 = 0,1,2\n",
    "p1 = sns.histplot(data=df1[column_name],\n",
    "             alpha=hist_alpha,\n",
    "             binwidth=bin_width,\n",
    "             edgecolor=\".3\",\n",
    "             log_scale=use_logscale,\n",
    "             linewidth=0.5,\n",
    "             label='Irrigation', fill=True, ax=ax, kde='true', \n",
    "             color=palette_colors[ix1], kde_kws={'cut': 5})\n",
    "p2 = sns.histplot(data=df2[column_name],\n",
    "             alpha=hist_alpha,\n",
    "             binwidth=bin_width,\n",
    "             log_scale=use_logscale,\n",
    "             edgecolor=\".3\",\n",
    "             linewidth=0.5,\n",
    "             label='Hydroelectric', fill=True, ax=ax, kde='true', \n",
    "             color=palette_colors[ix2], kde_kws={'cut': 5})\n",
    "p3 = sns.histplot(data=df3[column_name],\n",
    "             alpha=hist_alpha,\n",
    "             binwidth=bin_width,\n",
    "             log_scale=use_logscale,\n",
    "             edgecolor=\".3\",\n",
    "             linewidth=0.5,\n",
    "             label='Multipurpose', fill=True, ax=ax, kde='true', \n",
    "             color=palette_colors[ix3], kde_kws={'cut': 5}, line_kws={\"color\": \"black\"})\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "\n",
    "df_1_x, df_1_y = ax.lines[0].get_data()\n",
    "df_2_x, df_2_y = ax.lines[1].get_data()\n",
    "df_3_x, df_3_y = ax.lines[2].get_data()\n",
    "# use fill_between to demonstrate where the extracted curves lie\n",
    "ax.fill_between(df_1_x, 0, df_1_y, color=palette_colors[ix1], alpha=0.215)\n",
    "ax.fill_between(df_2_x, 0, df_2_y, color=palette_colors[ix2], alpha=0.215)\n",
    "ax.fill_between(df_3_x, 0, df_3_y, color=palette_colors[ix3], alpha=0.215)\n",
    "#sns.histplot(data=df_irr_no_treatment, x='total', hue='type', kde='true')\n",
    "l1 = ax.legend(title='Reservoir Type', fontsize=14, frameon=False, fancybox=True, shadow=True)\n",
    "ax.set_xlim(0,7_000)\n",
    "ax.xaxis.grid(False)\n",
    "ax.yaxis.grid(False)\n",
    "plt.xlabel(r'Total Net Emission, gCO$_{2e}$ m$^{-2}$ yr$^{-1}$', fontsize=16)\n",
    "plt.ylabel('Number of reservoirs', fontsize=16)\n",
    "plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "sns.despine(offset=15, trim=False)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"emission_histograms_per_reservoir_type.svg\"), format=\"svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b16eb029",
   "metadata": {},
   "source": [
    "## Plot emissions vs altitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d837a71b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the DataFrame into two based on 'soil type'\n",
    "df = df_irr_prim_treatment\n",
    "df_reg = df[df['fsl_masl']>300]\n",
    "df_hp = df[df['type']=='hydroelectric']\n",
    "df_multi = df[df['type']=='multipurpose']\n",
    "grid = sns.jointplot(\n",
    "    height=5,\n",
    "    data=df, x=\"fsl_masl\", y=\"total\", hue='type',\n",
    "    hue_order = ['hydroelectric', 'irrigation', 'multipurpose'],\n",
    "    joint_kws = {'edgecolor':'k'},\n",
    "    marginal_kws=dict(fill=False),\n",
    "    edgecolors='k',\n",
    "    ratio = 2,\n",
    "    palette=\"colorblind\", alpha=0.4, marginal_ticks=False)\n",
    "grid.plot_marginals(\n",
    "    sns.histplot, color='green', fill=True, alpha=0.3, kde=True, linewidth=0.5, edgecolor=\".3\",)\n",
    "grid.ax_joint.grid(False)\n",
    "sns.regplot(data=df_reg, x=\"fsl_masl\", y=\"total\", color='gray', ax=grid.ax_joint, scatter=False)\n",
    "# set the current axis to be the joint plot's axis\n",
    "plt.sca(grid.ax_joint)\n",
    "leg = plt.legend(title='Reservoir Type', fontsize=11, frameon=False)\n",
    "for lh in leg.legendHandles: \n",
    "    lh.set_alpha(0.4)\n",
    "new_labels = ['Hydroelectric', 'Irrigation', 'Multipurpose']\n",
    "for t, l in zip(leg.texts, new_labels):\n",
    "    t.set_text(l)\n",
    "sns.set_style(\"whitegrid\", {'axes.grid' : False})\n",
    "grid.ax_joint.set_ylim([0,7_000])\n",
    "grid.ax_joint.set_xlim([0,1_200]) \n",
    "grid.ax_marg_x.remove()\n",
    "grid.ax_marg_y.set_axis_off()\n",
    "plt.ylabel(r'Total Net Emission, gCO$_{2e}$/m$^{2}$/yr', fontsize=15)\n",
    "plt.xlabel('Elevation at FSL, m.a.s.l.', fontsize=15)\n",
    "plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.setp(grid.ax_joint.get_legend().get_texts(), fontsize='13') # for legend text\n",
    "plt.setp(grid.ax_joint.get_legend().get_title(), fontsize='13') # for legend title\n",
    "sns.despine(offset=10, trim=False)\n",
    "fig.tight_layout()\n",
    "plt.savefig(\"figures/ghg_visualisation/emission_elevation_plot.svg\", format='svg')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "829d421e",
   "metadata": {},
   "source": [
    "## Piece-wise distance matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48504a2c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Calculate distance matrices\n",
    "from geopy.distance import geodesic\n",
    "num_points = len(df)\n",
    "emm_dist_matrix = np.zeros((num_points, num_points))\n",
    "co2_emm_dist_matrix = np.zeros((num_points, num_points))\n",
    "ch4_emm_dist_matrix = np.zeros((num_points, num_points))\n",
    "geom_dist_matrix = emm_dist_matrix.copy()\n",
    "prop_factor = 0\n",
    "\n",
    "for i in range(num_points):\n",
    "    for j in range(num_points):\n",
    "        emm_dist_matrix[i, j] = abs(df['total'].iloc[i] - df['total'].iloc[j])\n",
    "        co2_emm_dist_matrix[i, j] = abs(df['co2_net'].iloc[i] - df['co2_net'].iloc[j])\n",
    "        ch4_emm_dist_matrix[i, j] = abs(df['ch4_net'].iloc[i] - df['ch4_net'].iloc[j])\n",
    "        geom_dist_matrix[i,j] = geodesic(\n",
    "            (df['coordinates_0'].iloc[i], df['coordinates_1'].iloc[i]),\n",
    "            (df['coordinates_0'].iloc[j], df['coordinates_1'].iloc[j])).meters/1_000 + \\\n",
    "        prop_factor * abs(df['fsl_masl'].iloc[i] - df['fsl_masl'].iloc[j])\n",
    "    \n",
    "# Put the distance matrices data into a dataframe, filter, cut into bins and plot\n",
    "conn_type_arr = np.empty((num_points, num_points), dtype='object')\n",
    "conn_type_arr.fill(\"\")\n",
    "# create short type names\n",
    "type_mapping = {'irrigation': 'irr', 'hydroelectric': 'hp', 'multipurpose': 'multi'}\n",
    "# Add a new column 'type_short' using the mapping\n",
    "df['type_short'] = df['type'].map(type_mapping)\n",
    "\n",
    "for i in range(num_points):\n",
    "    for j in range(num_points):\n",
    "        conn_type = (sorted(tuple((df['type_short'].iloc[i], df['type_short'].iloc[j]))))\n",
    "        conn_type_str = '-'.join(conn_type)\n",
    "        conn_type_arr[i,j] = conn_type_str\n",
    "        \n",
    "mid_point = int(np.floor(len(geom_dist_matrix.flatten())/2))\n",
    "df_dist = pd.DataFrame({\n",
    "    'Distance': geom_dist_matrix.flatten()[:mid_point],\n",
    "    'GHG Difference': emm_dist_matrix.flatten()[:mid_point],\n",
    "    'CO2 emission Difference': co2_emm_dist_matrix.flatten()[:mid_point],\n",
    "    'CH4 emission Difference': ch4_emm_dist_matrix.flatten()[:mid_point],\n",
    "    'Res Types': conn_type_arr.flatten()[:mid_point]\n",
    "})\n",
    "df_dist = df_dist.query(\"Distance < 100 & Distance != 0\")\n",
    "\n",
    "df_dist['Distance Category'] = pd.cut(\n",
    "    df_dist['Distance'], [0, 25, 75, 10000], labels = ['$\\mathrm{<}25$', '$25-75$', '$\\mathrm{>}75$'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5552191",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_dist.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b366f6e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.legend_handler import HandlerTuple\n",
    "fig, ax = plt.subplots(figsize=(5, 3.5))\n",
    "palette='colorblind'\n",
    "# make grouped boxplot and save it in a variable\n",
    "bp = sns.boxplot(y='GHG Difference', x='Distance Category', \n",
    "                 ax=ax,\n",
    "                 data=df_dist, \n",
    "                 fliersize=0,\n",
    "                 showfliers=False,\n",
    "                 showmeans=False,\n",
    "                 linewidth=0.5,\n",
    "                 palette=palette, notch=False,\n",
    "                 boxprops = dict(alpha=0.4, linestyle='-', linewidth=0.5),\n",
    "                 hue='Res Types')\n",
    " \n",
    "# make grouped stripplot and save it in a variable\n",
    "sns.stripplot(y='GHG Difference', x='Distance Category', \n",
    "                   data=df_dist, \n",
    "                   jitter=True,\n",
    "                   dodge=True, \n",
    "                   marker='o', \n",
    "                   size=4,\n",
    "                   linewidth=0.5,\n",
    "                   edgecolor='k',\n",
    "                   ax=ax,\n",
    "                   alpha=0.3,\n",
    "                   hue='Res Types')\n",
    "sns.despine(offset=15, trim=False)\n",
    "# get legend information from the plot object\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "n_handles = 6\n",
    "h_items = [(handles[i], handles[n_handles+i]) for i in range(0,n_handles)]\n",
    "ll = ax.legend(handles=handles,\n",
    "          labels=['irr-irr', 'hp-irr', 'irr-multi', 'hp-hp', 'hp-multi','multi-multi'],\n",
    "          ncol=2,\n",
    "          columnspacing=1,\n",
    "          alignment='left',\n",
    "          edgecolor='k',\n",
    "          handlelength=2.5,\n",
    "          handleheight=1,\n",
    "          title='Reservoir types', fontsize=10, frameon=True, fancybox=False, shadow=False, \n",
    "          framealpha=0.70,\n",
    "          loc=\"upper left\",\n",
    "          handler_map={tuple: HandlerTuple(ndivide=None)})\n",
    "ll.get_frame().set_linewidth(0.5)\n",
    "ll.get_frame().set_edgecolor(None)\n",
    "for i in range(0,6):\n",
    "    ll.legend_handles[i].set_alpha(0.3)\n",
    "    ll.legend_handles[i].set_width(22)\n",
    "    ll.legend_handles[i].set_height(8)\n",
    "    ll.legend_handles[i].set_linewidth(0.5)\n",
    "    ll.legend_handles[i].set(edgecolor='k', linewidth=0.7)\n",
    "plt.setp(ll.get_texts(), fontsize='13') # for legend text\n",
    "plt.setp(ll.get_title(), fontsize='13') # for legend title\n",
    "ax.set_ylabel(r'Emission Difference, gCO$_{2e}$/m$^{2}$/yr', fontsize=15)\n",
    "ax.set_xlabel('Distance between dams, km', fontsize=15)\n",
    "ax.set_ylim([0,6_000])\n",
    "plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "plt.savefig(\"figures/ghg_visualisation/distance_plot.svg\", format='svg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a64dfe60",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "geom_flat = geom_dist_matrix.flatten()\n",
    "emm_flat = emm_dist_matrix.flatten()\n",
    "mid_point = int(np.floor(len(geom_dist_matrix.flatten())/2))\n",
    "fig = plt.figure(figsize=(5,3))\n",
    "ax = plt.scatter(geom_flat[:mid_point],emm_flat[:mid_point], alpha=0.2, color='k', marker='o', s=5)\n",
    "plt.xlabel('Distance between reservoirs, km')\n",
    "plt.ylabel(\"Difference between total net emissions\")\n",
    "plt.xlim([0,100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60bfc01c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the DataFrame into two based on 'soil type'\n",
    "df = df_irr_prim_treatment\n",
    "fig, ax = plt.subplots(figsize=(5, 4))\n",
    "sns.set_style(\"whitegrid\", {'axes.grid' : False})\n",
    "fig.suptitle(\n",
    "    r'Share of methane emissions in total emissions of a reservoir', \n",
    "    fontsize=14)\n",
    "ax = sns.boxenplot(data=df, y=\"ch4_share\", x=\"type\")\n",
    "ax.xaxis.grid(False)\n",
    "ax.yaxis.grid(False)\n",
    "#l1 = ax.legend(\n",
    "#    title='Intake Type', fontsize=13, frameon=False, fancybox=True, shadow=True, loc=\"center right\")\n",
    "plt.ylabel(r'CH$_4$ to total emission ratio, \\%', fontsize=14)\n",
    "plt.xlabel('Reservoir Type', fontsize=14)\n",
    "plt.tick_params(axis='both', which='major', labelsize=13)\n",
    "sns.despine(offset=15, trim=False)\n",
    "ax.set_ylim(0, 100)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"ch4_emission_share_distributions.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e280cd9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Get the ch4 and co2 emission statistics\n",
    "summary_stats = df_irr_prim_treatment.groupby('type')['ch4_share'].describe()\n",
    "summary_stats.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32223e6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_irr_prim_treatment_filt = df_irr_prim_treatment[df_irr_prim_treatment['ch4_net']<12_000]\n",
    "summary_stats = df_irr_prim_treatment_filt.groupby('type')['ch4_net'].describe()\n",
    "summary_stats.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f791dde",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_stats = df_irr_prim_treatment.groupby('type')['co2_net'].describe()\n",
    "summary_stats.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "920c3b19",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_stats = df_irr_prim_treatment.groupby('type')['total'].describe()\n",
    "summary_stats.head(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0f1903d",
   "metadata": {},
   "source": [
    "### 3c. Quantify the difference in CH4 emissions for reservoirs with deep and shallow water uptakes\n",
    "\n",
    "The data is plotted for Mineral soil type, Low landuse intensity and primary level of wastewater treatment in the catchment"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46084e11",
   "metadata": {},
   "source": [
    "#### <span style=\"color: red;\">TODO: QUANTIFY THE LOSS OF VOLUME REQUIRED TO OPERATE THE RESERVOIRS WITHOUT TOUCHING THE HYPOLIMNION</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6413e042",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Merged the dataframe with shallow and deep intake data\n",
    "merged_df_deep_intake = merged_df.copy().drop(\"ch4_net_nodegassing\", axis=1)\n",
    "merged_df_deep_intake['Intake'] = 'deep'\n",
    "merged_df_shallow_intake = (merged_df.copy()\n",
    "    .drop(\"ch4_net\", axis=1)\n",
    "    .rename(columns={\"ch4_net_nodegassing\": \"ch4_net\"}))\n",
    "merged_df_shallow_intake['Intake'] = 'shallow'\n",
    "df_combined = pd.concat([merged_df_deep_intake, merged_df_shallow_intake], ignore_index=True)\n",
    "df_combined['total'] = df_combined['ch4_net'] + df_combined['co2_net']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87fde7de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the DataFrame into two based on 'soil type'\n",
    "df = df_combined.query('Soil == \"MIN\" & `Landuse intensity`==\"LOW\" & Treatment==\"PRIM\"')\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "fig.suptitle(\n",
    "    r'Total net GHG emissions per reservoir and intake type', \n",
    "    fontsize=14)\n",
    "ax = sns.boxenplot(data=df, x=\"total\", y=\"type\", hue=\"Intake\")\n",
    "ax.xaxis.grid(False)\n",
    "ax.yaxis.grid(False)\n",
    "l1 = ax.legend(\n",
    "    title='Intake Type', fontsize=13, frameon=False, fancybox=True, shadow=True, loc=\"center right\")\n",
    "plt.xlabel(r'Total Emission, gCO$_{2,eq}$ m$^{-2}$ yr$^{-1}$', fontsize=14)\n",
    "#plt.ylabel('Reservoir type', fontsize=14)\n",
    "plt.ylabel(None)\n",
    "plt.tick_params(axis='both', which='major', labelsize=13)\n",
    "sns.despine(offset=15, trim=False)\n",
    "ax.set_xlim(0, 8_000)\n",
    "fig.tight_layout()\n",
    "fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"total_net_emission_per_res_and_intake_type.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8b0dbd0",
   "metadata": {},
   "source": [
    "### 3d. Plot the distribution of percentage reductions of total emission when shallow vs deep intakes are employed\n",
    "#### <span style=\"color: red;\">TODO: Plot that versus loss of volume AND find out which reservoirs have the highest potential for emission reduction</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d3c29f",
   "metadata": {},
   "outputs": [],
   "source": [
    "deep = df.query('Intake==\"deep\"').reset_index()\n",
    "shallow = df.query('Intake==\"shallow\"').reset_index()\n",
    "deep['reduction'] = (deep['total']-shallow['total'])/deep['total'] * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccbd93ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "with sns.axes_style(\"white\"):\n",
    "    g = sns.FacetGrid(deep, hue=\"type\", col=\"type\", height=4)\n",
    "g.fig.suptitle(\n",
    "    r'Potential reduction in GHG emissions by employing shallow water uptakes', \n",
    "    fontsize=14)\n",
    "g.map(hexbin, \"reduction\", 'coordinates_0', extent=[0, 100, 10, 30]); \n",
    "g.set(ylabel=\"Latitude, deg\")\n",
    "g.set(xlabel=\"% reduction in total GHG emissions\")\n",
    "axes = g.axes.flat\n",
    "axes[0].set(title=\"Irrigation Reservoirs\")\n",
    "axes[1].set(title=\"Hydroelectric Reservoirs\")\n",
    "axes[2].set(title=\"Multipurpose Reservoirs\")\n",
    "g.fig.tight_layout()\n",
    "g.fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"red_by_shallow_water_uptakes.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afcd35a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "with sns.axes_style(\"white\"):\n",
    "    g = sns.FacetGrid(deep, hue=\"type\", col=\"type\", height=4)\n",
    "g.fig.suptitle(\n",
    "    r'Distribution of total net GHG emissions by reservoir type and latitude', \n",
    "    fontsize=14)\n",
    "g.map(hexbin, \"total\", 'coordinates_0', extent=[0, 6_000, 10, 30])\n",
    "g.set(ylabel=\"Latitude, deg\")\n",
    "g.set(xlabel=r'Total Emission, gCO$_{2,eq}$ m$^{-2}$ yr$^{-1}$')\n",
    "axes = g.axes.flat\n",
    "axes[0].set(title=\"Irrigation Reservoirs\")\n",
    "axes[1].set(title=\"Hydroelectric Reservoirs\")\n",
    "axes[2].set(title=\"Multipurpose Reservoirs\")\n",
    "g.fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"distribution_of_ghg_emissions_res_type_lat.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7e767b0",
   "metadata": {},
   "source": [
    "### 3e. Plot distributions of reductions on histograms with kde plots on top of them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30b57ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "palette_colors = sns.color_palette(\"Set2\") \n",
    "fig.suptitle(\n",
    "    r'Reductions in total GHG emissions via shallow water uptakes', \n",
    "    fontsize=14)\n",
    "df1 = deep.query('type==\"hydroelectric\"')\n",
    "df2 = deep.query('type==\"multipurpose\"')\n",
    "df3 = deep.query('type==\"irrigation\"')\n",
    "ix1, ix2, ix3 = 0, 4, 2\n",
    "p1 = sns.distplot(a=df1['reduction'],\n",
    "             rug=False,\n",
    "             label='Hydroelectric', \n",
    "             ax=ax, kde='true', \n",
    "             color=palette_colors[ix1], kde_kws={'cut': 5})\n",
    "p2 = sns.distplot(a=df2['reduction'],\n",
    "             rug=False,\n",
    "             label='Multipurpose',  \n",
    "             ax=ax, kde='true', \n",
    "             color=palette_colors[ix2], kde_kws={'cut': 5})\n",
    "p3 = sns.distplot(a=df3['reduction'],\n",
    "             rug=False,\n",
    "             label='Irrigation',  \n",
    "             ax=ax, kde='true', \n",
    "             color=palette_colors[ix3], kde_kws={'cut': 5})\n",
    "\n",
    "sns.set_theme(style=\"whitegrid\")\n",
    "\n",
    "df_1_x, df_1_y = ax.lines[0].get_data()\n",
    "df_2_x, df_2_y = ax.lines[1].get_data()\n",
    "df_3_x, df_3_y = ax.lines[2].get_data()\n",
    "# use fill_between to demonstrate where the extracted curves lie\n",
    "ax.fill_between(df_1_x, 0, df_1_y, color=palette_colors[ix1], alpha=0.32)\n",
    "ax.fill_between(df_2_x, 0, df_2_y, color=palette_colors[ix2], alpha=0.32)\n",
    "ax.fill_between(df_3_x, 0, df_3_y, color=palette_colors[ix3], alpha=0.32)\n",
    "#sns.histplot(data=df_irr_no_treatment, x='total', hue='type', kde='true')\n",
    "l1 = ax.legend(title='Reservoir Type', fontsize=13, frameon=False, fancybox=True, shadow=True,\n",
    "              loc=\"upper left\")\n",
    "ax.set_xlim(0, 100)\n",
    "ax.xaxis.grid(False)\n",
    "ax.yaxis.grid(False)\n",
    "plt.xlabel(r'Reduction in total emission %', fontsize=14)\n",
    "plt.ylabel('Density', fontsize=14)\n",
    "plt.tick_params(axis='both', which='major', labelsize=14)\n",
    "sns.despine(offset=15, trim=False)\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig(os.path.join(\n",
    "    \"figures\",\"ghg_visualisation\",\"reductions_gas_emissions.svg\"), format=\"svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f4c1eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find statistics\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(\"GHG emission reduction in hydroelectric reservoirs\")\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(deep.query('type==\"hydroelectric\"')['reduction'].describe())\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(\"GHG emission reduction in irrigation reservoirs\")\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(deep.query('type==\"irrigation\"')['reduction'].describe())\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(\"GHG emission reduction in multipurpose reservoirs\")\n",
    "print(\"-----------------------------------------------------------\")\n",
    "print(deep.query('type==\"multipurpose\"')['reduction'].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4814dd4",
   "metadata": {},
   "source": [
    "### 3f. Plot all emissions across all reservoirs dividing the results into HP and non-HP domain\n",
    "\n",
    "#### <span style=\"color: red;\">TODO: Plot the values in sorted order and create maps</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eef760d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_emission_split(\n",
    "        df: pd.DataFrame, col_indices, file_name: str | None = None,\n",
    "        xlimits = ((1000,7000), (100, 800), (1000,7000)),\n",
    "        height=20, aspect=.2):\n",
    "    \"\"\" \"\"\"\n",
    "    x_vars = []\n",
    "    sns.set_theme(style=\"whitegrid\")\n",
    "    g = sns.PairGrid(df,\n",
    "                     x_vars=df.columns[col_indices], y_vars=[\"Name\"], hue='Treatment',\n",
    "                     height=height, aspect=aspect)\n",
    "    g.map(sns.stripplot, size=10, orient=\"h\", jitter=True, alpha=0.8,\n",
    "          linewidth=1, edgecolor=\"w\", palette=\"flare_r\") #palette=[\"m\", \"g\", \"r\", \"y\"])\n",
    "    # Use the same x axis limits on all columns and add better labels\n",
    "    g.set(xlabel=r\"Emission, gCO$_{2,eq}$ m$^{-2}$ yr$^{-1}$\", ylabel=\"\") #, xscale='log')\n",
    "    g = g.add_legend()\n",
    "    # Use semantically meaningful titles for the columns\n",
    "    titles = [\n",
    "        \"Total emission\", \n",
    "        \"CO$_{2}$ net emission\", \n",
    "        \"CH$_{4}$ net emission\"]\n",
    "    axes = g.axes.flat\n",
    "    axes[0].set_xlim(*xlimits[0])\n",
    "    axes[1].set_xlim(*xlimits[1])\n",
    "    axes[2].set_xlim(*xlimits[2])\n",
    "    for ax, title in zip(g.axes.flat, titles):\n",
    "        # Set a different title for each axes\n",
    "        ax.set(title=title)\n",
    "        # Make the grid horizontal instead of vertical\n",
    "        ax.xaxis.grid(False)\n",
    "        ax.yaxis.grid(True)\n",
    "        \n",
    "    sns.despine(left=True, bottom=True)\n",
    "    plt.savefig(file_name, format=\"svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0cdab0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate total emission\n",
    "df_plot['total'] = df_plot['co2_net'] + df_plot['ch4_net']\n",
    "# Filter the dataframe to remove repeated entries\n",
    "df_filt = df_plot.query(' `Landuse intensity`==\"LOW\" & Soil==\"MIN\" ')\n",
    "col_indices = df_plot.columns.get_indexer(['total', 'co2_net', 'ch4_net'])\n",
    "df_1 = df_filt.query('type==\"irrigation\"')\n",
    "df_sorted = df_1.sort_values(\"total\", ascending=False)\n",
    "\n",
    "plot_emission_split(\n",
    "    df_sorted.iloc[:int(len(df_sorted)/2)+4,:], \n",
    "    col_indices, \n",
    "    os.path.join(\"figures\",\"ghg_visualisation\",\"tot_irr_plot_irrig_1.svg\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee097cc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_emission_split(\n",
    "    df_sorted.iloc[(int(len(df_sorted)/2)+5):,:], \n",
    "    col_indices, \n",
    "    os.path.join(\"figures\",\"ghg_visualisation\",\"tot_irr_plot_irrig_2.svg\"),\n",
    "    xlimits = ((0,2500), (0, 700), (0,2500)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5dbfceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2 = df_filt.query('type==\"hydroelectric\" | type==\"multipurpose\"')\n",
    "df_sorted = df_2.sort_values(\"total\", ascending=False)\n",
    "plot_emission_split(\n",
    "    df_sorted, \n",
    "    col_indices, \n",
    "    os.path.join(\"figures\",\"ghg_visualisation\",\"tot_irr_plot_hp_multi_1.svg\"),\n",
    "    xlimits = ((0,4000), (0, 700), (0,4000)),\n",
    "    height = 20, aspect=.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb227fe9",
   "metadata": {},
   "source": [
    "### END"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}