Notebook_6_create_reservoir_maps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ca1763f8",
   "metadata": {},
   "source": [
    "# Create maps of Myanmar with reservoirs and reservoir emissions\n",
    "### T. JANUS, email: tomasz.k.janus@gmail.com , tomasz.janus@manchester.ac.uk\n",
    "#### Last modified: 27-June-2024\n",
    "\n",
    "## ALSO:\n",
    "### 1. Processes discounted emissions from multipurpose reservoirs by taking into account emissions associated with agricultural use\n",
    "### 2. Incorporates emissions of RoR hydropower by introducing a common emission intensity per generated energy\n",
    "### 3. Creates a combined dataframe for visualisation purposes and as a source of information for subsequent optimal dam selection.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15c3629d",
   "metadata": {},
   "source": [
    "### Load libraries and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90180aaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load libraries\n",
    "library(tidyverse)\n",
    "library(dplyr)\n",
    "library(sf)\n",
    "library(Hmisc)\n",
    "library(ggplot2)\n",
    "library(colorspace)\n",
    "library(ggspatial)\n",
    "library(latex2exp)\n",
    "library(readxl)\n",
    "library(ggnewscale)\n",
    "library(writexl) \n",
    "library(ggthemes)\n",
    "library(ggExtra)\n",
    "library(ggbreak)\n",
    "\n",
    "mya_epsg_code <- \"EPSG:32646\"\n",
    "googlecrs <- \"EPSG:4326\"\n",
    "\n",
    "text_fontsize <- 12\n",
    "axis_label_fontsize <- 13\n",
    "legend_text_fontsize <- 10\n",
    "legend_title_fontsize <- 12\n",
    "plot_title_fontsize <- 15\n",
    "strip_fontsize <- 13\n",
    "\n",
    "# Assumed emission intensity (per generated energy) of RoR hydroelectric plants\n",
    "em_intensity <- 3 # gCO2e/kWh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09ec2e9b",
   "metadata": {},
   "source": [
    "## UNITS\n",
    "1. Net emission fluxes (aka. net aerial emissions) in gCO$_{2e}$/m$^2$/year\n",
    "2. Emission intensities of HP in gCO$_{2e}$/kWh, e.g. 3~gCO$_{2e}$/kWh\n",
    "3. Total (annual) emissions in ktCO$_{2e}$year\n",
    "\n",
    "## VARIABLES\n",
    "* `co2_net`, `ch4_net`, `n2o_mean` are aerial emissions (fluxes) in gCO$_{2e}$/m$^2$/yr\n",
    "* `tot_em_int` represents fluxes (misnamed intensity) in the same unit as `co2_net` etc.\n",
    "* `tot_em` are total emissions in ktCO$_{2e}$/year\n",
    "* `em_intensity` is the emission intensity of hydropower in (gCO$_{2e}$/kWh) - same as tCO$_{2e}$/GWh\n",
    "* `ann_gen` is annual hydropower generation in (GWh/yr)\n",
    "* `mean` is the mean hydropower generation in (MW)\n",
    "\n",
    "## CONVERSIONS\n",
    "- `tot_em` (ktCO$_{2e}$/yr)     <- `em_intensity` (gCO$_{2e}$kWh) * `ann_gen` (GWh/yr) / 1000\n",
    "- `ann_gen` (GWh/yr)            <- `mean` (MW) * 24 * 365.25 / 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d591769b",
   "metadata": {},
   "source": [
    "## 1. PROCESS RESERVOIRS ONLY (no RoRs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2cb4b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the Myanmar country outline shape file and create a box (margins) around it\n",
    "mya_outline <- sf::st_read(\n",
    "    dsn = file.path(\"bin\",\"gis_layers\",\"myanmar_outline\",\"Myanmar_outline.shp\"), \n",
    "    quiet = TRUE, geometry_column = 'geometry') %>% \n",
    "    st_transform(mya_epsg_code)\n",
    "\n",
    "df_sf <- sf::st_as_sf(mya_outline, coords=c(\"Lon\", \"Lat\"), crs=googlecrs, agr = \"identity\")\n",
    "df_xy <- sf::st_transform(df_sf, crs=mya_epsg_code)\n",
    "bbox <- sf::st_bbox(df_xy)\n",
    "# Expand box by 20% to give a little extra room\n",
    "Dx <- (bbox[[\"xmax\"]]-bbox[[\"xmin\"]])*0.1\n",
    "Dy <- (bbox[[\"ymax\"]]-bbox[[\"ymin\"]])*0.1\n",
    "bbox[\"xmin\"] <- bbox[\"xmin\"] - Dx\n",
    "bbox[\"xmax\"] <- bbox[\"xmax\"] + Dx\n",
    "bbox[\"ymin\"] <- bbox[\"ymin\"] - Dy\n",
    "bbox[\"ymax\"] <- bbox[\"ymax\"] + Dy\n",
    "# \n",
    "bb <- c(bbox[\"xmin\"], bbox[\"ymin\"], bbox[\"xmax\"], bbox[\"ymax\"])\n",
    "\n",
    "# Load and filter the shape data with reservoirs, rivers and neighbouring countries\n",
    "mya_reservoirs <- sf::st_read(\n",
    "    dsn = file.path(\"bin\", \"heet_outputs_MIN_LOW_PRI\",\"reservoirs_updated.shp\"), quiet = TRUE) %>% \n",
    "    st_transform(mya_epsg_code) %>% \n",
    "    filter(!is.na(tot_em)) %>% \n",
    "    mutate(type = replace(type, type == \"unknown\", \"hydroelectric\")) %>%\n",
    "    rename(reservoir_type = type) %>%\n",
    "    # Get some names in mya_reservoirs to the IFC format (some names are off as they were not \n",
    "    # coded properly in HEET)\n",
    "    mutate(name = case_when(\n",
    "        name == \"Paung Laung (lower)\" ~ \"Lower Paunglaung\", \n",
    "        name == \"Paung Laung (upper)\" ~ \"Upper Paunglaung\",\n",
    "        name == \"Thapanzeik\" ~ \"Thaphanseik\",\n",
    "        TRUE ~ name\n",
    "    ))\n",
    "# Arange rows in res_centroids (reemission output data with geometry) from highest to lowest area\n",
    "res_centroids <- sf::st_centroid(mya_reservoirs, of_largest_polygon = TRUE) %>% arrange(desc(r_area_km2))\n",
    "# Load shape files with rivers and neigbouring countries\n",
    "mya_rivers <- sf::st_read(dsn = file.path(\"bin\",\"gis_layers\",\"mya_rivers.shp\"), quiet = TRUE)  %>% \n",
    "    st_transform(mya_epsg_code)\n",
    "mya_neighbours <- sf::st_read(\n",
    "    dsn = file.path(\"bin\",\"gis_layers\",\"World_Countries_Generalized.shp\"), quiet = TRUE)  %>% \n",
    "    st_transform(mya_epsg_code) %>%\n",
    "    filter(COUNTRY != \"Myanmar\")\n",
    "formatted_string <- paste(\n",
    "    \"Getting a cropped map with \",\n",
    "    nrow(res_centroids), \n",
    "    \"hydroelectric, multipurpose and irrigation reservoris.\")\n",
    "print(formatted_string)\n",
    "mya_neighbours_crop <- st_crop(mya_neighbours, y = bbox)\n",
    "num_hydroelectric = nrow(res_centroids %>% filter(reservoir_type == \"hydroelectric\"))\n",
    "num_multipurpose = nrow(res_centroids %>% filter(reservoir_type == \"multipurpose\"))\n",
    "num_irrigation = nrow(res_centroids %>% filter(reservoir_type == \"irrigation\"))\n",
    "formatted_string <- paste(\n",
    "    \"Number of hydroelectric reservoirs: \",\n",
    "    num_hydroelectric, \n",
    "    \". Number of multipurpose reservoirs: \",\n",
    "    num_multipurpose,\n",
    "    \". Number of irrigation reservoirs: \",\n",
    "    num_irrigation)\n",
    "print(formatted_string)\n",
    "# Calculate total emissions (in thousands of tonnes CO2e/year)\n",
    "res_centroids$tot_em_int = res_centroids$tot_em * res_centroids$r_area_km2 / 1000\n",
    "print(\"Reservoir types:\")\n",
    "print(unique(res_centroids$reservoir_type))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ac60d6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Hmisc::describe(res_centroids$tot_em_int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7ae7814",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(res_centroids,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d6472c7",
   "metadata": {},
   "source": [
    "# 2. CREATE A FACET MAP OF EMISSION INTENSITIES WITH ALL THREE RESERVOIR TYPES"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b093ac3a",
   "metadata": {},
   "source": [
    "#### Hints\n",
    "* to remove basemap layer from the map, comment out `ggspatial::annotation_map_tile(\n",
    "        type = \"osm\", zoom = 6, alpha=0.4, interpolate = TRUE, quiet = TRUE) +`\n",
    "* increase alpha in `geom_sf(data=mya_neighbours_crop, fill=\"lightgrey\", col=\"black\", alpha=0.1) +` from 0.1 to 0.5\n",
    "* increase alpha in `geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=0.8, size = 5.5)` to 0.7 to 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c81eca6c",
   "metadata": {},
   "source": [
    "## 2A. Maps with net (per area) GHG emissions for all three types of reservoirs - hydroelectric, multipurpose and irrigation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f257c6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin net (per-surface-area) emissions\n",
    "em_ranges <- c(\n",
    "    \"< 500\", \n",
    "    \"500 - 1000\", \n",
    "    \"1000 - 2000\", \n",
    "    \"2000 - 3000\",\n",
    "    \"3000 - 5000\",\n",
    "    \"> 5000\"\n",
    ")\n",
    "res_centroids <- \n",
    "    res_centroids %>%\n",
    "    mutate(\n",
    "        em_range = factor(case_when(\n",
    "            `tot_em` < 500 ~ em_ranges[1],\n",
    "            between(`tot_em`, 500, 1000) ~ em_ranges[2],\n",
    "            between(`tot_em`, 1000, 2000) ~ em_ranges[3],\n",
    "            between(`tot_em`, 2000, 3000) ~ em_ranges[4],\n",
    "            between(`tot_em`, 3000, 5000) ~ em_ranges[5],\n",
    "            `tot_em` >= 5000 ~ em_ranges[6]\n",
    "        ), levels = em_ranges)\n",
    "    )\n",
    "colors_map <- colorspace::diverging_hcl(6, palette = \"Tropic\", rev = F) # OK: Cork\n",
    "# Plot with ggplot2\n",
    "mya_gg <- ggplot(lims_method = \"geometry_bbox\")\n",
    "mya_gg +\n",
    "    geom_sf(fill = \"lightgrey\", alpha=0.9) +\n",
    "    geom_sf(data=mya_neighbours_crop, fill=\"lightgrey\", col=\"black\", alpha=0.5) +\n",
    "    geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=0.8, size = 5.5) +\n",
    "    geom_sf_label(data=mya_neighbours_crop, aes(label = COUNTRY), fill = \"white\") + \n",
    "    geom_sf(data=mya_rivers, col=\"lightblue3\") +\n",
    "    geom_sf(data = res_centroids %>% arrange(desc(r_area_km2)),\n",
    "          pch = 21,\n",
    "          aes(size = r_area_km2, fill = em_range),\n",
    "          col = \"Gray60\") +\n",
    "    scale_size(breaks = c(10, 50, 100, 300, 500),\n",
    "        range = c(1, 15),\n",
    "        guide = guide_legend(\n",
    "          direction = \"horizontal\",\n",
    "          ncol = 5)) +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map, \n",
    "        name=TeX(\n",
    "            r\"(\\overset{\\normalsize{Net aerial emission}}{\\overset{\\normalsize{gCO$_{2e}$/m$^2$/yr}}})\")) +\n",
    "    labs(#title = \"GHG emission intensities per reservoir type\",\n",
    "           sub = \"FutureDAMS Project\",\n",
    "           #caption = \"T. Janus, 2023\",\n",
    "           x = NULL, y = NULL,\n",
    "           size = TeX(r\"(\\overset{\\normalsize{Reservoir Area}}{\\overset{\\normalsize{km$^2$}}})\")) +\n",
    "    xlab(NULL) + ylab(NULL) +\n",
    "    annotate(geom = \"text\", x = 820000, y = 1500000, label = \"Andaman\\nSea\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    annotate(geom = \"text\", x = 505000, y = 1950000, label = \"Bay\\nof\\nBengal\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    ggspatial::annotation_scale(\n",
    "        location = \"tl\",\n",
    "        bar_cols = c(\"grey60\", \"white\"),\n",
    "        text_family = \"Avenir Next Condensed\") +\n",
    "    ggspatial::annotation_north_arrow(\n",
    "        location = \"tl\", which_north = \"true\",\n",
    "        pad_x = unit(0.12, \"in\"), pad_y = unit(0.3, \"in\"),\n",
    "        style = ggspatial::north_arrow_fancy_orienteering(\n",
    "          fill = c(\"grey40\", \"white\"),\n",
    "          line_col = \"grey20\",\n",
    "          text_family = \"Avenir Next Condensed\"))+\n",
    "    facet_wrap(~reservoir_type) +\n",
    "    xlim(bbox[\"xmin\"]+Dx, bbox[\"xmax\"]-Dx) + ylim(bbox[\"ymin\"]+Dy, bbox[\"ymax\"]-Dy) +\n",
    "    coord_sf(crs = mya_epsg_code) +\n",
    "    theme_minimal() +\n",
    "    theme(\n",
    "        text = element_text(\"Avenir Next Condensed\", size=text_fontsize),\n",
    "        legend.text = element_text(size = legend_text_fontsize),\n",
    "        legend.title = element_text(size=legend_title_fontsize),\n",
    "        axis.text.x = element_text(size = axis_label_fontsize),\n",
    "        axis.text.y = element_text(size = axis_label_fontsize),\n",
    "        strip.text = element_text(face=\"bold\", hjust = 0, size=strip_fontsize),\n",
    "        plot.title = element_text(hjust = 0, size=plot_title_fontsize),\n",
    "        plot.caption = element_text(face = \"italic\"),\n",
    "        panel.grid.major = element_line(\n",
    "            colour = gray(0.5), linetype = \"dashed\", size = 0.5), \n",
    "        panel.background = element_rect(fill = \"aliceblue\"), \n",
    "        panel.border = element_rect(fill = NA),\n",
    "        panel.spacing = unit(1, \"lines\"),\n",
    "        legend.direction = \"horizontal\", legend.position = \"bottom\",\n",
    "        panel.ontop = FALSE\n",
    "    )\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"emission_intensities_per_res_type.svg\"), \n",
    "    width = 27, height = 20, units = \"cm\")\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"emission_intensities_per_res_type.png\"), \n",
    "    width = 27, height = 20, units = \"cm\")\n",
    "#ggsave(\"map-us-ggdraw.pdf\", width = 5, height = 8) # Save to a pdf\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a54cc848",
   "metadata": {},
   "source": [
    "## 2B. Maps with total emissions in ktCO$_{2e}$/yr for all three types of reservoirs - hydroelectric, multipurpose and irrigation - RoRs not included"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14aa5c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot histogram of total emissions of reservoirs (includes hydroelectric, multipurpose and irrigation reservoirs)\n",
    "# tot_em_int are in ktCO2e/year\n",
    "hist(res_centroids$tot_em_int, breaks=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "536afe3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a categorical total emission variable by binning total emissions (in thousands of tonnes CO2e/yr)\n",
    "# into ranges of values defined in the vector 'tot_em_ranges'\n",
    "tot_em_ranges <- c(\n",
    "    \"< 3\", \n",
    "    \"3 - 30\",  \n",
    "    \"30 - 50\",\n",
    "    \"50 - 200\",\n",
    "    \"200 - 400\",\n",
    "    \"> 400\"\n",
    ")\n",
    "\n",
    "res_centroids <- \n",
    "    res_centroids %>%\n",
    "    mutate(\n",
    "        tot_em_range = factor(case_when(\n",
    "            `tot_em_int` < 3 ~ tot_em_ranges[1],\n",
    "            between(`tot_em_int`, 3, 30) ~ tot_em_ranges[2],\n",
    "            between(`tot_em_int`, 30, 50) ~ tot_em_ranges[3],\n",
    "            between(`tot_em_int`, 50, 200) ~ tot_em_ranges[4],\n",
    "            between(`tot_em_int`, 200, 400) ~ tot_em_ranges[5],\n",
    "            `tot_em_int` >= 400 ~ tot_em_ranges[6]\n",
    "        ), levels = tot_em_ranges)\n",
    "    )\n",
    "\n",
    "colors_map <- colorspace::diverging_hcl(6, palette = \"Tropic\", rev = F) # OK: Cork\n",
    "\n",
    "mya_gg <- ggplot(lims_method = \"geometry_bbox\")\n",
    "mya_gg +\n",
    "    geom_sf(fill = \"lightgrey\", alpha=0.9) +\n",
    "    geom_sf(data=mya_neighbours_crop, fill=\"lightgrey\", col=\"black\", alpha=0.5) +\n",
    "    geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=0.8, size = 5.5) +\n",
    "    geom_sf_label(data=mya_neighbours_crop, aes(label = COUNTRY), fill = \"white\") + \n",
    "    geom_sf(data=mya_rivers, col=\"lightblue3\") +\n",
    "    geom_sf(data = res_centroids %>% arrange(desc(r_area_km2)),\n",
    "          pch = 21,\n",
    "          aes(size = r_area_km2, fill = tot_em_range),\n",
    "          col = \"Gray60\") +\n",
    "    scale_size(breaks = c(10, 50, 100, 300, 500),\n",
    "        range = c(1, 15),\n",
    "        guide = guide_legend(\n",
    "          direction = \"vertical\",\n",
    "          nrow = 1)) +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map, \n",
    "        guide = guide_legend(\n",
    "            direction = \"vertical\",\n",
    "            nrow = 3),\n",
    "        name=TeX(r\"(\\overset{\\normalsize{Total emission}}{\\overset{\\normalsize{ktCO$_{2e}$/yr}}})\")) +\n",
    "    labs(\n",
    "        #title = \"GHG emissions per reservoir type\",\n",
    "        sub = \"FutureDAMS Project\",\n",
    "        #caption = \"T. Janus, 2023\",\n",
    "        size = TeX(r\"(\\overset{\\normalsize{Reservoir Area}}{\\overset{\\normalsize{km$^2$}}})\")) +\n",
    "    xlab(NULL) + ylab(NULL) +\n",
    "    theme_minimal() +\n",
    "    theme(\n",
    "        text = element_text(\"Avenir Next Condensed\", size=text_fontsize),\n",
    "        legend.text = element_text(size = legend_text_fontsize),\n",
    "        legend.title = element_text(size=legend_title_fontsize),\n",
    "        legend.key.size = unit(1,\"line\"),\n",
    "        axis.text.x = element_text(size = axis_label_fontsize),\n",
    "        axis.text.y = element_text(size = axis_label_fontsize),\n",
    "        strip.text = element_text(face=\"bold\", hjust = 0, size=strip_fontsize),\n",
    "        plot.title = element_text(hjust = 0, size=plot_title_fontsize),\n",
    "        plot.caption = element_text(face = \"italic\"),\n",
    "        panel.grid.major = element_line(\n",
    "            colour = gray(0.5), linetype = \"dashed\", size = 0.5), \n",
    "        panel.background = element_rect(fill = \"aliceblue\"), \n",
    "        panel.border = element_rect(fill = NA),\n",
    "        legend.direction = \"hozirontal\", legend.position = \"bottom\",\n",
    "        panel.spacing = unit(0.7, \"lines\"),\n",
    "        panel.ontop = FALSE\n",
    "    ) +\n",
    "    annotate(geom = \"text\", x = 820000, y = 1500000, label = \"Andaman\\nSea\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    annotate(geom = \"text\", x = 505000, y = 1950000, label = \"Bay\\nof\\nBengal\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    ggspatial::annotation_scale(\n",
    "        location = \"tl\",\n",
    "        bar_cols = c(\"grey60\", \"white\"),\n",
    "        text_family = \"ArcherPro Book\") +\n",
    "    ggspatial::annotation_north_arrow(\n",
    "        location = \"tl\", which_north = \"true\",\n",
    "        pad_x = unit(0.12, \"in\"), pad_y = unit(0.3, \"in\"),\n",
    "        style = ggspatial::north_arrow_fancy_orienteering(\n",
    "          fill = c(\"grey40\", \"white\"),\n",
    "          line_col = \"grey20\",\n",
    "          text_family = \"ArcherPro Book\"))+\n",
    "    facet_wrap(~reservoir_type) +\n",
    "    xlim(bbox[\"xmin\"]+Dx, bbox[\"xmax\"]-Dx) + ylim(bbox[\"ymin\"]+Dy, bbox[\"ymax\"]-Dy) +\n",
    "    coord_sf(crs = mya_epsg_code)\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\", \"total_emissions_per_res_type.svg\"), width = 27, height = 20, units = \"cm\")\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\", \"total_emissions_per_res_type.png\"), width = 27, height = 20, units = \"cm\")\n",
    "#ggsave(\"map-us-ggdraw.pdf\", width = 5, height = 8)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8875e5ed",
   "metadata": {},
   "source": [
    "## 2C. Plot a minimal map for visual abstract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec8b8d85",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors_map <- colorspace::diverging_hcl(6, palette = \"Tropic\", rev = F) # OK: Cork\n",
    "\n",
    "mya_gg <- ggplot(lims_method = \"geometry_bbox\")\n",
    "mya_gg +\n",
    "    geom_sf(fill = \"lightgrey\", alpha=0.9) +\n",
    "    geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=1, size = 5.5) +\n",
    "    geom_sf(data=mya_rivers, col=\"lightblue3\") +\n",
    "    geom_sf(data = res_centroids,\n",
    "          pch = 24,\n",
    "          size=2.5,\n",
    "          alpha = 0.7,\n",
    "          aes(fill = tot_em_range),\n",
    "          col = \"black\") +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map, \n",
    "        name=TeX(r\"(\\overset{\\normalsize{Total emission}}{\\overset{\\normalsize{ktCO$_{2e}$/yr}}})\")) +\n",
    "    labs(#caption = \"T. Janus, 2023\",\n",
    "           size = \"Reservoir Area\\nkm2\") +\n",
    "    xlab(NULL) + ylab(NULL) +\n",
    "    theme_void() +\n",
    "    theme(\n",
    "        text = element_text(\"Avenir Next Condensed\"),\n",
    "        strip.text = element_text(face=\"bold\", hjust = 0),\n",
    "        plot.title = element_text(hjust = 0),\n",
    "        plot.caption = element_text(face = \"italic\"),\n",
    "        panel.background = element_blank(),\n",
    "        panel.border = element_blank(),\n",
    "        legend.direction = \"horizontal\", legend.position = \"none\",\n",
    "        panel.ontop = FALSE\n",
    "    ) +\n",
    "    xlim(bbox[\"xmin\"]+Dx, bbox[\"xmax\"]-Dx) + ylim(bbox[\"ymin\"]+Dy, bbox[\"ymax\"]-Dy) +\n",
    "    coord_sf(crs = mya_epsg_code)\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_vis_abstract.png\"), width = 27, height = 20, units = \"cm\")\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_vis_abstract.svg\"), width = 27, height = 20, units = \"cm\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59a5d9bc",
   "metadata": {},
   "source": [
    "## 2D. Calculate total emissions by status (existing and future) including RoR hydroelectric\n",
    "### NOTE (reminder): \n",
    "**tot_em_int** is a total emission of a reservoir in ktCO$_{2e}$/yr\n",
    "\n",
    "**tot_em** is a unit per area emission in gCO$_{2e}$/m$^2$/yr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02881c85",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Import ifc database that contains information about statuses of hydroelectric and multipurpose dams\n",
    "# Why are we filtering out Bult dams?\n",
    "# Change statuses to Built and Future\n",
    "ifc_dams <- sf::st_read(\n",
    "    dsn = file.path(\"bin\", \"gis_layers\", \"ifc_database\", \"all_dams_replaced_refactored.shp\"), \n",
    "    quiet = TRUE)  %>% \n",
    "    st_transform(mya_epsg_code) %>%\n",
    "    filter(Status != \"Built\") %>%\n",
    "    st_drop_geometry() %>%\n",
    "    select(\"IFC_ID\", \"Status\", \"DAM_NAME\") %>%\n",
    "    mutate(Status = ifelse(\n",
    "        Status %in% c(\n",
    "            'Construction', 'MOU', 'LocMoU', 'Identified', 'GOM Plan', 'JVA', 'MOA',\n",
    "            'Suspended', 'Covenant'), 'Future', Status))\n",
    "# Set status in res_centroids to Existing by default.\n",
    "res_centroids$Status = \"Existing\"\n",
    "\n",
    "# Join the centroids dataframe with ifc_dams and replace the status to Future for Future IFC dams\n",
    "res_centroids_joined = dplyr::left_join(res_centroids, ifc_dams, by = join_by(id == IFC_ID)) %>%\n",
    "  st_transform(mya_epsg_code) %>%\n",
    "  mutate(Status = ifelse(is.na(Status.y), Status.x, Status.y)) %>%\n",
    "  select(-Status.x, -Status.y) %>%\n",
    "  mutate(Status = ifelse(\n",
    "      Status == \"Future\", \"Future Hydroelectric\", Status)) %>%\n",
    "  mutate(Status = ifelse(\n",
    "      Status == \"Existing\", \"Existing Hydroelectric\", Status)) %>%\n",
    "  arrange(desc(r_area_km2))\n",
    "\n",
    "# Find a list of future RoR hydroelectric from the IFC database\n",
    "# Future ROR\n",
    "ifc_dams_ror_future <- sf::st_read(\n",
    "    dsn = file.path(\"bin\",\"gis_layers\",\"ifc_database\",\"all_dams_replaced_refactored.shp\"), \n",
    "    quiet = TRUE) %>%\n",
    "    st_transform(mya_epsg_code) %>%\n",
    "    filter(Status != \"Built\") %>%\n",
    "    mutate(`RoR.or.Sto` = tolower(`RoR.or.Sto`)) %>%\n",
    "    filter(`RoR.or.Sto` == 'ror') %>%\n",
    "    mutate(Status = \"Future Hydroelectric\")\n",
    "# Existing ROR\n",
    "ifc_dams_ror_existing <- sf::st_read(\n",
    "    dsn = file.path(\"bin\",\"gis_layers\",\"ifc_database\",\"all_dams_replaced_refactored.shp\"), \n",
    "    quiet = TRUE) %>%\n",
    "    st_transform(mya_epsg_code) %>%\n",
    "    filter(Status == \"Built\") %>%\n",
    "    mutate(`RoR.or.Sto` = tolower(`RoR.or.Sto`)) %>%\n",
    "    filter(`RoR.or.Sto` == 'ror') %>%\n",
    "    mutate(Status = \"Existing Hydroelectric\")\n",
    "\n",
    "ifc_dams_ror_all <- rbind(ifc_dams_ror_future, ifc_dams_ror_existing) \n",
    "\n",
    "# Find if any future RoRs are already included in the reservoirs list for GHG emission estimation \n",
    "# (some RoRs were misclassified as RoR whereas in fact they're small storage dams)\n",
    "doubled_names_all <- dplyr::inner_join(res_centroids_joined %>% \n",
    "    st_drop_geometry(), ifc_dams_ror_all, by=join_by(name == DAM_NAME)) %>%\n",
    "    select(name) %>%\n",
    "    as.list()\n",
    "\n",
    "doubled_names_future <- dplyr::inner_join(res_centroids_joined %>% \n",
    "    st_drop_geometry(), ifc_dams_ror_future, by=join_by(name == DAM_NAME)) %>%\n",
    "    select(name) %>%\n",
    "    as.list()\n",
    "\n",
    "print(\"Doubled names of future dams\")\n",
    "print(doubled_names_future)\n",
    "\n",
    "ifc_dams_ror_all_filt <- ifc_dams_ror_all %>% filter(!(DAM_NAME %in% doubled_names_all$name))\n",
    "ifc_dams_ror_future_filt <- ifc_dams_ror_future %>% filter(!(DAM_NAME %in% doubled_names_future$name))\n",
    "\n",
    "# Units in the IFC data\n",
    "# Inst_cap - MW\n",
    "# Annual.Gen - GWh\n",
    "# Des_Head - m\n",
    "head(ifc_dams_ror_all_filt %>% select('DAM_NAME', 'Inst_cap', 'Annual.Gen', 'Des_Head'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19e7af1b",
   "metadata": {},
   "source": [
    "### `ifc_dams_ror_all_filt` is a dataframe of RoR hydroelectric unit with removed doubled entries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d3cc613",
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(res_centroids_joined)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0e8986c",
   "metadata": {},
   "outputs": [],
   "source": [
    "res_centroids_joined %>% filter(name==\"Thaphanseik\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f844b34",
   "metadata": {},
   "source": [
    "# 3. Calculate splits between emissions allocated to hydroelectric and to agricultural production in multipurpose reservoirs\n",
    "### ** The split usage is based on the amount of water going to turbine vs irrigation supply **\n",
    "- **TODO:** Check where the irrigation water should be abstracted from - Sedawgyi or Sedawgyi (upper)\n",
    "- **NOTE:** There are no MultiPurpose reservoirs in the Salween river basin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e8b1bc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Calculate flow proportions HP / Irrigation splits\n",
    "sittaung_averages <- read_excel(\n",
    "    file.path(\"outputs/pywr_sim_aligned_models/outputs_sittaung_hist_with_flow_rec.xlsx\")) %>%\n",
    "    slice(-c(1:2)) %>%       # Remove first two rows\n",
    "    select((ncol(.)-9):ncol(.)-1) %>% # Select last six columns\n",
    "    summarise_all(mean, na.rm = TRUE)\n",
    "pair_indices <- seq(1, ncol(sittaung_averages), by = 2)\n",
    "turbine_flows <- colSums(sittaung_averages[, pair_indices], na.rm = TRUE)\n",
    "irr_flows <- colSums(sittaung_averages[, pair_indices + 1], na.rm = TRUE)\n",
    "sittaung_hp_ratios <- turbine_flows / (irr_flows + turbine_flows)\n",
    "irrawaddy_averages <- read_excel(\n",
    "    file.path(\"outputs/pywr_sim_aligned_models/outputs_irrawaddy_hist_with_flow_rec.xlsx\")) %>%\n",
    "    slice(-c(1:2)) %>%       # Remove first two rows\n",
    "    select((ncol(.)-11):ncol(.)) %>% # Select last sixteen columns\n",
    "    summarise_all(mean, na.rm = TRUE)\n",
    "pair_indices <- seq(1, ncol(irrawaddy_averages), by = 2)\n",
    "turbine_flows <- colSums(irrawaddy_averages[, pair_indices], na.rm = TRUE)\n",
    "irr_flows <- colSums(irrawaddy_averages[, pair_indices + 1], na.rm = TRUE)\n",
    "irrawaddy_hp_ratios <- turbine_flows / (irr_flows + turbine_flows)\n",
    "\n",
    "df1 <- data.frame(sittaung_hp_ratios) %>% rename(hp_ratio = sittaung_hp_ratios)\n",
    "df2 <- data.frame(irrawaddy_hp_ratios) %>% rename(hp_ratio = irrawaddy_hp_ratios)\n",
    "hp_ratios <- rbind(df1, df2)\n",
    "new_row_names <- sub(\"_Turbine_Flow\", \"\", rownames(hp_ratios))\n",
    "rownames(hp_ratios) <- new_row_names\n",
    "hp_ratios$irr_ratio = 1 - hp_ratios$hp_ratio\n",
    "hp_ratios$name <- rownames(hp_ratios)\n",
    "\n",
    "print(\"HP/total flow proportions for  multipurpose reservoirs\")\n",
    "hp_ratios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df897142",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Fill the proportions manually and update the multipurpose reservoir table\n",
    "# The following three vectors are filled in manually from the values calculated in the previous step.\n",
    "# This step can be automated but due to time-constraints, we speed up by putting in the numbers by hand.\n",
    "reservoir_names <- c(\n",
    "    \"Kinda\", \"Myittha\", \"Myogyi\", \"Zawgyi II\", \"Sedawgyi\", \"Thaphanseik\",\n",
    "    \"Buywa (upper)\", \"Mone Chaung\", \"Kyee Ohn Kyee Wa\", \"Buywa\", \"Sedawgyi (upper)\",\n",
    "    \"Lower Paunglaung\", \"Kabaung\", \"Phyu Chaung\",\"Kun Chaung\", \"Yenwe\"\n",
    "    )\n",
    "river_basin <- c(\n",
    "    \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \n",
    "    \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \"irrawaddy\", \"irrawaddy\",\n",
    "    \"sittaung\", \"sittaung\", \"sittaung\", \"sittaung\", \"sittaung\"\n",
    "    )\n",
    "# Values calculated manually from the outputs (excel files from pywer simulations)\n",
    "percent_hp <- c(\n",
    "    0.67, 0.94, 0.83, 0.35, 0.76, 0.53,\n",
    "    1.0, 1.0, 1.0, 1.0, 1.0,\n",
    "    0.94, 0.69, 0.74, 0.82, 0.83)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaad838c",
   "metadata": {},
   "source": [
    "### Create a new dataframe `total_emissions_multi_merged` with emissions of multipurpose reservoirs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae2c5ba7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mone Chaung, Buywa, Buywa (upper) and Kyee Ohn Kyee Wa store water for mazali irrigation. \n",
    "# They have volumes of 0.61, 0.99, 4.25, and 0.21 (6.06 km3 in total)\n",
    "# let's water abstraction for irrigation based on relative volumes of each reservoirs and taking\n",
    "# irrigation water consumption from Mazali irrigation abstraction\n",
    "# CORRECTION: the irrigation flow is 5% of total flow. Therefore we have not used \n",
    "# The fractions computed are:\n",
    "# --------------------------------------------\n",
    "# vol_fractions <- c(0.61, 0.99, 4.25, 0.21)\n",
    "# vol_fractions / sum(vol_fractions)\n",
    "# --------------------------------------------\n",
    "# Result: 0.1, 0.163, 0.70, 0.035\n",
    "hp_fractions_multi <- data.frame(name = reservoir_names, basin = river_basin, hp_fraction = percent_hp)\n",
    "# Get the data from the IFC database that only lists multipurpose reservoirs\n",
    "hp_fractions_multi <- data.frame(name = reservoir_names, basin = river_basin, hp_fraction = percent_hp)\n",
    "# Get the data from the IFC database that only lists multipurpose reservoirs\n",
    "total_emissions_multi <- res_centroids_joined %>%\n",
    "  filter(!(DAM_NAME %in% doubled_names_all$name)) %>%\n",
    "  filter(reservoir_type == \"multipurpose\") %>%\n",
    "  select(name, Status, reservoir_type, r_volume_m, r_area_km2, r_mean_dep, tot_em, tot_em_int) %>%\n",
    "  mutate(volume_km3 = r_volume_m / 10^9) %>%\n",
    "  select(!r_volume_m)\n",
    "# Find emission splits between hp and irrigation for multipurpose reservoirs\n",
    "total_emissions_multi_merged <- merge(total_emissions_multi, hp_fractions_multi, by=\"name\") %>%\n",
    "    mutate(tot_em_hp=tot_em_int*hp_fraction, tot_em_irr=tot_em_int*(1-hp_fraction)) %>%\n",
    "    select(\n",
    "        name, Status, reservoir_type, r_area_km2, r_mean_dep, tot_em, volume_km3, \n",
    "        hp_fraction, tot_em_int, tot_em_hp, tot_em_irr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22bb24ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(total_emissions_multi_merged, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d1d8e75",
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_summary <- total_emissions_multi_merged %>%\n",
    "    group_by(Status) %>%\n",
    "    summarise(\n",
    "        emissions_total_hp = sum(tot_em_hp, na.rm = TRUE),\n",
    "        emissions_total_irr = sum(tot_em_irr, na.rm = TRUE),\n",
    "        emissions_total = sum(tot_em_int, na.rm = TRUE),\n",
    "        total_area = sum(r_area_km2, na.rm = TRUE),\n",
    "        total_volume = sum(volume_km3, na.rm = TRUE),\n",
    "        mean_depth = mean(r_mean_dep, na.rm = TRUE),\n",
    "        num = n()\n",
    "  ) \n",
    "# total_emissions_multi_merged provides unit net emissions of MULTIPURPOSE reservoirs \n",
    "# (per surface area) divided between hydroelectric and agricultural uses\n",
    "multi_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78fe1034",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Find emission statistics of MULTIPURPOSE RESERVOIRS\n",
    "emissions_multi_by_utility <- multi_summary %>%\n",
    "  select(-emissions_total) %>%\n",
    "  pivot_longer(cols = starts_with(\"emissions_total\"),\n",
    "               names_to = \"utility\",\n",
    "               values_to = \"emissions_total\") %>%\n",
    "  mutate(utility = ifelse(utility == \"emissions_total_hp\", \"hp\", \"irr\")) %>%\n",
    "  group_by(Status, utility) %>%\n",
    "  summarise(\n",
    "      n = num,\n",
    "      total_area = sum(total_area, na.rm = TRUE),\n",
    "      total_volume = sum(total_volume, na.rm = TRUE),\n",
    "      mean_depth = mean(mean_depth, na.rm = TRUE),\n",
    "      total_emission = sum(emissions_total, na.rm = TRUE)\n",
    "      )\n",
    "emissions_multi_by_utility"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dd44bfc",
   "metadata": {},
   "source": [
    "<b style='color:red;'>NOTE</b> the dataframe with discounted emissions of multipurpose reservoirs, i.e. `total_emissions_multi_merged` **has not been saved** <b style='color:red;'>!!!!!!!!!</b>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a6c25f1",
   "metadata": {},
   "source": [
    "# 4. Calculate total emissions and emisssions per GWh - include RoR and plot as triangles"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1eff126",
   "metadata": {},
   "source": [
    "## 4A. Calculate emissions of RoR plants using a constant emission intensity - see `em_intensity` variable and merged the calculated emissions and emission intensities into a dataframe with reservoir-based hydroelectric dams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d13bdaac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get emissions of RoR plants from pywr output excel file(s)\n",
    "ror_em_intensity <- em_intensity\n",
    "\n",
    "# Merge HP production figures into res_centroids data\n",
    "emissions_data <- readxl::read_excel(\n",
    "    file.path(\"intermediate\", \"merged_table.xlsx\")) %>%\n",
    "    distinct(ifc_id, .keep_all = TRUE) %>%\n",
    "    select(\n",
    "        ifc_id, dam_name, status_2_ifc, ro_r_or_sto_ifc, hp_type_reem, \n",
    "        res_area, mean, pctile_2, pctile_3, pctile_5, pctile_10) %>%\n",
    "    rename(name=dam_name) %>%\n",
    "    mutate(ann_gen = mean * 24 * 365.25 / 1000)\n",
    "    # Annual generation in GWh.year\n",
    "\n",
    "# Geospatial data from heet and reemission\n",
    "hp_dams <- res_centroids_joined %>%\n",
    "    filter(reservoir_type!='irrigation')\n",
    "\n",
    "# Merge the geospatial dataframe from heet and reemission with hp production estimates from \n",
    "# pywr (emissions data)\n",
    "# Represents only the delineated hydroelectric reservoirs, not ror plants\n",
    "joined_data <- merge(hp_dams, emissions_data, by = \"name\", all.x=TRUE) %>%\n",
    "    mutate(em_intensity = tot_em_int / ann_gen * 1000) %>%\n",
    "    mutate(\n",
    "        Status = ifelse(Status == 'Existing Reservoirs', 'Existing Hydroelectric', Status))\n",
    "#Hmisc::describe(joined_data$em_intensity)\n",
    "\n",
    "# Find HP generation for RoR plants from emissions_data from pywr filtered by hp type (sto/ror)\n",
    "emissions_ror <- emissions_data %>%\n",
    "    filter(hp_type_reem == 'ror') %>%\n",
    "    mutate(emission_tot_int = ann_gen * ror_em_intensity / 1000)\n",
    "\n",
    "# Create a dataframe of RoR hydropower with emissions and emission intensity information\n",
    "# x [GWh/year] * y [gCO2/kWh] * 1,000,000 kWh / GWh / 1,000,000 gCO2/tonneCO2 -> X[tCO2/year]\n",
    "# Divide by 1000 to get emission in ktCO2/year\n",
    "ror_with_emissions_gis <- \n",
    "    merge(x=ifc_dams_ror_all_filt, y=emissions_ror, by.x=\"IFC_ID\", by.y=\"ifc_id\") %>%\n",
    "    mutate(\n",
    "        tot_em_range = factor(case_when(\n",
    "            `emission_tot_int` < 3 ~ tot_em_ranges[1],\n",
    "            between(`emission_tot_int`, 3, 30) ~ tot_em_ranges[2],\n",
    "            between(`emission_tot_int`, 30, 50) ~ tot_em_ranges[3],\n",
    "            between(`emission_tot_int`, 50, 200) ~ tot_em_ranges[4],\n",
    "            between(`emission_tot_int`, 200, 400) ~ tot_em_ranges[5],\n",
    "            `emission_tot_int` >= 400 ~ tot_em_ranges[6]\n",
    "        ), levels = tot_em_ranges)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87f19323",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(ror_with_emissions_gis,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38bb224b",
   "metadata": {},
   "source": [
    "## 4B. Plot total emissions by status (existing and future) including RoR hydroelectric"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e7e3288",
   "metadata": {},
   "source": [
    "**NOTE** The map contains two \"dam\" layers - reservoirs, including multipurpose, hydroelectric and irrigation and a layer with RoR hydropower that is plotted with triangles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ef5320",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors_map <- colorspace::diverging_hcl(7, palette = \"Tropic\", rev = F) # OK: Cork\n",
    "\n",
    "mya_gg <- ggplot(lims_method = \"geometry_bbox\")\n",
    "mya_gg +\n",
    "    geom_sf(fill = \"lightgrey\", alpha=0.9) +\n",
    "    geom_sf(data=mya_neighbours_crop, fill=\"lightgrey\", col=\"black\", alpha=0.5) +\n",
    "    geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=0.8, size = 5.5) +\n",
    "    geom_sf_label(data=mya_neighbours_crop, aes(label = COUNTRY), fill = \"white\") + \n",
    "    geom_sf(data=mya_rivers, col=\"lightblue3\") +\n",
    "    geom_sf(data = res_centroids_joined,\n",
    "          pch = 21,\n",
    "          aes(size = r_area_km2, fill = tot_em_range),\n",
    "          col = \"grey31\") +\n",
    "    guides(fill = guide_legend(override.aes = list(size = 3))) +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map, \n",
    "        guide = guide_legend(ncol = 1),\n",
    "        name=TeX(r\"(\\overset{\\normalsize{Total emission}}{\\overset{\\normalsize{ktCO$_{2e}$/yr}}})\")) +\n",
    "    new_scale_fill() + \n",
    "    geom_sf(\n",
    "        data=ror_with_emissions_gis, pch = 24, color=\"grey31\", #fill=\"white\", \n",
    "        aes(fill=tot_em_range),\n",
    "        show_guide = FALSE,\n",
    "        alpha=0.8, size=2, lwd=0) +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map) + \n",
    "    scale_size(breaks = c(10, 50, 100, 200, 300, 500),\n",
    "        range = c(1, 15),\n",
    "        guide = guide_legend(\n",
    "          direction = \"vertical\",\n",
    "          ncol = 1)) +\n",
    "    labs(\n",
    "        #title = \"GHG emissions per construction status\",\n",
    "        sub = \"FutureDAMS Project\",\n",
    "        #caption = \"T. Janus, 2023\",\n",
    "        size = TeX(r\"(\\overset{\\normalsize{Reservoir Area}}{\\overset{\\normalsize{km$^2$}}})\")) +\n",
    "    xlab(NULL) + ylab(NULL) +\n",
    "    theme_minimal() +\n",
    "    theme(\n",
    "        text = element_text(\"Avenir Next Condensed\", size=text_fontsize),\n",
    "        legend.text = element_text(size = legend_text_fontsize),\n",
    "        legend.title = element_text(size=legend_title_fontsize),\n",
    "        legend.key.size = unit(1.3,\"line\"),\n",
    "        axis.text.x = element_text(size = axis_label_fontsize),\n",
    "        axis.text.y = element_text(size = axis_label_fontsize),\n",
    "        strip.text = element_text(face=\"bold\", hjust = 0, size=strip_fontsize),\n",
    "        plot.title = element_text(hjust = 0, size=plot_title_fontsize),\n",
    "        plot.caption = element_text(face = \"italic\"),\n",
    "        panel.grid.major = element_line(\n",
    "            colour = gray(0.5), linetype = \"dashed\", size = 0.5), \n",
    "        panel.background = element_rect(fill = \"aliceblue\"), \n",
    "        panel.border = element_rect(fill = NA),\n",
    "        legend.direction = \"vertical\", legend.position = \"right\",\n",
    "        panel.spacing = unit(1, \"lines\"),\n",
    "        panel.ontop = FALSE\n",
    "    ) +\n",
    "    annotate(geom = \"text\", x = 820000, y = 1500000, label = \"Andaman\\nSea\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    annotate(geom = \"text\", x = 505000, y = 1950000, label = \"Bay\\nof\\nBengal\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    ggspatial::annotation_scale(\n",
    "        location = \"tl\",\n",
    "        bar_cols = c(\"grey60\", \"white\"),\n",
    "        text_family = \"ArcherPro Book\") +\n",
    "    ggspatial::annotation_north_arrow(\n",
    "        location = \"tl\", which_north = \"true\",\n",
    "        pad_x = unit(0.12, \"in\"), pad_y = unit(0.3, \"in\"),\n",
    "        style = ggspatial::north_arrow_fancy_orienteering(\n",
    "          fill = c(\"grey40\", \"white\"),\n",
    "          line_col = \"grey20\",\n",
    "          text_family = \"ArcherPro Book\"))+\n",
    "    facet_wrap(~Status) +\n",
    "    xlim(bbox[\"xmin\"]+Dx, bbox[\"xmax\"]-Dx) + ylim(bbox[\"ymin\"]+Dy, bbox[\"ymax\"]-Dy) +\n",
    "    coord_sf(crs = mya_epsg_code)\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_per_status.png\"), width = 27, height = 20, units = \"cm\")\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_per_status.svg\"), width = 27, height = 20, units = \"cm\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80c44664",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(res_centroids_joined %>% select(reservoir_type, name, Status), 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4aec59d",
   "metadata": {},
   "source": [
    "**Emission intensity of HP :**\n",
    "* tot_emisssion $E_t$ [ktCO$_{2e}$/yr]\n",
    "* Flow, $Q$ in Mm$^3$/d\n",
    "* Head $H$ in metres\n",
    "\n",
    "Emission intensity is in tonnesCO$_{2e}$/GWh\n",
    "\n",
    "P = $\\eta \\rho Q \\, g \\, H$ \n",
    "\n",
    "Emission intensity $\\propto E_t / P$ = $ (E_{CO_2} + E_{CH_4}) / (Q \\, H)$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfc4ff14",
   "metadata": {},
   "source": [
    "## 4C. Print statistics about emissions of existing and future hydroelectric, irrigation and multipurpose reservoirs whilst taking into account partitioning of emissions into agricultural and hydroelectric use in multipurpose reservoirs via the split of water directed towards the two water uses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "306d4da3",
   "metadata": {},
   "outputs": [],
   "source": [
    "emissions_multi_by_utility <- res_centroids_joined %>%\n",
    "  group_by(Status, reservoir_type) %>%\n",
    "  summarise(\n",
    "      number = n(),\n",
    "      total_emission = sum(tot_em_int, na.rm = TRUE),\n",
    "      total_area = sum(r_area_km2, na.rm = TRUE),\n",
    "      total_volume = sum(r_volume_m, na.rm = TRUE) / 1e9,\n",
    "      mean_depth = mean(r_mean_dep, na.rm = TRUE)\n",
    "      ) %>%\n",
    "    select(!geometry)\n",
    "# Emissions of multipurpose reservoirs by utility\n",
    "# Legend:\n",
    "# emissions_total in ktCO2e/yr\n",
    "# total_area in km2\n",
    "# total_volume in km3\n",
    "# mean_depth in m\n",
    "print(\"Emissions of multipurpose reservoirs by utility\")\n",
    "head(emissions_multi_by_utility)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3708f85e",
   "metadata": {},
   "source": [
    "## 4D. Emissions per generated hydroelectric power of future and existing hydroelectric dams\n",
    "\n",
    "###  1. Combine RoR Multipurpose and HP data with name/id/area/volume/mean_depth/generation/em_intensity\n",
    "### 2. Visualise emission intensities on maps using intensity as colour and generation as size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea45e3c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform spatial left join with an sf dataframe containing multipurpose reservoirs with\n",
    "# hp_fraction that quantified the percentage of 'usage' of the reservoir to generate\n",
    "# hydroelectric energy. The rest of the water in the reservoir is used for agriculture\n",
    "# since the multipurpose reservoirs in Myanmar are HP/irrigation only. We discount emission\n",
    "# contributions to HP generation by the amount of water from the reservoir used for\n",
    "# agriculture\n",
    "# Set all multipurpose reservoirs to storage type\n",
    "total_emissions_multi_merged$hp_type_reem = \"sto\"\n",
    "\n",
    "ror_emissions_trimmed <- ror_with_emissions_gis %>% select(\n",
    "    c('name','status_2_ifc','ro_r_or_sto_ifc','hp_type_reem','res_area',\n",
    "       'mean','pctile_2','pctile_3','pctile_5','pctile_10','ann_gen',\n",
    "       'emission_tot_int','geometry','DAM_NAME', 'Status')\n",
    ") %>% rename(tot_em_int_sc=emission_tot_int) %>%\n",
    "mutate(reservoir_type = \"hydroelectric\")\n",
    "\n",
    "\n",
    "joined_data_with_hp_fraction <- joined_data %>% st_join(\n",
    "    total_emissions_multi_merged %>% select(name, hp_fraction),\n",
    "    by=c(\"name\", \"name\")) %>%\n",
    "    select(-ends_with(\".y\")) %>%\n",
    "    rename(name = name.x)\n",
    "# Scale total emissions for reservoirs in which hp_fraction is not <NA>\n",
    "joined_data_with_hp_fraction <- joined_data_with_hp_fraction %>%\n",
    "    mutate(\n",
    "        co2_net_sc = ifelse(\n",
    "            !is.na(hp_fraction), co2_net * hp_fraction, co2_net)\n",
    "    ) %>%\n",
    "    mutate(\n",
    "        ch4_net_sc = ifelse(\n",
    "            !is.na(hp_fraction), ch4_net * hp_fraction, ch4_net)\n",
    "    ) %>%\n",
    "    mutate(\n",
    "        tot_em_sc = ifelse(\n",
    "            !is.na(hp_fraction), tot_em * hp_fraction, tot_em)\n",
    "    ) %>%\n",
    "    mutate(\n",
    "        tot_em_int_sc = ifelse(\n",
    "            !is.na(hp_fraction), tot_em_int * hp_fraction, tot_em_int))\n",
    "\n",
    "# Include RoR in joined data\n",
    "joined_data_hp_multi_ror <- \n",
    "    bind_rows(joined_data_with_hp_fraction, ror_emissions_trimmed) %>%\n",
    "    select(-c(\n",
    "        'r_volume_m','r_maximum_','r_mean_dep','r_msocs_kg',\n",
    "        'r_mghr_all','r_mghr_may','r_mghr_nov','r_mean_ann','c_mar_mm',\n",
    "        'c_area_km2','n_populati','c_mean_slo','c_map_mm','c_mpet_mm',\n",
    "        'c_masm_mm','c_biome','c_soil_typ','c_mean_ols','ms_length',\n",
    "        'co2_net', 'ch4_net', 'n2o_mean', 'tot_em', 'tot_em_int',\n",
    "        'ifc_id', 'id', 'DAM_NAME', 'em_range','tot_em_range')) %>%\n",
    "    mutate(\n",
    "        reservoir_type = case_when(\n",
    "          is.na(reservoir_type) ~ \"hydroelectric\",\n",
    "          TRUE ~ reservoir_type)) %>%\n",
    "    mutate(\n",
    "        hp_fraction = case_when(\n",
    "          is.na(hp_fraction) ~ 1.0,\n",
    "          TRUE ~ hp_fraction)) %>%\n",
    "    mutate(Status = case_when(\n",
    "        is.na(Status) & status_2_ifc == \"E\" ~ \"Existing Hydroelectric\",\n",
    "        is.na(Status) ~ \"Future Hydroelectric\",\n",
    "        TRUE ~ Status\n",
    "      )) %>%\n",
    "    rename(\n",
    "        co2_net = co2_net_sc,\n",
    "        ch4_net = ch4_net_sc,\n",
    "        tot_em_net = tot_em_sc,\n",
    "        tot_em = tot_em_int_sc) %>% # Back-calculate emission intensities in tCO2e/GWh\n",
    "    mutate(em_intensity = tot_em / ann_gen * 1000) %>% # Re-evaluate em_intensity_range\n",
    "    # Fill res_area from r_area_km2 if r_area_km2 is not NA but res_area is NA\n",
    "    mutate(res_area = ifelse(!is.na(r_area_km2) & is.na(res_area), r_area_km2, res_area))\n",
    "\n",
    "# Correct some errors in data\n",
    "#joined_data_hp_multi_ror[joined_data_hp_multi_ror$name==\"Thaphanseik\", 'hp_type_reem'] <- \"sto\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370161fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(joined_data_hp_multi_ror,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a14550",
   "metadata": {},
   "source": [
    "`joined_data_hp_multi_ror` contains storage and RoR hydroelectric unit data including emissions and emission intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa0241d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "head(joined_data_hp_multi_ror,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ae05b80",
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(joined_data_hp_multi_ror)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12565365",
   "metadata": {},
   "source": [
    "**Descriptions of variables:**\n",
    "* `co2_net`, `ch4_net` and `tot_em_net` are net per area emissions in gCO$_{2e}$/m$^2$/yr for hydroelectric and multipurpose reservoirs, i.e. rows with non-zero (non-NA) res_area.\n",
    "* `ann_gen` is annual hydroelectric power generation in GWh/yr.\n",
    "* `mean` is the average daily hydropower production obtained from the water resources model implemented in Pywr. It is given in MW.\n",
    "* `tot_em_int` (or `tot_em` depending on the dataframe) is total emission in ktCO$_{2e}$/yr and is calculated from net (per-area) emission for reservoirs and from `ann_gen` for RoR hydropower by applying a constant emission intensity of a given value. Check `em_intensity` column for values. `em_intensity` is given in tCO$_{2e}$/GWh.\n",
    "\n",
    "**Conversions:**\n",
    "- `ann_gen` (GWh/yr)            <- `mean` (MW) * 24 * 365.25 / 1000\n",
    "- `tot_em` (ktCO$_{2e}$/yr)     <- `em_intensity` (tCO$_{2e}$/GWh) * `ann_gen` (GWh/yr) / 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67a5e9c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined_data_hp_multi_ror %>% group_by(Status, hp_type_reem) %>%\n",
    "  summarise(\n",
    "      number = n(),\n",
    "      total_generation = sum(ann_gen, na.rm = TRUE),\n",
    "      total_emission = sum(tot_em, na.rm = TRUE),\n",
    "      )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c27b5943",
   "metadata": {},
   "source": [
    "## 2G. Plot emission intensities of hydroelectric reservoirs in Myanmar including Hydroelectric Dams, Multipurpose Dams and RoR hydropower"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2473cde0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List multipurpose reservoirs\n",
    "joined_data_hp_multi_ror %>% \n",
    "    select(name, reservoir_type, hp_type_reem) %>% \n",
    "    filter(reservoir_type==\"multipurpose\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65540d63",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined_data_hp_multi_ror <- joined_data_hp_multi_ror %>% \n",
    "    mutate(\n",
    "        #tot_em_int = ifelse(hp_type_reem==\"ror\" & is.na(tot_em_int), em_intensity, tot_em_int),\n",
    "        tot_em = ifelse(hp_type_reem==\"ror\" & is.na(tot_em), em_intensity * ann_gen / 1000, tot_em)\n",
    "    )\n",
    "\n",
    "#filter(hp_type_reem==\"ror\") %>%  \n",
    "#    mutate(tot_em = ann_gen * em_intensity / 1000) %>%\n",
    "#    mutate(tot_em_int = em_intensity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4284857",
   "metadata": {},
   "outputs": [],
   "source": [
    "joined_data_hp_multi_ror"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f1d092b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin emissions intensities of hydroelectric plants in tonnesCO2eq/GWh\n",
    "em_intensity_ranges <- c(\n",
    "    \"< 5\",\n",
    "    \"5 - 20\", \n",
    "    \"20 - 50\", \n",
    "    \"50 - 150\", \n",
    "    \"150 - 300\",\n",
    "    \"300 - 450\",\n",
    "    \"450 - 1500\",\n",
    "    \"> 1500\"\n",
    ")\n",
    "\n",
    "joined_data_hp_multi_ror <- \n",
    "    joined_data_hp_multi_ror %>%\n",
    "    mutate(\n",
    "        em_intensity_range = factor(case_when(\n",
    "            `em_intensity` < 5 ~ em_intensity_ranges[1],\n",
    "            between(`em_intensity`, 5, 20) ~ em_intensity_ranges[2],\n",
    "            between(`em_intensity`, 20, 50) ~ em_intensity_ranges[3],\n",
    "            between(`em_intensity`, 50, 150) ~ em_intensity_ranges[4],\n",
    "            between(`em_intensity`, 150, 300) ~ em_intensity_ranges[5],\n",
    "            between(`em_intensity`, 300, 450) ~ em_intensity_ranges[6],\n",
    "            between(`em_intensity`, 450, 1500) ~ em_intensity_ranges[7],\n",
    "            `em_intensity` >= 1500 ~ em_intensity_ranges[8]\n",
    "        ), levels = em_intensity_ranges)\n",
    "    )\n",
    "\n",
    "unique(joined_data_hp_multi_ror$Status)\n",
    "hist(joined_data_hp_multi_ror$em_intensity, breaks = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f336157",
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(joined_data_hp_multi_ror$ann_gen, breaks = 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe571930",
   "metadata": {},
   "outputs": [],
   "source": [
    "write_xlsx(\n",
    "    joined_data_hp_multi_ror,\n",
    "    file.path(\"intermediate\", \"hp_multi_ror_emissions_and_generation.xlsx\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4a4e8bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "colors_map <- colorspace::diverging_hcl(8, palette = \"Tropic\", rev = F, power = 0.5) # OK: Cork\n",
    "joined_data_hp_multi_ror_plot <- joined_data_hp_multi_ror %>%\n",
    "     mutate(Status = replace(Status, Status == \"Existing Reservoirs\", \"Existing Hydroelectric\"))\n",
    "\n",
    "xy_coords <- st_coordinates(joined_data_hp_multi_ror_plot$geometry)\n",
    "\n",
    "mya_gg <- ggplot(lims_method = \"geometry_bbox\")\n",
    "mya_gg +\n",
    "    geom_sf(fill = \"lightgrey\", alpha=0.9) +\n",
    "    geom_sf(data=mya_neighbours_crop, fill=\"lightgrey\", col=\"black\", alpha=0.5) +\n",
    "    geom_sf(data=mya_outline, fill=\"antiquewhite1\", col=\"black\", alpha=0.8, size = 5.5) +\n",
    "    geom_sf_label(data=mya_neighbours_crop, aes(label = COUNTRY), fill = \"white\") + \n",
    "    geom_sf(data=mya_rivers, col=\"lightblue3\") +\n",
    "    geom_sf(\n",
    "        data = joined_data_hp_multi_ror_plot %>% \n",
    "            filter(hp_type_reem != 'ror') %>%\n",
    "            arrange(desc(ann_gen)),\n",
    "          pch = 21,\n",
    "          aes(\n",
    "              size = ann_gen, \n",
    "              fill = em_intensity_range,\n",
    "          ),\n",
    "          col = \"grey31\", alpha=0.8) +\n",
    "\n",
    "    ggrepel::geom_text_repel(\n",
    "        data = joined_data_hp_multi_ror_plot,\n",
    "        aes(label = name, x=xy_coords[,'X'], y=xy_coords[,'Y']), \n",
    "        size = 2, max.overlaps = 30, color = 'black',\n",
    "        segment.size=0.3) +\n",
    "\n",
    "    geom_sf(\n",
    "         data = joined_data_hp_multi_ror_plot %>% \n",
    "            filter(hp_type_reem == 'ror') %>%\n",
    "            arrange(desc(ann_gen)),\n",
    "          pch = 24,\n",
    "          aes(\n",
    "              size = ann_gen, \n",
    "              fill = em_intensity_range,\n",
    "          ),\n",
    "          legend=FALSE,\n",
    "          col = \"grey31\", alpha=0.8) +\n",
    "    guides(fill = guide_legend(override.aes = list(size = 3))) +\n",
    "    scale_fill_manual(\n",
    "        values = colors_map, \n",
    "        name=TeX(r\"(\\overset{\\normalsize{Emission intensity}}{\\overset{\\normalsize{gCO$_{2e}$/kWh}}})\"),\n",
    "        #name=\"Emission footprint\\ntonne CO2eq/GWh\",\n",
    "        guide = guide_legend(\n",
    "          direction = \"horizontal\",\n",
    "          ncol = 2)) +\n",
    "    scale_size(breaks = c(100, 500, 2000, 7000, 10000, 15000, 35000),\n",
    "        range = c(1, 12),\n",
    "        guide = guide_legend(\n",
    "          direction = \"horizontal\",\n",
    "          ncol = 3)) +\n",
    "    scale_shape_manual(\n",
    "        name=TeX(r\"(\\overset{\\normalsize{HP Type}}{\\overset{\\normalsize{Sto/RoR}}})\"),\n",
    "        values = c(\"sto\" = 21, \"ror\" = 24)) +\n",
    "    labs(\n",
    "        #title = expression(paste(\"GHG emissions per generated energy in dams over 5 M\", m^3)),\n",
    "        sub = \"FutureDAMS Project\",\n",
    "        fill=\"\",\n",
    "        #caption = \"T. Janus, 2023\",\n",
    "        size = TeX(r\"(\\overset{\\normalsize{Annual Generation}}{\\overset{\\normalsize{GWh/yr}}})\")) +\n",
    "    #guides(\n",
    "    #    fill = guide_legend(\n",
    "    #        override.aes = list(size = c(2.5,2.5), pch = 21)),\n",
    "    #    size = guide_legend(\n",
    "    #        override.aes = list(pch = 21))) +\n",
    "    xlab(NULL) + ylab(NULL) +\n",
    "    theme_minimal() +\n",
    "    theme(       \n",
    "        text = element_text(\"Avenir Next Condensed\", size=text_fontsize),\n",
    "        legend.text = element_text(size = legend_text_fontsize),\n",
    "        legend.title = element_text(size=legend_title_fontsize),\n",
    "        axis.text.x = element_text(size = axis_label_fontsize),\n",
    "        axis.text.y = element_text(size = axis_label_fontsize),\n",
    "        strip.text = element_text(face=\"bold\", hjust = 0, size=strip_fontsize),\n",
    "        plot.title = element_text(hjust = 0, size=plot_title_fontsize),\n",
    "        plot.caption = element_text(face = \"italic\"),\n",
    "        panel.grid.major = element_line(\n",
    "            colour = gray(0.5), linetype = \"dashed\", size = 0.5), \n",
    "        panel.background = element_rect(fill = \"aliceblue\"), \n",
    "        panel.border = element_rect(fill = NA),\n",
    "        legend.direction = \"vertical\", legend.position = \"right\",\n",
    "        panel.spacing = unit(1, \"lines\"),\n",
    "        panel.ontop = FALSE\n",
    "    ) +\n",
    "    annotate(geom = \"text\", x = 820000, y = 1500000, label = \"Andaman\\nSea\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    annotate(geom = \"text\", x = 505000, y = 1950000, label = \"Bay\\nof\\nBengal\", \n",
    "     fontface = \"italic\", color = \"grey22\", size = 4) +\n",
    "    ggspatial::annotation_scale(\n",
    "        location = \"tl\",\n",
    "        bar_cols = c(\"grey60\", \"white\"),\n",
    "        text_family = \"ArcherPro Book\") +\n",
    "    ggspatial::annotation_north_arrow(\n",
    "        location = \"tl\", which_north = \"true\",\n",
    "        pad_x = unit(0.12, \"in\"), pad_y = unit(0.3, \"in\"),\n",
    "        style = ggspatial::north_arrow_fancy_orienteering(\n",
    "          fill = c(\"grey40\", \"white\"),\n",
    "          line_col = \"grey20\",\n",
    "          text_family = \"ArcherPro Book\"))+\n",
    "    facet_wrap(~Status) +\n",
    "    xlim(bbox[\"xmin\"]+Dx, bbox[\"xmax\"]-Dx) + ylim(bbox[\"ymin\"]+Dy, bbox[\"ymax\"]-Dy) +\n",
    "    coord_sf(crs = mya_epsg_code)\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_per_GWh.png\"), width = 27, height = 20, units = \"cm\")\n",
    "ggsave(\n",
    "    file.path(\"figures\",\"maps\",\"total_emissions_per_GWh.svg\"), width = 27, height = 20, units = \"cm\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "448f9507",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same but from plot data\n",
    "joined_data_hp_multi_ror_plot %>%\n",
    "  filter(hp_type_reem == \"ror\") %>%\n",
    "  group_by(Status, hp_type_reem) %>%\n",
    "  summarise(\n",
    "      number = n(),\n",
    "      total_emission = sum(tot_em, na.rm = TRUE)\n",
    "      )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b78daf19",
   "metadata": {},
   "source": [
    "## 3. CREATE A STATISTICS PLOT TO COMPLEMENT THE EMISSION INTENSITY PLOT OF HYDROELECTRIC RESERVOIRS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "200b336e",
   "metadata": {},
   "source": [
    "It would be useful to :\n",
    "* **A.** check the relationship between per-area emissions and per HP emissions\n",
    "* **B.** Plot histograms of per GWh emissions of Existing/Future HP Units. Overlay a line at emission intensity of RoR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd0d235",
   "metadata": {},
   "outputs": [],
   "source": [
    "len <- length(joined_data_hp_multi_ror_plot$name)\n",
    "shapes <- rep_len(x = c(21, 22, 24), length.out = len)\n",
    "font_size_axis_xy <- 15\n",
    "font_size_labels_xy <- 17\n",
    "text_fontsize_xy <- 14\n",
    "tick_size <- 15\n",
    "library(scico)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7d06e0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plot_stat$hp_type_reem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "614914b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plot_stat <- joined_data_hp_multi_ror_plot %>% \n",
    "    filter(hp_type_reem == 'sto' | name == \"Upper Paunglaung\" | name == \"Lower Paunglaung\") %>% \n",
    "    arrange(desc(tot_em))\n",
    "\n",
    "data_plot_stat_ror <- joined_data_hp_multi_ror_plot %>% filter(hp_type_reem == 'ror') %>% arrange(desc(tot_em))\n",
    "df_ror <-data.frame(\n",
    "    tot_em_net=0,em_intensity=3, tot_em = 400, res_area = 100,\n",
    "    reservoir_type = \"ror\")\n",
    "\n",
    "# Select reservoirs to annotate and inspect\n",
    "# Points with extraordinarily large emission intensities\n",
    "points_large_em_int <- data_plot_stat %>% filter(em_intensity > 1000)\n",
    "\n",
    "points_large_net_emissions <- data_plot_stat %>% filter(tot_em_net>3000)\n",
    "\n",
    "points_interm_area_low_net_emissions <- data_plot_stat %>% filter(\n",
    "    tot_em_net < 800 & res_area > 200 & \n",
    "    res_area < 700 & em_intensity > 10\n",
    ")\n",
    "\n",
    "points_small_tot_em_small_area_interm_net_emissions <- data_plot_stat %>%\n",
    "    filter(\n",
    "        tot_em_net > 1700 & tot_em_net < 2000 & em_intensity < 1000 &\n",
    "        em_intensity > 10\n",
    "    )\n",
    "\n",
    "label_df <- bind_rows(\n",
    "    points_large_em_int, \n",
    "    points_large_net_emissions, \n",
    "    points_interm_area_low_net_emissions,\n",
    "    points_small_tot_em_small_area_interm_net_emissions\n",
    ")\n",
    "\n",
    "# Plot points\n",
    "base <- ggplot(\n",
    "    data = data_plot_stat, \n",
    "    aes(\n",
    "        tot_em_net, em_intensity, \n",
    "        size = tot_em, \n",
    "        fill=res_area, \n",
    "        shape=factor(reservoir_type),\n",
    "        color=factor(reservoir_type)\n",
    "        #order=-as.numeric(factor(joined_data_hp_multi_ror$hp_type_reem))\n",
    "    ), alpha=0.4\n",
    ") +\n",
    "geom_point(stroke=0.8) + \n",
    "# Add point representing emission intensity of RoR\n",
    "geom_point(data = df_ror, size = 8, legend=FALSE, shape=24, color=\"#2471a3\", fill=\"white\") +\n",
    "\n",
    "#ggforce::geom_mark_ellipse(aes(fill = res_area, filter = size > 1000)) +\n",
    "# Set labels\n",
    "xlab(TeX(\"Net aerial emission, gCO$_{2e}$/m$^2$/yr\")) +\n",
    "ylab(TeX(\"Emission intensity, gCO$_{2e}$/kWh\")) +\n",
    "# Add text labels to the selected points\n",
    "ggrepel::geom_text_repel(\n",
    "    seed=40,\n",
    "    data = label_df,\n",
    "    box.padding = 1.6,\n",
    "    aes(label = name), size = 4, max.overlaps = 30, color = 'grey35', fill=\"white\") +\n",
    "# Set scales\n",
    "scale_colour_manual(values = c(\"#d8b365\", \"#5ab4ac\", \"#2471a3\"), name=\"Reservoir Type\") +\n",
    "scale_size(\n",
    "    range = c(0, 15),\n",
    "    TeX(r\"(\\overset{Total Emission}{\\overset{\\normalsize{ktCO$_{2e}$/yr}}})\")) +\n",
    "scale_shape_manual(values = shapes, name=\"Reservoir Type\") +\n",
    "scale_fill_gradient(\n",
    "    low = \"white\", high = \"black\",\n",
    "    limits=c(0,1000),\n",
    "    name=TeX(r\"(\\overset{Reservoir Area}{\\overset{\\normalsize{km$^2$}}})\"))+\n",
    "scale_y_log10(limits = c(1,10000)) +  \n",
    "\n",
    "# Set the theme\n",
    "theme_bw() +\n",
    "theme(\n",
    "  text = element_text(size=text_fontsize_xy),\n",
    "  #panel.border = element_rect(fill=NA, colour = \"black\", size=0.1),\n",
    "  panel.border = element_blank(),\n",
    "  panel.grid.major = element_blank(),\n",
    "  panel.grid.minor = element_blank(),\n",
    "  strip.text.x = element_blank(),\n",
    "  strip.background = element_rect(colour=\"white\", fill=\"white\"),\n",
    "  axis.line=element_line(size=0.7),\n",
    "  axis.text = element_text(size = font_size_axis_xy),\n",
    "  axis.ticks=element_line(size=0.7),\n",
    "  axis.ticks.length=unit(.12, \"cm\"),\n",
    "  axis.title=element_text(size=font_size_labels_xy),\n",
    "    legend.direction = \"vertical\", legend.position = \"right\"\n",
    ") + \n",
    "guides(\n",
    "    shape = guide_legend(\n",
    "        override.aes = list(size = 5, width=12)),\n",
    "    linetype = guide_legend(override.aes = list(size = 5)),\n",
    "    color = guide_legend(\n",
    "        override.aes = list(size=5, width=12)),\n",
    "    size = guide_legend(\n",
    "        override.aes = list(color=\"gray88\"))\n",
    ")\n",
    "\n",
    "xbox <- cowplot::axis_canvas(base, axis = \"x\", coord_flip = TRUE) + \n",
    "  geom_boxplot(\n",
    "      data = data_plot_stat,\n",
    "      aes(\n",
    "          y = tot_em_net, \n",
    "          x = factor(reservoir_type), \n",
    "          color = factor(reservoir_type))) + \n",
    "    scale_x_discrete() +\n",
    "    scale_y_continuous(limits = c(0,4000)) +\n",
    "    scale_colour_manual(values = c(\"#d8b365\", \"#5ab4ac\", \"#2471a3\"), name=\"Reservoir Type\") + \n",
    "    coord_flip()\n",
    "ybox <- cowplot::axis_canvas(base, axis = \"y\") + \n",
    "  geom_boxplot(\n",
    "      data = data_plot_stat, \n",
    "      aes(\n",
    "          y = em_intensity, \n",
    "          x = factor(reservoir_type), \n",
    "          color = factor(reservoir_type))) +\n",
    "    scale_y_log10(limits = c(1,10000)) +\n",
    "    scale_x_discrete() +\n",
    "    scale_colour_manual(values = c(\"#d8b365\", \"#5ab4ac\", \"#2471a3\"), name=\"Reservoir Type\")\n",
    "# Add box plots visualising distributions across x and y axis\n",
    "p1 <- cowplot::insert_xaxis_grob(base, xbox, grid::unit(0.7, \"in\"), position = \"top\")\n",
    "p2 <- cowplot::insert_yaxis_grob(p1, ybox, grid::unit(0.7, \"in\"), position = \"right\")\n",
    "cowplot::ggdraw(p2)\n",
    "cowplot::save_plot('figures/maps/stats.png', p2, base_height = 6, ncol = 1, base_asp = 1.6)\n",
    "cowplot::save_plot('figures/maps/stats.svg', p2, base_height = 6, ncol = 1, base_asp = 1.6)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}