Notebook_12_run_MOO_and_visualise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d19b77a3",
   "metadata": {},
   "source": [
    "## Visualise results from the selection of dam porfolios using multiobjective optimization\n",
    "### T. Janus\n",
    "### 15/01/2024"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f010746f",
   "metadata": {},
   "source": [
    "## TODO:\n",
    "1. Visualise different dam scenarios on maps (borrow the maps from one of the previous notebooks\n",
    "2. Create a composite figure with tiles using facetgrid etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ca8e2b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from dataclasses import dataclass, field\n",
    "from typing import ClassVar, Dict, List, Any, Tuple, Set, Tuple, Sequence\n",
    "from typing import TypeAlias, TypeVar, Generic\n",
    "import subprocess\n",
    "import pathlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import string\n",
    "from datetime import datetime\n",
    "from parse import parse\n",
    "from ttp import ttp\n",
    "import json\n",
    "import gc\n",
    "import bson\n",
    "import pprint\n",
    "import re\n",
    "import ast\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go\n",
    "import networkx as nx\n",
    "\n",
    "import pygmo as pg\n",
    "from tqdm import tqdm\n",
    "\n",
    "import seaborn as sns # for Data visualization\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt # for Data visualization\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "from jinja2 import Template\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "T = TypeVar(\"T\")\n",
    "GenericCollection: TypeAlias = Set[T] | Tuple[T] | List[T]\n",
    "NumType= TypeVar('NumType', int, float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1f68912",
   "metadata": {},
   "outputs": [],
   "source": [
    "def reduce_mem_usage(df: pd.DataFrame, verbose: bool = True) -> pd.DataFrame:\n",
    "    \"\"\"Reduce memory taken up by a dataframe by setting and downgrading the data types of columns\n",
    "    based on their content.\n",
    "    Allows to load larger dataframes for a given size of memory without chunking or using external\n",
    "    libraries such as Dask, Ray, Modin or Vaex.\"\"\"\n",
    "    start_mem = df.memory_usage().sum() / 1024**2\n",
    "    print('Memory usage of dataframe is {:.2f} MB'.format(start_mem))\n",
    "\n",
    "    for col in df.columns:\n",
    "        col_type = df[col].dtype\n",
    "    if col_type != object:\n",
    "            c_min = df[col].min()\n",
    "            c_max = df[col].max()\n",
    "            if str(col_type)[:3] == 'int':\n",
    "                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n",
    "                    df[col] = df[col].astype(np.int8)\n",
    "                elif c_min > np.iinfo(np.uint8).min and c_max < np.iinfo(np.uint8).max:\n",
    "                    df[col] = df[col].astype(np.uint8)\n",
    "                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n",
    "                    df[col] = df[col].astype(np.int16)\n",
    "                elif c_min > np.iinfo(np.uint16).min and c_max < np.iinfo(np.uint16).max:\n",
    "                    df[col] = df[col].astype(np.uint16)\n",
    "                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n",
    "                    df[col] = df[col].astype(np.int32)\n",
    "                elif c_min > np.iinfo(np.uint32).min and c_max < np.iinfo(np.uint32).max:\n",
    "                    df[col] = df[col].astype(np.uint32)\n",
    "                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n",
    "                    df[col] = df[col].astype(np.int64)\n",
    "                elif c_min > np.iinfo(np.uint64).min and c_max < np.iinfo(np.uint64).max:\n",
    "                    df[col] = df[col].astype(np.uint64)\n",
    "            elif str(col_type)[:5] == 'float':\n",
    "                if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n",
    "                    df[col] = df[col].astype(np.float16)\n",
    "                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n",
    "                    df[col] = df[col].astype(np.float32)\n",
    "                else:\n",
    "                    df[col] = df[col].astype(np.float64)\n",
    "\n",
    "    end_mem = df.memory_usage().sum() / 1024**2\n",
    "    if verbose:\n",
    "        print('Memory usage after optimization is: {:.2f} MB'.format(end_mem))\n",
    "        print('Decreased by {:.1f}%'.format(100 * (start_mem - end_mem) / start_mem))\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ebdc95",
   "metadata": {},
   "outputs": [],
   "source": [
    "#matplotlib.rcParams['font.family'] = ['monospace', 'sans-serif']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9fcae9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxiliary functions\n",
    "def read_id_ifc_map(file_name: pathlib.Path = pathlib.Path('outputs/moo/id_to_ifc.json')):\n",
    "    \"\"\"Read the dictionary showing mapping between optimization dam IDs and IFC dam IDs\"\"\"\n",
    "    with open(file_name, 'r') as file:\n",
    "        _id_ifc_map = json.load(file)\n",
    "        return {int(key) : value for key, value in _id_ifc_map.items()}\n",
    "    \n",
    "def set_remap(\n",
    "        value_list: Set[NumType], value_map, missing_val_id: int = -99, \n",
    "        safe: bool = False) -> Set[NumType]:\n",
    "    \"\"\"Map values in a set to new values using a dictionary given in `value_map`\"\"\"\n",
    "    if safe:\n",
    "        id_set = set([value_map.get(value, missing_val_id) for value in value_list])\n",
    "        return id_set\n",
    "    return {value_map[value] for value in value_list}\n",
    "\n",
    "def get_every_n_row(df: pd.DataFrame, n: int) -> pd.DataFrame:\n",
    "    \"\"\"Return every n-th row of a dataframe\n",
    "    Used for 'skimming-down' large chunks of data to help with post-processing\n",
    "    and prototyping visualisation\"\"\"\n",
    "    return df.iloc[::n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdcff9a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class SolutionFileParser:\n",
    "    \"\"\"Class for parsing optimization solution .sol files returned from the optimization\n",
    "    algorithm in `moo_solver_CPAIOR`.\n",
    "    Uses TTP - a Python module that allows relatively fast performance parsing of \n",
    "        semi-structured text data using templates.\n",
    "    \"\"\"\n",
    "    file_path: str | pathlib.Path\n",
    "    # Template of the header text\n",
    "    header_template: ClassVar[str] = \"\"\"\n",
    "Date/time: {{ date_time }}\n",
    "Data file: {{ file_name }}\n",
    "Wall time: {{ wall_time | to_float}} seconds.\n",
    "CPU time: {{ cpu_time | to_float}} seconds.\n",
    "seed: {{ seed | to_int}}\n",
    "num_solutions: {{ num_solutions | to_int}}\n",
    "# pruning steps (# nodes): {{ num_pruning_steps | to_int}}\n",
    "Max policies considered: {{ max_policies | to_int}}\n",
    "Policies considered: {{ num_policies | to_int}}\n",
    "Pruned policies: {{ pruned_policies | to_int}}\n",
    "epsilon: {{ epsilon | to_float}}\n",
    "batch size: {{ batch_size | to_int }}\n",
    "criteria: {{ criteria | ORPHRASE | split(',')}}\n",
    "\"\"\"\n",
    "    # Template for the solutions part of the .sol file \n",
    "    # (generated dynamically by `_create_solution_template` method)\n",
    "    solution_template: str = \"\"\n",
    "    header: str = \"\"\n",
    "    solutions: str = \"\"\n",
    "    data: Dict[str, Any] = field(default_factory=dict)\n",
    "    \n",
    "    def __post_init__(self) -> None:\n",
    "        \"\"\"Load data from the file which path is provided as an attribute.\"\"\"\n",
    "        with open(self.file_path, 'r') as file:\n",
    "            raw_data = file.readlines()\n",
    "            self.header = \"\".join(raw_data[0:13])\n",
    "            self.solutions = \"\".join(raw_data[14:])\n",
    "            \n",
    "    def _create_solution_template(self) -> None:\n",
    "        \"\"\"Dynamically creates solution template for the ttp package.\"\"\"\n",
    "        sol_template = \"\"\n",
    "        crit_template = \"{{{{ {} | to_float }}}}\"\n",
    "        num_dams_template = \"{{{{ {} | to_int }}}}\"\n",
    "        for criterion in self.data['header']['criteria']:\n",
    "            sol_template += crit_template.format(criterion) + \", \"\n",
    "        sol_template += num_dams_template.format('num_dams') + \", \"\n",
    "        sol_template += \"{{ dam_ids | ORPHRASE | split(' ')}}\"\n",
    "        self.solution_template = sol_template\n",
    "        \n",
    "    def parse(self) -> None:\n",
    "        \"\"\"Parse the file contents given the contents of the sections of the sol file,\n",
    "        i.e. header and solutions and their templates. Note that the template for the\n",
    "        header is manually input while the template for the solution section is\n",
    "        generated dynamically.\n",
    "        The parsed file contents is then stored in `self.data` dictionary\n",
    "        \"\"\"\n",
    "        # 1. Parse header\n",
    "        parser_header = ttp(self.header, self.header_template)\n",
    "        parser_header.parse()\n",
    "        header_data = parser_header.result(structure=\"flat_list\")[0]\n",
    "        header_data['criteria'] = [criterion.strip() for criterion in header_data['criteria']]\n",
    "        date_time_formatted = re.sub(r'_+', ',', header_data['date_time'])\n",
    "        header_data['date_time'] = datetime.strptime(\n",
    "            date_time_formatted, \"%a,%b,%d,%H,%M,%S,%Y\")\n",
    "        self.data['header'] = header_data\n",
    "        # 2. Dynamically create a solution template\n",
    "        self._create_solution_template()\n",
    "        # 3. Parse solutions\n",
    "        parser_sol = ttp(self.solutions, self.solution_template)\n",
    "        parser_sol.parse()\n",
    "        solution_data = parser_sol.result(structure=\"flat_list\")\n",
    "        for solution in solution_data:\n",
    "            solution['dam_ids'] = [int(dam_id) for dam_id in solution['dam_ids']]\n",
    "        self.data['solutions'] = solution_data\n",
    "        \n",
    "    @property\n",
    "    def solutions_df(self) -> pd.DataFrame:\n",
    "        \"\"\"Converts solutions into a dataframe and compute additional variables: \n",
    "            1. land_loss := [lost forest area] + [lost agricultural area]\n",
    "            2. ghg_intensity := [ghg emissions] / [energy produced] * [unit conversion factor]\n",
    "        Returns a dataframe with new columns\n",
    "        \"\"\"\n",
    "        df = pd.DataFrame(self.data['solutions'])\n",
    "        if set(['loss_agri', 'loss_forest']).issubset(set(df.columns)):\n",
    "            df['land_loss'] = df['loss_agri'] + df['loss_forest']\n",
    "        # Calculate ghg intensity, for ghg in tonneCO2eq/year and energy in MWh/d\n",
    "        # GHG intensity needs to be in gCO2eq/kWh\n",
    "        df['ghg_intensity'] = df['ghg'] / df['energy'] * 1_000 / 365.25 / 24\n",
    "        return reduce_mem_usage(df)\n",
    "        \n",
    "        \n",
    "    def to_json(self, json_file: str | pathlib.Path) -> None:\n",
    "        \"\"\"Save output data to a json file.\"\"\"  \n",
    "        def serialize_datetime(obj: Any) -> str:\n",
    "            \"\"\"Custom serialization function for datetime objects\"\"\"\n",
    "            if isinstance(obj, datetime):\n",
    "                return obj.isoformat()\n",
    "        # Save to json using custom datetime object serialization\n",
    "        with open(json_file, 'w') as file:\n",
    "            json_string = json.dumps(self.data, default=serialize_datetime, indent=4)\n",
    "            file.write(json_string)\n",
    "            \n",
    "    def to_bson(self, file_path: str | pathlib.Path) -> None:\n",
    "        \"\"\"Save output date to a bson (binary) file.\"\"\"\n",
    "        with open(file_path, 'wb') as bson_file:\n",
    "            serialized_data = bson.dumps(self.data)\n",
    "            bson_file.write(serialized_data)\n",
    "            \n",
    "    def to_csv(self, csv_file: str | pathlib.Path) -> None:\n",
    "        \"\"\"Save solutions to a csv file.\n",
    "        Dam IDs are saved as a string representation of a list of integers. To retrieve\n",
    "        the list of interegers, parse the dam_ids column\n",
    "        df_read.dam_ids = df_read.dam_ids.map(ast.literal_eval) \n",
    "        (requires `import ast`)\"\"\"\n",
    "        self.solutions_df.to_csv(csv_file, encoding='utf-8', index=False)\n",
    "\n",
    "\n",
    "@dataclass\n",
    "class OutputVisualiser:\n",
    "    \"\"\"Visualisation of optimization results from `moo_solver_CPAIOR`\n",
    "    Attributes:\n",
    "        data: pd.DataFrame Optimization outputs\n",
    "    \"\"\"\n",
    "    data: pd.DataFrame\n",
    "        \n",
    "    @property\n",
    "    def columns(self) -> List[str]:\n",
    "        \"\"\"Return the list of column names\"\"\"\n",
    "        return list(self.data.columns)\n",
    "    \n",
    "    def plot_parallel(\n",
    "            self, columns: GenericCollection[str], color_col: str, title: str | None = None,\n",
    "            color_limits: Tuple[float, float] = (0,200), **kwargs) -> go.Figure:\n",
    "        \"\"\"Plot a parallel coordinate plot using plotly\n",
    "        Args:\n",
    "            columns: list of columns (variables) to be plotted\n",
    "            color_col: column which values should be associated with line color\n",
    "            title: plot title\n",
    "            color_limits: tuple with minimum and maximum value associated with the minimum\n",
    "                and the maximum color value in the color pallette\n",
    "            kwargs: Additional keyword arguments for plotly express `parallel_coordinates`.\n",
    "        Returns:\n",
    "            plotly.graph_objects.Figure object\n",
    "        \"\"\"\n",
    "        # Other scales:\n",
    "        # px.colors.diverging.Tealrose\n",
    "        # px.colors.sequential.Blues\n",
    "        # px.colors.sequential.Oranges\n",
    "        # px.colors.diverging.RdYlBu\n",
    "        # color_continuous_scale=px.colors.diverging.Armyrose\n",
    "        # color_continuous_midpoint=2\n",
    "        fig = px.parallel_coordinates(self.data, color=color_col, dimensions=columns,\n",
    "                              color_continuous_scale=px.colors.diverging.Tealrose, width=1000,\n",
    "                              title=title, range_color=color_limits, **kwargs)\n",
    "        fig.show()\n",
    "        return fig\n",
    "        \n",
    "    def plot_scatter_2D(\n",
    "            self, x_col: str, y_col: str, hue: str | None = None, \n",
    "            size: str | None = None, palette: str = \"hot\", \n",
    "            xlabel: str | None = None, ylabel: str | None = None) -> go.Figure:\n",
    "        \"\"\"Plot a 2D scatter plot with optimization outputs.\n",
    "        Args:\n",
    "            x_col: data (column) to be placed on the x-axis\n",
    "            y_col: data (column) to be placed on the y-axis\n",
    "            hue: data (column) associated with hue\n",
    "            size: data (column) associated with marker size\n",
    "            palette: plotly.express color palette object \n",
    "            xlabel: x-axis label text\n",
    "            ylabel: y-axis label text\n",
    "        Returns:\n",
    "            plotly.graph_objects.Figure object\n",
    "        \"\"\"\n",
    "        kwargs  =   {\n",
    "             'edgecolor' : \"k\",\n",
    "             'facecolor' : \"w\",\n",
    "             'linewidth' : 0.2,\n",
    "             'linestyle' : '-',\n",
    "            }\n",
    "        \n",
    "        fig, ax = plt.subplots(figsize=(8, 5))\n",
    "        sns.set_style('white')\n",
    "        sns.set_context(\"paper\", font_scale = 1)\n",
    "        sns.despine(right = True)\n",
    "        sns.scatterplot(\n",
    "            x = x_col, y = y_col, data = self.data, hue=hue, palette=palette, size=size,\n",
    "            marker = 'o', **kwargs, alpha = 0.95)\n",
    "        ax.legend(title='Development scenario / Firm Energy, MW', fontsize=10, \n",
    "                  title_fontsize=12, frameon=False,\n",
    "                  ncol=3)\n",
    "        plt.xlabel(xlabel, fontsize=12)\n",
    "        plt.ylabel(ylabel, fontsize=12)\n",
    "        fig.show()\n",
    "        return fig\n",
    "        \n",
    "        \n",
    "@dataclass\n",
    "class ObjectiveCalculator:\n",
    "    \"\"\" [what does this class do?]\n",
    "    \n",
    "    Attributes:\n",
    "        obj_df: Pandas datframe with objective values\n",
    "        ids: list of dam IDs from the IFC database of dams (IFC_ID)s that are \n",
    "            used in processing objectives\n",
    "        obj_names: names of objective values\n",
    "    \"\"\"\n",
    "    obj_df: pd.DataFrame\n",
    "    ids: List[NumType] = field(default_factory = list)\n",
    "    obj_names: ClassVar[List[str]] = [\n",
    "        'HP_mean', 'HP_firm', 'tot_em', 'crop_area_loss_km2', \n",
    "        'forest_area_loss_km2']\n",
    "    \n",
    "    def _filter_df(self) -> Tuple[pd.DataFrame, List[NumType]]:\n",
    "        \"\"\"Select dams from the list of ids and return a tuple with filtered dataframe and \n",
    "            missing indices. If no missing indices found, the returned list will be empty.\n",
    "        \"\"\"\n",
    "        filtered_df = self.obj_df[self.obj_df.index.isin(self.ids)]\n",
    "        found_indices = filtered_df.index.to_list()\n",
    "        missed_indices = set(self.ids) - set(found_indices)\n",
    "        return filtered_df, list(missed_indices)\n",
    "    \n",
    "    @property\n",
    "    def objectives(self) -> pd.Series:\n",
    "        \"\"\"Returns the sums of objective values for the dams with IFC ids provided in the\n",
    "            ids attribute.\n",
    "        \"\"\"\n",
    "        return self._filter_df()[0][self.obj_names].sum()\n",
    "\n",
    "    \n",
    "def map_ids(moo_ids: Sequence[int], id_map: Dict[int, int]) -> Sequence[int]:\n",
    "    \"\"\"Takes a sequence of values, e.g. dam IDs and returns a mapped sequence where\n",
    "    the original values are mapped to new values using a dictionary\"\"\"\n",
    "    return [id_map.get(item, item) for item in moo_ids]\n",
    "\n",
    "\n",
    "def find_solution_by_dam_numbers(df: pd.DataFrame, dam_ids: Set[int]) -> pd.DataFrame:\n",
    "    \"\"\"Find row(s) of dataframe matching a specified selected dams by dam id\n",
    "    Args:\n",
    "        df: Pandas dataframe with optimization results\n",
    "        dam_ids: Set of dam IDs that need to be displayed\n",
    "    Returns:\n",
    "        Pandas dataframe of optimization values for selected dam IDs\n",
    "    \"\"\"\n",
    "    return df[df['dam_ids']] == dam_ids\n",
    "\n",
    "def return_row_by_criterion(df: pd.DataFrame, criterion: str, value: Any) -> pd.Series:\n",
    "    \"\"\" \n",
    "    Args:\n",
    "        df: pandas dataframe\n",
    "        criterion: column name for which we want to find a row with the closest value.\n",
    "        value: value for which we're looking for a dam which criterion is the closest to, e.g.\n",
    "            if we have three dams with \"criterion\" = 2.4, 5.6, 6.7 and we set the value= 6.0\n",
    "            the function will return a dam for which criterion=5.6 as the difference between\n",
    "            the criterion and 6.0 is the smallest for criterion equal to 5.6\n",
    "    Returns:\n",
    "        pd.Series with the dam entry.\n",
    "    \"\"\"\n",
    "    # Calculate the absolute differences\n",
    "    differences = (df[criterion] - value).abs()\n",
    "    # Find the index of the minimum absolute difference\n",
    "    closest_index = differences.idxmin()\n",
    "    # Return the row using the index\n",
    "    closest_row = df.loc[closest_index]\n",
    "    return closest_row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "952c4059",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read file with all hydroelectric plants for optimization and their objective values\n",
    "dam_data: pd.DataFrame = pd.read_csv(pathlib.Path(\"outputs/moo/all_hp.csv\"))\n",
    "# Find IFC ids of built dams\n",
    "built_dam_ifc_ids: Set[int] = set(\n",
    "    dam_data[dam_data['status'] == 'Existing']['ifc_id'].to_list())\n",
    "# Read the ID (water model identifiers) to IFC map\n",
    "id_ifc_map: Dict[int, int] = read_id_ifc_map()\n",
    "# List IFC IDs in a sorted order of dams included in the analysis\n",
    "ifc_ids: List[int] = sorted([id_ifc_map[_id+1] for _id, _ in enumerate(id_ifc_map) ])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "732d6319",
   "metadata": {},
   "source": [
    "### Define constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "152a42f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optimization solution files with optimizations from two scenarios:\n",
    "# * nobuilt - all dams are considered for optimization as if no dams that currently exist had been built\n",
    "# * built - current state with built dams in place\n",
    "mya_nobuilt_5obj_filename = 'mya_5_obj_nobuilt.sol'\n",
    "mya_built_5obj_filename = 'mya_5_obj_built.sol'\n",
    "# Paths to output files from the algorithm with expansion / compression\n",
    "sol_file_folder = pathlib.Path('moo_solver_CPAIOR/outputs/epsilon2_5obj')\n",
    "mya_nobuilt_5obj_path_cpaior = sol_file_folder / pathlib.Path(mya_nobuilt_5obj_filename)\n",
    "mya_built_5obj_path_cpaior = sol_file_folder / pathlib.Path(mya_built_5obj_filename)\n",
    "# Load dataframe with dam ids and objective values for each dam\n",
    "dam_data_filename = pathlib.Path(\"outputs/moo/all_hp.csv\")\n",
    "# Load the mapping between ids used in the MOO algorithm and the IDs in the IFC database\n",
    "map_file_path = pathlib.Path('outputs/moo/id_to_ifc.json')\n",
    "with open(map_file_path, 'r') as file:\n",
    "    id_map = json.load(file)\n",
    "id_map: Dict[int, int] = {int(key): value for key, value in id_map.copy().items()} # Maps optim ids to ifc ids\n",
    "ifc_to_id_map: Dict[int, int] = {value: key for key, value in id_map.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b8f69ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "dam_df = pd.read_csv(dam_data_filename, index_col=0).set_index('ifc_id')\n",
    "built_dam_ids_ifc: Set[int] = set(dam_df[dam_df['status_int'] == 1].index.to_list())\n",
    "# Set of dams used for optimization\n",
    "built_dam_ids_opt: Set[int] = set_remap(built_dam_ids_ifc, ifc_to_id_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65151bbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print a number of existing dams\n",
    "num_existing_dams = dam_df[dam_df['status'] == 'Existing']['name'].count()\n",
    "print(f\"Number of existing dams: {num_existing_dams}\")\n",
    "# Print a set of existing dam ids\n",
    "existing_dam_ids = set(dam_df[dam_df['status'] == 'Existing'].index.to_list())\n",
    "print(\"Existing dam ids:\", existing_dam_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f92f7a9",
   "metadata": {},
   "source": [
    "## Call the optimizer by calling external script using subprocess (not recommended)\n",
    "### Better to read pre-calculated solutions as the optimization takes a long time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "208f2dc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# RUNS optimization as a subprocess from within this notebook\n",
    "# Switch the rerun_CPAIOR flag to True if you want to run the optimization yourself\n",
    "# (may take 1hr+ to execute). Otherwise keep it as False and read pre-saved optimization\n",
    "# results in .sol files.\n",
    "rerun_CPAIOR: bool = False\n",
    "\n",
    "# Use subprocess to call and execute the shell script\n",
    "if rerun_CPAIOR:\n",
    "    # Specify the path to your shell script\n",
    "    script_path = 'moo_solver_CPAIOR/run_myanmar_dam_selection.sh'\n",
    "    try:\n",
    "        subprocess.run(['bash', script_path], check=True)\n",
    "        print(\"Optimization runs successful.\")\n",
    "    except subprocess.CalledProcessError as e:\n",
    "        print(f\"Error executing script: {e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43dbdd2b",
   "metadata": {},
   "source": [
    "## Parse outputs - 5 objective optimization\n",
    "<h2> NOTE: The code below might take a long time to execut. If you want to plot the figure using saved data, move down to section \n",
    "<font color=\"red\">\n",
    "<b>PLOT: Plot from saved data</b>\n",
    "    </font>  \n",
    "    and set \n",
    "    <pre><code data-trim>\n",
    "<div id=\"awesomecpp\"></div>\n",
    "load_from_file: bool = True\n",
    "</code></pre>\n",
    "</h2>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95415ea8",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_rows_to_remove: int = 3  # Remove n rows per n+1 rows, e.g. if value is 3 every 3 out of four rows \n",
    "                            # will be removed\n",
    "filter_dataframe: bool = True\n",
    "# Save to json and/or csv if boolean flags for each are set to True\n",
    "save_to_json, save_to_csv = False, False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9280ca96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parse solutions from optimizations with non-built dam scenarios\n",
    "# Parse solutions with no built dams\n",
    "mya_nobuilt_parser_cpaior_5obj = SolutionFileParser(mya_nobuilt_5obj_path_cpaior)\n",
    "mya_nobuilt_parser_cpaior_5obj.parse()\n",
    "if save_to_json:\n",
    "    mya_nobuilt_parser_cpaior_5obj.to_json(sol_file_folder / pathlib.Path('mya_5_obj_nobuilt.json'))\n",
    "if save_to_csv:\n",
    "    mya_nobuilt_parser_cpaior_5obj.to_csv(sol_file_folder / pathlib.Path('mya_5_obj_nobuilt.csv'))\n",
    "# Concetenate dataframes into `merged_df`\n",
    "df_nobuilt = reduce_mem_usage(mya_nobuilt_parser_cpaior_5obj.solutions_df)\n",
    "# Remove the unwanted objects from memory\n",
    "del mya_nobuilt_parser_cpaior_5obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf9c869c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parse solutions with built dams\n",
    "mya_built_parser_cpaior_5obj = SolutionFileParser(mya_built_5obj_path_cpaior)\n",
    "mya_built_parser_cpaior_5obj.parse()\n",
    "if save_to_json:\n",
    "    mya_built_parser_cpaior_5obj.to_json(sol_file_folder / pathlib.Path('mya_5_obj_built.json'))\n",
    "if save_to_csv:\n",
    "    mya_built_parser_cpaior_5obj.to_csv(sol_file_folder / pathlib.Path('mya_5_obj_built.csv'))\n",
    "df_built = reduce_mem_usage(mya_built_parser_cpaior_5obj.solutions_df)\n",
    "# Remove the unwanted objects from memory\n",
    "del mya_built_parser_cpaior_5obj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d7570b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter the dataframes before merging in the next step\n",
    "df_built = df_built.sort_values(by=['energy'], ascending=True)\n",
    "df_nobuilt = df_nobuilt.sort_values(by=['energy'], ascending=True)\n",
    "# Filter the dataframe if filtering is selected\n",
    "if filter_dataframe:\n",
    "    df_built_filt = get_every_n_row(df_built, no_rows_to_remove+1)\n",
    "    df_nobuilt_filt = get_every_n_row(df_nobuilt, no_rows_to_remove+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c90a399",
   "metadata": {},
   "source": [
    "### Combine dataframes with results with 'built' and 'notbuilt' scenarios, sorf values, add new columns and reduce memory usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "050430d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nobuilt_filt['Scenario'] = \"Not Built\"\n",
    "df_built_filt['Scenario'] = \"Built\"\n",
    "# Merge the built and nobuilt dataframes. They are concatenated by the row dimension in\n",
    "# order \"Not Built\" -> \"Built\"\n",
    "merged_df = pd.concat([df_built_filt, df_nobuilt_filt], ignore_index=True)\n",
    "old_new_col_map = {\n",
    "    'energy': \"Mean annual HP, [MW]\",\n",
    "    'ghg': 'GHG emissions [tonne CO<sub>2,eq</sub>/year]', \n",
    "    'firm_energy': 'Firm HP, [MW]',\n",
    "    'loss_agri': 'Agricultural land loss, [km<sup>2</sup>]',\n",
    "    'loss_forest': 'Deforestation, [km<sup>2</sup>]',\n",
    "    'num_dams': 'No. of selected dams',\n",
    "    'dam_ids': 'Dam IDs',\n",
    "    'land_loss': 'Land loss, [km<sup>2</sup>]',\n",
    "    'ghg_intensity': 'GHG intensity [gCO<sub>2,eq</sub>/kWh]'}\n",
    "merged_df.rename(columns=old_new_col_map, inplace=True)\n",
    "merged_df['Dam IDs'] = merged_df['Dam IDs'].apply(set)\n",
    "merged_df['Firm Power Ratio, [%]'] = \\\n",
    "    merged_df['Firm HP, [MW]'] / merged_df['Mean annual HP, [MW]'] * 100\n",
    "merged_df['Scenario, [1/0]'] = merged_df['Scenario'].map({'Built': 1, 'Not Built': 0})\n",
    "\n",
    "# Define bin edges for land loss\n",
    "bins = [0, 300, 600, 1000, 1500, 2000]\n",
    "# Define labels for the bins\n",
    "labels = ['0-300 km2', '300-500 km2', '500-1000 km2', '1000-1500 km2', '1500-2000 km2']\n",
    "merged_df[\"Loss of Land [km<sup>2</sup>]\"] = pd.cut(\n",
    "    merged_df['Land loss, [km<sup>2</sup>]'], bins=bins, labels=labels, right=False)\n",
    "# Arrange by status and energy in ascending order\n",
    "merged_df = merged_df.sort_values(by=['Scenario', 'Mean annual HP, [MW]'], ascending=True)\n",
    "# Introduce new columns\n",
    "merged_df['HP Production [GWh/year]'] = merged_df[\"Mean annual HP, [MW]\"] * 365.25 * 24 / 1_000\n",
    "merged_df['Mean HP [GWh/d]'] = merged_df[\"Mean annual HP, [MW]\"] * 24 / 1_000\n",
    "merged_df['Firm HP [GWh/d]'] = merged_df['Firm HP, [MW]'] * 24 / 1_000\n",
    "\n",
    "# Reduce size of some data in merged_df\n",
    "merged_df['Scenario, [1/0]'] = merged_df['Scenario, [1/0]'].astype('uint8')\n",
    "merged_df[\"Loss of Land [km<sup>2</sup>]\"] = merged_df[\"Loss of Land [km<sup>2</sup>]\"].astype('category')\n",
    "merged_df['No. of selected dams'] = merged_df['No. of selected dams'].astype('int8')\n",
    "# Use an automated method\n",
    "merged_df = reduce_mem_usage(merged_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46f6cd54",
   "metadata": {},
   "source": [
    "### Perform nondominated sorting in 2D in order to visualise the lower boundary of the Pareto front"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e7135cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform non-dominated sorting in 2D for two 5D nondominated fronts: for Built and NotBuilt scenarios\n",
    "# Create a pareto dominant front for not built data\n",
    "col_names = [\n",
    "    \"Mean annual HP, [MW]\", 'GHG emissions [tonne CO<sub>2,eq</sub>/year]', \n",
    "    'GHG intensity [gCO<sub>2,eq</sub>/kWh]', 'Firm Power Ratio, [%]', \n",
    "    'Land loss, [km<sup>2</sup>]', 'Dam IDs']\n",
    "xy_pairs_built = merged_df\\\n",
    "    .loc[merged_df['Scenario'] =='Built', col_names]\n",
    "xy_pairs_built[\"Mean annual HP, [MW]\"] = xy_pairs_built[\"Mean annual HP, [MW]\"] * -1\n",
    "xy_pairs_built_list = xy_pairs_built.to_numpy().tolist()\n",
    "\n",
    "xy_pairs_nobuilt = merged_df\\\n",
    "    .loc[merged_df['Scenario'] =='Not Built', col_names]\n",
    "xy_pairs_nobuilt[\"Mean annual HP, [MW]\"] = xy_pairs_nobuilt[\"Mean annual HP, [MW]\"] * -1\n",
    "xy_pairs_nobuilt_list = xy_pairs_nobuilt.to_numpy().tolist()\n",
    "\n",
    "xy_pair_array_built_np = np.array(xy_pairs_built_list)\n",
    "xy_pair_array_nobuilt_np = np.array(xy_pairs_nobuilt_list)\n",
    "\n",
    "# Find non-dominated fronts, i.e. indexes of nondominated points. List them in the order of\n",
    "# increasing HP production\n",
    "non_dom_front_built = pg.non_dominated_front_2d(points=xy_pair_array_built_np[:,:2])[::-1]\n",
    "non_dom_front_nobuilt = pg.non_dominated_front_2d(points=xy_pair_array_nobuilt_np[:,:2])[::-1]\n",
    "# Convert back from negative to positive values\n",
    "xy_pair_array_built_np[:,0] = xy_pair_array_built_np[:,0] * -1\n",
    "xy_pair_array_nobuilt_np[:,0] = xy_pair_array_nobuilt_np[:,0] * -1\n",
    "# Select nondominated points\n",
    "xy_nondom_built_np = xy_pair_array_built_np[non_dom_front_built]\n",
    "xy_nondom_nobuilt_np = xy_pair_array_nobuilt_np[non_dom_front_nobuilt]\n",
    "xy_nondom_all = np.concatenate((xy_nondom_built_np, xy_nondom_nobuilt_np), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dc5eb3a",
   "metadata": {},
   "source": [
    "### Convert nondominated solutions back to dataframes and add new columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c31c637",
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_nondom_built_df = pd.DataFrame(xy_nondom_built_np, columns=col_names)\n",
    "xy_nondom_built_df['Scenario'] = \"Built\"\n",
    "xy_nondom_nobuilt_df = pd.DataFrame(xy_nondom_nobuilt_np, columns=col_names)\n",
    "xy_nondom_nobuilt_df['Scenario'] = \"Not Built\"\n",
    "xy_nondom_df = pd.concat([xy_nondom_built_df, xy_nondom_nobuilt_df], ignore_index=True)\n",
    "xy_nondom_df['HP Production [GWh/year]'] = xy_nondom_df[\"Mean annual HP, [MW]\"] * 365.25 * 24 / 1_000\n",
    "xy_nondom_df[\"Loss of Land [km<sup>2</sup>]\"] = pd.cut(\n",
    "    xy_nondom_df['Land loss, [km<sup>2</sup>]'], bins=bins, labels=labels, right=False)\n",
    "em_int_nondom_df = xy_nondom_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c1cfd55",
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_nondom_df.head(3)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "cf678af0",
   "metadata": {},
   "source": [
    "# Find the solution ids of the 2D nondominated song in the order of increasing annual energy production\n",
    "# Plot the results (for built and not built scenarios) as a scatterplot\n",
    "em_int_nondom_built = merged_df.loc[merged_df['Scenario'] =='Built'].iloc[non_dom_front_built]\n",
    "em_int_nondom_nobuilt = merged_df.loc[merged_df['Scenario'] =='Not Built'].iloc[non_dom_front_nobuilt]\n",
    "em_int_nondom_df = pd.concat(\n",
    "    [em_int_nondom_built, em_int_nondom_nobuilt], ignore_index=True).sort_values(\"Mean annual HP, [MW]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "708a6e25",
   "metadata": {},
   "outputs": [],
   "source": [
    "save_to_file: bool = False\n",
    "if save_to_file:\n",
    "    em_int_nondom_df.to_csv('em_int_nondom_df.csv', index=False)\n",
    "    merged_df.to_csv('merged_df.csv', index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05cd7809",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Statistics\n",
    "number_of_solutions = len(df_nobuilt) + len(df_built)\n",
    "print(f\"Total number of solutions : {number_of_solutions}\")\n",
    "print(f\"Scenario with built constructed dams {len(df_built)} solutions\")\n",
    "print(f\"Scenario with zero constructed dams {len(df_nobuilt)} solutions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6141b68c",
   "metadata": {},
   "source": [
    "<h2>\n",
    "<font color=\"red\">\n",
    "<b>PLOT: Plot from saved data</b>\n",
    "    </font>    \n",
    "</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f40a51ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "load_from_file: bool = False\n",
    "if load_from_file:\n",
    "    merged_df = pd.read_csv('merged_df.csv')\n",
    "    em_int_nondom_df = pd.read_csv('em_int_nondom_df.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7df4acec",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90cd1228",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find a solution with built dams with minumum energy generation (solution with built dams only)\n",
    "built_min = merged_df[merged_df['Scenario']=='Built'].iloc[0]\n",
    "built_min_damids = built_min['Dam IDs']\n",
    "built_min_energy = built_min[\"Mean annual HP, [MW]\"]\n",
    "built_min_ghg_intensity = built_min['GHG intensity [gCO<sub>2,eq</sub>/kWh]']\n",
    "built_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "895535f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find total energy from input data using selected dams as an input\n",
    "oc_built = ObjectiveCalculator(dam_df, ids=set_remap(built_min['Dam IDs'], id_map))\n",
    "print(\n",
    "    f\"Min HP from optimization: {built_min_energy} MW, min HP from input data: {oc_built.objectives['HP_mean']} MW\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0dcdddf",
   "metadata": {},
   "outputs": [],
   "source": [
    "#em_int_nondom_df[em_int_nondom_df['Scenario']==\"Built\"].sort_values(by='No. of selected dams').head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f9362ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check retrieval of solutions from the dataframe (not used)\n",
    "dam_ids_7848 = merged_df.loc[7848, 'Dam IDs']\n",
    "dam_ids_7599 = merged_df.loc[7599, 'Dam IDs']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21f03929",
   "metadata": {},
   "outputs": [],
   "source": [
    "em_int_nondom_df.head()"
   ]
  },
  {
   "cell_type": "raw",
   "id": "90f75491",
   "metadata": {},
   "source": [
    "plot_data = pd.to_numeric(em_int_nondom_df, errors='ignore')\n",
    "plot_data.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9482800",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get selected data points\n",
    "# Separate DAM IDs and No. of Dams from em_int_nondom_df and later join with plot data\n",
    "#plot_data = em_int_nondom_df.drop(columns=['Dam IDs'])\n",
    "#dam_ids = em_int_nondom_df[['Scenario','Dam IDs']]\n",
    "em_int_nondom_df['HP Production [TWh/year]'] = em_int_nondom_df['HP Production [GWh/year]'] / 1_000\n",
    "plot_data = em_int_nondom_df.apply(pd.to_numeric, errors='ignore')\n",
    "#\n",
    "\n",
    "current_solution = (built_min['HP Production [GWh/year]']/1_000, built_min_ghg_intensity)\n",
    "notbuilt_current = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario']==\"Not Built\"], \n",
    "    'HP Production [TWh/year]', value=built_min['HP Production [GWh/year]']/1_000)\n",
    "\n",
    "# Get selected data points\n",
    "current_solution = (built_min['HP Production [GWh/year]']/1_000, built_min_ghg_intensity)\n",
    "notbuilt_current = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario'] == \"Not Built\"], 'HP Production [TWh/year]', \n",
    "    value=built_min['HP Production [GWh/year]']/1_000)\n",
    "# Find points (soluions) corresponding to HP production of 100 and 200 TWh/year respectively.\n",
    "built_100 = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario'] == \"Built\"], 'HP Production [TWh/year]', value=100)\n",
    "notbuilt_100 = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario'] == \"Not Built\"], 'HP Production [TWh/year]', value=100)\n",
    "built_200 = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario'] == \"Built\"], 'HP Production [TWh/year]', value=200)\n",
    "notbuilt_200 = return_row_by_criterion(\n",
    "    plot_data[plot_data['Scenario'] == \"Not Built\"], 'HP Production [TWh/year]', value=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba07792e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_data.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ef1b44a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_data.head(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08090fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(notbuilt_current['Dam IDs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "732caba9",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(built_min['Dam IDs'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b98dd953",
   "metadata": {},
   "outputs": [],
   "source": [
    "notbuilt_current"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbb32a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define bin edges for firm power ratio\n",
    "firm_power_bins2 = [0, 25, 45, 60, 75, 95]\n",
    "firm_power_labels2 = ['0%-25%', '25%-45%', '45%-60%', '60%-75%', '75%-95%']\n",
    "em_int_nondom_df['Firm Power Ratio Cat'] = pd.cut(\n",
    "    em_int_nondom_df['Firm Power Ratio, [%]'], bins=firm_power_bins2, \n",
    "    labels=firm_power_labels2, right=False)\n",
    "merged_df['Firm Power Ratio Cat'] = pd.cut(\n",
    "    merged_df['Firm Power Ratio, [%]'], bins=firm_power_bins2, \n",
    "    labels=firm_power_labels2, right=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "758715f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_x_y_intensity_plot(data: pd.Series) -> Tuple[float, float]:\n",
    "    \"\"\"Return a tuple consisting of:\n",
    "        (a) annual HP generation and (b) GHG intensity\"\"\"\n",
    "    return (data['HP Production [TWh/year]'], data['GHG intensity [gCO<sub>2,eq</sub>/kWh]'])\n",
    "\n",
    "def get_x_y_emissions_plot(data: pd.Series) -> Tuple[float, float]:\n",
    "    \"\"\"Return a tuple consisting of:\n",
    "        (a) annual HP generation and (b) total GHG emissions\"\"\"\n",
    "    return (\n",
    "        data['HP Production [TWh/year]'], \n",
    "        data['GHG emissions [tonne CO<sub>2,eq</sub>/year]'] / 1_000_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6b59bd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_ghg_intensity_figure(\n",
    "        ax: matplotlib.axes.Axes,\n",
    "        annotation_font_size: int = 17, \n",
    "        cur_sol_marker_size: int = 20,\n",
    "        tick_label_size: int = 15,\n",
    "        label_font_size: int = 15,\n",
    "        sol_marker_size: int = 14) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    marker_edge_width = 0.7\n",
    "    marker_alpha = 0.9\n",
    "    bbox = dict(boxstyle=\"round\", pad=0.15, facecolor='none', edgecolor='none')\n",
    "    annotation_shrink = 0.03\n",
    "\n",
    "    kwargs  =   {\n",
    "        'edgecolor':'grey',\n",
    "        'marker': 'o',\n",
    "        'facecolor':'none',\n",
    "        'linewidth':0.05,\n",
    "        'linestyle':'-',\n",
    "        'alpha': 0.05\n",
    "    }\n",
    "    \n",
    "    kwargs_ghg_intensity  =   {\n",
    "        'edgecolor':'k',\n",
    "        #'marker': 'o',\n",
    "        'alpha': 0.05,\n",
    "        'linewidth':0.05,\n",
    "        'linestyle':'-',\n",
    "    }\n",
    "    custom_palette = [\"#FFFFFF\", \"#808080\", \"#000000\", \"#3A0CA3\", \"#4361EE\", \"#4CC9F0\"]\n",
    "    merged_df['HP Production [TWh/year]'] = merged_df['HP Production [GWh/year]'] / 1_000\n",
    "    \n",
    "    sns.scatterplot(\n",
    "        x = 'HP Production [TWh/year]', \n",
    "        y = 'GHG intensity [gCO<sub>2,eq</sub>/kWh]', \n",
    "        ax = ax,\n",
    "        data = merged_df, \n",
    "        #palette='Set2', #'YlOrRd', \n",
    "        #hue='Firm Power Ratio Cat',\n",
    "        s = 50,\n",
    "        **kwargs)\n",
    "    \n",
    "    sns.set_style('white')\n",
    "    sns.set_context(\"paper\", font_scale = 1)\n",
    "    sns.despine(right = True)\n",
    "    sns.scatterplot(\n",
    "        x = 'HP Production [TWh/year]', y = 'GHG intensity [gCO<sub>2,eq</sub>/kWh]', \n",
    "        data = em_int_nondom_df, palette='Set2',\n",
    "        ax =ax,\n",
    "        size='Firm Power Ratio Cat', hue=\"Loss of Land [km<sup>2</sup>]\",\n",
    "        sizes = {'0%-25%': 10, '25%-45%': 50, '45%-60%': 100, '60%-75%': 150, '75%-95%': 200},\n",
    "        #size=\"Loss of Land [km<sup>2</sup>]\", hue='Firm Power Ratio Cat', \n",
    "        #sizes = {'0-400 km2': 10, '400-800 km2': 30, '800-1200 km2': 60, '1200-1600 km2': 120},\n",
    "        **kwargs_ghg_intensity)\n",
    "    legend1 = ax.legend(loc=\"upper right\", fontsize='large', ncols=2)\n",
    "    legend1.set_frame_on(False)\n",
    "    updated_text = [\n",
    "        \"Loss of Land\",\n",
    "        \"$0-300$ km$^2$\",\n",
    "        \"$300-600$ km$^2$\",\n",
    "        \"$600-1000$ km$^2$\",\n",
    "        \"$1000-1500$ km$^2$\",\n",
    "        \"$1500-2000$ km$^2$\",\n",
    "        \"Firm Power Ratio\",\n",
    "        \"$0\\% - 25\\%$\",\n",
    "        \"$25\\% - 45\\%$\",\n",
    "        \"$45\\% - 60\\%$\",\n",
    "        \"$60\\% - 75\\%$\",\n",
    "        \"$75\\% - 95\\%$\"\n",
    "    ]\n",
    "    for ix, text in enumerate(updated_text):\n",
    "        legend1.get_texts()[ix].set_text(text)\n",
    "    ax.axvline(\n",
    "        x=current_solution[0], color='grey', linestyle='--', linewidth=2,\n",
    "        label='Current mean annual power production, GWh/year')\n",
    "    ax.text(current_solution[0]-8, 34, 'Current mean annual HP production', \n",
    "             rotation=90, va='center', ha='center', fontsize=12)\n",
    "\n",
    "    # Plot chosen solution scenarios\n",
    "    ax.plot(\n",
    "        current_solution[0], current_solution[1], \n",
    "        marker='*', markersize=cur_sol_marker_size, color='yellow', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate(\n",
    "        '$I_{b}$', \n",
    "        xy=(current_solution[0], current_solution[1]), \n",
    "        xytext=(current_solution[0] + 20, current_solution[1] - 7.5),\n",
    "        arrowprops=dict(\n",
    "         facecolor='black', \n",
    "         shrink=annotation_shrink, \n",
    "         width=2, headwidth = 8 ),\n",
    "        fontsize=annotation_font_size, bbox=bbox)\n",
    "    # Current not built\n",
    "    curr_nbuilt_x, curr_nbuilt_y = get_x_y_intensity_plot(notbuilt_current)\n",
    "    ax.annotate(\n",
    "        '$I_{nb}$',\n",
    "        xy=(curr_nbuilt_x, curr_nbuilt_y), \n",
    "        xytext=(curr_nbuilt_x + 13, curr_nbuilt_y + 10),\n",
    "        arrowprops=dict(\n",
    "         facecolor='black', \n",
    "         shrink=annotation_shrink, \n",
    "         width=2, \n",
    "         headwidth = 8 ),\n",
    "        fontsize=annotation_font_size, bbox=bbox)\n",
    "    ax.plot(\n",
    "        curr_nbuilt_x, curr_nbuilt_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    # Other values\n",
    "    built_100_x, built_100_y = get_x_y_intensity_plot(built_100)\n",
    "    ax.plot(\n",
    "        built_100_x, built_100_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$II_b$', xy=(built_100_x, built_100_y), xytext=(built_100_x + 20, built_100_y + 13),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    notbuilt_100_x, notbuilt_100_y = get_x_y_intensity_plot(notbuilt_100)\n",
    "    ax.plot(\n",
    "        notbuilt_100_x, notbuilt_100_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$II_{nb}$', xy=(notbuilt_100_x, notbuilt_100_y), xytext=(notbuilt_100_x + 38, notbuilt_100_y + 19),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    built_200_x, built_200_y = get_x_y_intensity_plot(built_200)\n",
    "    ax.plot(\n",
    "        built_200_x, built_200_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$III_{b}$', xy=(built_200_x, built_200_y), xytext=(built_200_x + 15, built_200_y + 9),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    notbuilt_200_x, notbuilt_200_y = get_x_y_intensity_plot(notbuilt_200)\n",
    "    ax.plot(\n",
    "        notbuilt_200_x, notbuilt_200_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$III_{nb}$', xy=(notbuilt_200_x, notbuilt_200_y), \n",
    "                 xytext=(notbuilt_200_x + 25, notbuilt_200_y -7),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "\n",
    "    ax.tick_params(axis='x', labelsize=tick_label_size)\n",
    "    ax.tick_params(axis='y', labelsize=tick_label_size)\n",
    "    ax.set_ylim(-5,70)\n",
    "    ax.set_xlim(0,250)\n",
    "    ax.set_xlabel(\"HP Production, TWh/year\", fontsize=label_font_size)\n",
    "    ax.set_ylabel(\"Biogenic emission intensity, gCO$_{2,eq}$/kWh\", fontsize=label_font_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc13af8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df.sort_values(by=['Mean annual HP, [MW]'], ascending=True).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c897558",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(merged_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c003207d",
   "metadata": {},
   "outputs": [],
   "source": [
    "em_int_nondom_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f32f26f",
   "metadata": {},
   "outputs": [],
   "source": [
    "built_max_all = merged_df[merged_df['Scenario']=='Built'].iloc[-1]\n",
    "built_max_nondom = em_int_nondom_df[em_int_nondom_df['Scenario']=='Built'].iloc[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8694a9df",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_ghg_emissions_figure(\n",
    "        ax: matplotlib.axes.Axes, \n",
    "        annotation_font_size: int = 18, \n",
    "        tick_label_size: int = 15,\n",
    "        sol_marker_size: int = 14,\n",
    "        cur_sol_marker_size: int = 20,\n",
    "        label_font_size: int = 15,\n",
    "        vline_font_size: int = 12) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    marker_edge_width = 0.7\n",
    "    marker_alpha = 1\n",
    "    bbox = dict(boxstyle=\"round\", pad=0.15, facecolor='none', edgecolor='none')\n",
    "    annotation_shrink = 0.03\n",
    "    \n",
    "    merged_df['HP Production [TWh/year]'] = merged_df['HP Production [GWh/year]'] / 1_000\n",
    "    merged_df['GHG emissions [Mt CO2/year]'] = \\\n",
    "        merged_df['GHG emissions [tonne CO<sub>2,eq</sub>/year]'] / 1_000_000\n",
    "    em_int_nondom_df['GHG emissions [Mt CO2/year]'] = \\\n",
    "        em_int_nondom_df['GHG emissions [tonne CO<sub>2,eq</sub>/year]'] / 1_000_000\n",
    "\n",
    "    kwargs  =   {\n",
    "        'edgecolor':'grey',\n",
    "        'marker': 'o',\n",
    "        'facecolor':'none',\n",
    "        'linewidth':0.05,\n",
    "        'linestyle':'-',\n",
    "        'alpha': 0.05\n",
    "    }\n",
    "    kwargs_nondom  =   {\n",
    "        'edgecolor':'k',\n",
    "        'marker': 'o',\n",
    "        'linewidth':0.05,\n",
    "        'linestyle':'-',\n",
    "        'alpha': 0.05\n",
    "    }\n",
    "    custom_palette = [\"#FFFFFF\", \"#808080\", \"#000000\", \"#3A0CA3\", \"#4361EE\", \"#4CC9F0\"]\n",
    "\n",
    "    sns.set_style('white')\n",
    "    sns.set_context(\"paper\", font_scale = 1)\n",
    "    sns.despine(right = True)\n",
    "    # Plot all results\n",
    "    \n",
    "    sns.scatterplot(\n",
    "        x = 'HP Production [TWh/year]', \n",
    "        y = 'GHG emissions [Mt CO2/year]', \n",
    "        ax = ax,\n",
    "        data = merged_df, \n",
    "        #palette='Set2', #'YlOrRd', \n",
    "        #hue='Firm Power Ratio Cat',\n",
    "        s = 50,\n",
    "        **kwargs)\n",
    "\n",
    "    #ax.legend(title='Development scenario / Firm Energy, MW', fontsize=10, \n",
    "    #          title_fontsize=12, frameon=False,\n",
    "    #          ncol=3)\n",
    "\n",
    "    # Plot nondominated results\n",
    "    sns.scatterplot(\n",
    "        x = 'HP Production [TWh/year]', y = 'GHG emissions [Mt CO2/year]', \n",
    "        data = em_int_nondom_df, \n",
    "        palette='Set2',\n",
    "        size='Firm Power Ratio Cat', hue=\"Loss of Land [km<sup>2</sup>]\", \n",
    "        ax = ax,\n",
    "        #sizes = {'0-400 km2': 10, '400-800 km2': 30, '800-1200 km2': 60, '1200-1600 km2': 120},\n",
    "        sizes = {'0%-25%': 10, '25%-45%': 50, '45%-60%': 100, '60%-75%': 150, '75%-95%': 200},\n",
    "        #sizes = (1,200),\n",
    "        **kwargs_nondom)\n",
    "    \n",
    "\n",
    "    legend2 = ax.legend(loc=\"upper left\", fontsize='large', bbox_to_anchor=(0.08, 1.05), ncol=2)\n",
    "    legend2.set_frame_on(False)\n",
    "    updated_text = [\n",
    "        \"Loss of Land\",\n",
    "        \"$0-300$ km$^2$\",\n",
    "        \"$300-600$ km$^2$\",\n",
    "        \"$600-1000$ km$^2$\",\n",
    "        \"$1000-1500$ km$^2$\",\n",
    "        \"$1500-2000$ km$^2$\",\n",
    "        \"Firm Power Ratio\",\n",
    "        \"$0\\% - 25\\%$\",\n",
    "        \"$25\\% - 45\\%$\",\n",
    "        \"$45\\% - 60\\%$\",\n",
    "        \"$60\\% - 75\\%$\",\n",
    "        \"$75\\% - 95\\%$\"\n",
    "    ]\n",
    "    for ix, text in enumerate(updated_text):\n",
    "        legend2.get_texts()[ix].set_text(text)\n",
    "\n",
    "    # Plot solution points\n",
    "    current_solution = (\n",
    "        built_min['HP Production [GWh/year]']/1_000, \n",
    "        built_min['GHG emissions [tonne CO<sub>2,eq</sub>/year]']/1_000/1_000)\n",
    "    ax.plot(\n",
    "        current_solution[0], current_solution[1], \n",
    "        marker='*', markersize=cur_sol_marker_size, color='yellow', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$I_{b}$', xy=(current_solution[0], current_solution[1]), \n",
    "                 xytext=(current_solution[0] + 10, current_solution[1] + 0.55),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    \n",
    "    # Current not built (point)\n",
    "    curr_nbuilt_x, curr_nbuilt_y = get_x_y_emissions_plot(notbuilt_current)\n",
    "    ax.plot(\n",
    "        curr_nbuilt_x, curr_nbuilt_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$I_{nb}$', xy=(curr_nbuilt_x, curr_nbuilt_y), xytext=(curr_nbuilt_x + 25, curr_nbuilt_y + 1.5),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    # Other point values\n",
    "    built_100_x, built_100_y = get_x_y_emissions_plot(built_100)\n",
    "    ax.plot(\n",
    "        built_100_x, built_100_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$II_b$', xy=(built_100_x, built_100_y), xytext=(built_100_x + 11, built_100_y + 0.64),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    notbuilt_100_x, notbuilt_100_y = get_x_y_emissions_plot(notbuilt_100)\n",
    "    ax.plot(\n",
    "        notbuilt_100_x, notbuilt_100_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$II_{nb}$', xy=(notbuilt_100_x, notbuilt_100_y), xytext=(notbuilt_100_x + 30, notbuilt_100_y + 1.6),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    built_200_x, built_200_y = get_x_y_emissions_plot(built_200)\n",
    "    ax.plot(\n",
    "        built_200_x, built_200_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$III_{b}$', xy=(built_200_x, built_200_y), xytext=(built_200_x - 15, built_200_y + 0.7),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox)\n",
    "    notbuilt_200_x, notbuilt_200_y = get_x_y_emissions_plot(notbuilt_200)\n",
    "    ax.plot(\n",
    "        notbuilt_200_x, notbuilt_200_y, \n",
    "        marker='o', markersize=sol_marker_size, color='white', markeredgecolor='k', \n",
    "        markeredgewidth=marker_edge_width, alpha=marker_alpha)\n",
    "    ax.annotate('$III_{nb}$', xy=(notbuilt_200_x, notbuilt_200_y), \n",
    "                 xytext=(notbuilt_200_x + 26, notbuilt_200_y -0.25),\n",
    "                 arrowprops=dict(facecolor='black', shrink=annotation_shrink, width=2, headwidth = 8 ),\n",
    "                 fontsize=annotation_font_size, bbox=bbox) \n",
    "    \n",
    "    # Add vertical line denoting the current HP production\n",
    "    ax.axvline(\n",
    "        x=current_solution[0], color='grey', linestyle='--', linewidth=2,\n",
    "        label='Current mean annual power production, GWh/year')\n",
    "    ax.text(current_solution[0]-8, 2.5, 'Current mean annual HP production', \n",
    "             rotation=90, va='center', ha='center', fontsize=vline_font_size)\n",
    "\n",
    "    arrow_x = current_solution[0] + 20\n",
    "    arrow_y_start = 0  # y-coordinate for the starting point of the arrow\n",
    "    arrow_y_end = 0.95  # y-coordinate for the ending point of the arrow\n",
    "    arrow_text = 'Forgone opportunity'  # Text to be displayed next to the arrow\n",
    "    # Plot the line with arrows\n",
    "    ax.annotate('', xy=(arrow_x, arrow_y_end), xytext=(arrow_x, arrow_y_start),\n",
    "                 arrowprops=dict(arrowstyle='<->', color='k', lw=1.5), annotation_clip=False)\n",
    "    # Add text annotation next to the arrow\n",
    "    ax.text(arrow_x + 3, (arrow_y_start + arrow_y_end) / 2 - 0.05, arrow_text, color='k', fontsize=12)\n",
    "\n",
    "    ax.tick_params(axis='x', labelsize=tick_label_size)\n",
    "    ax.tick_params(axis='y', labelsize=tick_label_size)\n",
    "    ax.set_ylim(-0.2,5)\n",
    "    ax.set_xlim(0,250)\n",
    "    ax.set_xlabel(\"HP Production, TWh/year\", fontsize=label_font_size)\n",
    "    ax.set_ylabel(\"GHG Emissions, Mt CO$_2$ / annum\", fontsize=label_font_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad7ff65b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot unit emissions and total emission figures (Pareto fronts)\n",
    "fig, (ax1, ax2) = plt.subplots(2,1,figsize=(8, 10))\n",
    "plt.subplots_adjust(hspace=0.3)\n",
    "make_ghg_intensity_figure(ax1)\n",
    "make_ghg_emissions_figure(ax2)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "9d5ed56e",
   "metadata": {},
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "make_ghg_emissions_figure(\n",
    "        ax, \n",
    "        annotation_font_size = 18, \n",
    "        tick_label_size = 17,\n",
    "        sol_marker_size = 15,\n",
    "        cur_sol_marker_size = 22,\n",
    "        vline_font_size = 15,\n",
    "        label_font_size = 16)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "75177e3b",
   "metadata": {},
   "source": [
    "# Can take lots of resources so currently commented out\n",
    "fig.savefig(pathlib.Path('figures/moo/GHG_intensity_emissions_HP_vis_abstract_smaller.svg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0938bc28",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(pathlib.Path('figures/moo/GHG_intensity_emissions_HP_vis_abstract.pdf'), format=\"pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32f1d201",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig.savefig(pathlib.Path('figures/moo/GHG_intensity_emissions_HP_plot.png'), transparent=True, dpi=600)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87ad29cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a dictionary with scenarios as keys and sets of constructed dams as values\n",
    "def map_ids_to_ifc(scenario_ids: Dict[str, Set[int]], id_map = id_ifc_map) -> Dict[str, List[int]]:\n",
    "    \"\"\" \"\"\"\n",
    "    return {key : list(set_remap(optim_ids, id_ifc_map)) for key, optim_ids in scenario_ids.items()}\n",
    "\n",
    "sc_dams: Dict[str, Set[int]] = map_ids_to_ifc({\n",
    "    \"Ib\": built_min['Dam IDs'],\n",
    "    \"Inb\": notbuilt_current['Dam IDs'],\n",
    "    \"IIb\": built_100['Dam IDs'],\n",
    "    \"IInb\": notbuilt_100['Dam IDs'],\n",
    "    \"IIIb\": built_200['Dam IDs'],\n",
    "    \"IIInb\": notbuilt_200['Dam IDs']\n",
    "})\n",
    "with open(pathlib.Path('intermediate/optim_scenarios.json'), 'w') as file:  \n",
    "    json_string = json.dumps(sc_dams, indent=4)\n",
    "    file.write(json_string)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "5e765e28",
   "metadata": {},
   "source": [
    "# TODO: Add labels\n",
    "labels (dict with str keys and str values (default {})) – By default, column names are used in the figure for axis titles, legend entries and hovers. This parameter allows this to be overridden. The keys of this dict should correspond to column names, and the values should correspond to the desired label to be displayed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c780ff20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This only works if previous figures have not been created\n",
    "merged_df['GHG emissions<br>[Mt CO2,eq/year]'] = \\\n",
    "    merged_df['GHG emissions [tonne CO<sub>2,eq</sub>/year]']/1000_000\n",
    "merged_df.rename(columns={\n",
    "    'Agricultural land loss, [km<sup>2</sup>]' : 'Agricultural land<br>loss, [km2]',\n",
    "    'Deforestation, [km<sup>2</sup>]' : 'Deforestation<br>[km2]',\n",
    "    'Land loss, [km<sup>2</sup>]' : 'Land loss<br>[km2]',\n",
    "    'GHG intensity [gCO<sub>2,eq</sub>/kWh]' : 'Biogenic GHG intensity<br>[gCO2,eq/kWh]',\n",
    "    'Mean HP [GWh/d]' : 'Mean HP<br>[GWh/d]',\n",
    "    'Firm HP [GWh/d]' : 'Firm HP<br>[GWh/d]'\n",
    "}, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39042bbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "merged_df['Scenario, [1/0]'] = merged_df['Scenario, [1/0]'].astype('float')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "8e2336d1",
   "metadata": {},
   "source": [
    "merged_df.rename(columns={\n",
    "    'Land loss, [km2]' : 'Land loss<br>[km2]'\n",
    "}, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58d60ed8",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(merged_df.columns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06383545",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b168f827",
   "metadata": {},
   "outputs": [],
   "source": [
    "mya_nobuilt_vis_simple = OutputVisualiser(merged_df[merged_df['Scenario, [1/0]'] == 1])\n",
    "fig1 = mya_nobuilt_vis_simple.plot_parallel(\n",
    "    columns = (\n",
    "        'Mean annual HP<br>[GWh/d]', \n",
    "        'Firm HP<br>[GWh/d]', \n",
    "        'GHG emissions<br>[Mt CO2,eq/year]',\n",
    "        'Agricultural land<br>loss, [km2]', \n",
    "        'Deforestation<br>[km2]'),\n",
    "    labels = {\n",
    "        'Mean HP<br>[GWh/d]' : 'Mean HP',\n",
    "        'Firm HP<br>[GWh/d]' : 'Firm HP',\n",
    "        'GHG emissions<br>[Mt CO2,eq/year]': 'GHG emissions',\n",
    "        'Agricultural land<br>loss, [km2]' : 'Agricultural land loss',\n",
    "        'Deforestation<br>[km2]' : 'Deforestation'},\n",
    "    color_col = 'Biogenic GHG intensity<br>[gCO2,eq/kWh]', color_limits=(0,60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de2874af",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1.write_image(pathlib.Path('figures/moo/parallel_plot.png'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "850c9856",
   "metadata": {},
   "outputs": [],
   "source": [
    "mya_nobuilt_vis = OutputVisualiser(merged_df) #merged_df)\n",
    "fig = mya_nobuilt_vis.plot_parallel(\n",
    "    columns = (\n",
    "        'Mean HP<br>[GWh/d]', \n",
    "        'Firm HP<br>[GWh/d]', \n",
    "        'GHG emissions<br>[Mt CO2,eq/year]',\n",
    "        'Agricultural land<br>loss, [km2]', \n",
    "        'Deforestation<br>[km2]', \n",
    "        'Land loss<br>[km2]', \n",
    "        'GHG intensity<br>[gCO2,eq/kWh]',\n",
    "        'Firm Power Ratio, [%]',\n",
    "        'Scenario, [1/0]'), \n",
    "    labels = {\n",
    "        'Mean HP<br>[GWh/d]' : 'Mean HP',\n",
    "        'Firm HP<br>[GWh/d]' : 'Firm HP',\n",
    "        'GHG emissions<br>[Mt CO2,eq/year]': 'GHG emissions',\n",
    "        'Agricultural land<br>loss, [km2]' : 'Agricultural land loss',\n",
    "        'Deforestation<br>[km2]' : 'Deforestation',\n",
    "        'Land loss<br>[km2]' : 'Land loss',\n",
    "        'GHG intensity<br>[gCO2,eq/kWh]' : 'GHG intensity',\n",
    "        'Firm Power Ratio, [%]' : 'Firm power ratio',\n",
    "        'Scenario, [1/0]' : 'Scenario, [1/0]'\n",
    "    },\n",
    "    color_col = 'Biogenic GHG intensity<br>[gCO2,eq/kWh]', color_limits=(0,60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4996f7cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_html_path=r\"figures/moo/parallel_plot.html\"\n",
    "input_template_path = r\"figures/moo/parallel_plot_template.html\"\n",
    "\n",
    "plotly_jinja_data = {\"fig\":fig.to_html(full_html=False)}\n",
    "#consider also defining the include_plotlyjs parameter to point to an external Plotly.js as described above\n",
    "\n",
    "with open(output_html_path, \"w\", encoding=\"utf-8\") as output_file:\n",
    "    with open(input_template_path) as template_file:\n",
    "        j2_template = Template(template_file.read())\n",
    "        output_file.write(j2_template.render(plotly_jinja_data))"
   ]
  }
 ],
 "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
}