Functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6cc83484",
   "metadata": {},
   "source": [
    "# Biomass Energy Potential Mapping and Analysis Tool (BEPMAT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5dc3645",
   "metadata": {},
   "source": [
    "### This notebook contains all the functions used in the project along with the data. We request you to kindly go through the supporting text on the GitHub repository and the article, available as a preprint (Insert Zenodo doi) to get an idea of the objectives and the methodology."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "918ce575",
   "metadata": {},
   "source": [
    "## Loading all the CSVs containing the required raster files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "823c91bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing a few important libraries essential to the work.\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "# Importing the Geoprocessing libraries \n",
    "import rasterio\n",
    "from rasterio.mask import mask\n",
    "from rasterio.plot import show\n",
    "from rasterio.crs import CRS\n",
    "from rasterio.transform import from_origin\n",
    "from rasterio.transform import Affine\n",
    "from rasterio import transform\n",
    "from rasterio.enums import Resampling\n",
    "from rasterio.io import MemoryFile\n",
    "\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "\n",
    "# Importing required libraries for plotting interactive rasters \n",
    "from bokeh.plotting import figure\n",
    "from bokeh.plotting import show as bokeh_show\n",
    "from bokeh.models import LinearColorMapper, ColorBar, HoverTool, GeoJSONDataSource\n",
    "import bokeh.palettes as bp\n",
    "import bokeh.plotting as bpl\n",
    "import matplotlib.patches as mpatches\n",
    "import plotly\n",
    "import plotly.graph_objects as go\n",
    "import plotly.io as pio\n",
    "import seaborn as sns\n",
    "from bokeh.palettes import magma , Blues8 , plasma\n",
    "\n",
    "# Importing required libraries to obtain shapefiles\n",
    "import gadm\n",
    "from gadm import GADMDownloader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4bde9e07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uploading all the CSVs into pandas DataFrmes\n",
    "\n",
    "potential_yield = pd.read_csv(\"./dataset/potentialyield.csv\")\n",
    "harvested_area = pd.read_csv(\"./dataset/harvest_data_actual.csv\") #The year is integer type\n",
    "production_values= pd.read_csv(\"./dataset/productionvalues_actualyield.csv\")\n",
    "exclusion_areas= pd.read_csv(\"./dataset/Exclusionareas.csv\")\n",
    "tree_cover_share= pd.read_csv(\"./dataset/Treecover_share_GAEZ.csv\")\n",
    "aez_classification = pd.read_csv(\"./dataset/Classificationzones57.csv\")\n",
    "\n",
    "\n",
    "# The production_values defined here will be used for calculating crop residue from the year 2000 and 2010 \n",
    "# under different conditions. The harvested area will be used for future residue from cropland calculations\n",
    "# and the potential yield will be used for future residue from cropland as well as biomass potential from future\n",
    "# marginal land. The remaining will be used to find the Total Available Land."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c0f8311",
   "metadata": {},
   "source": [
    "## Creating a shapefile generator which can generate the gadm shapefile for any region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "baeaf2bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def shapefile_generator(country, province=None):\n",
    "    downloader = GADMDownloader(version=\"4.0\")\n",
    "    \n",
    "    if province:\n",
    "        # Download shapefile for a specific province\n",
    "        ad_level = 1\n",
    "        country_name = country\n",
    "        gdf_country = downloader.get_shape_data_by_country_name(country_name=country_name, ad_level=ad_level)\n",
    "        gdf_province = gdf_country[gdf_country[\"NAME_1\"] == province]\n",
    "        return gdf_province\n",
    "    else:\n",
    "        # Download shapefile for the entire country\n",
    "        ad_level = 0\n",
    "        gdf_country = downloader.get_shape_data_by_country_name(country_name=country, ad_level=ad_level)\n",
    "        return gdf_country"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63714350",
   "metadata": {},
   "source": [
    "## Notebook work flow: \n",
    "- We start with the simplest calculations of the raw biomass energy we could have obtained from the harvests in the past using the data for years 2000 and 2010 (This is only applicable for cropland since this was the land that was actually harvested and we are keeping the marginal for maximizing the energy output in the future).\n",
    "- Next we will calculate the residue and the raw biomass energy potential from the cropland in the future \n",
    "- Finally we will calculate the residue and the raw biomass energy potential from the marginal land in the future\n",
    "- We will keep defining helper functions along the way wherever required."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b2ae977",
   "metadata": {},
   "source": [
    "## I. Raw Biomass Energy Potential from Agricultural Residues using Actual Yields and Production (2000 and 2010) [GAEZv4 Theme 5 : Actual Yields and Productions]\n",
    "Theme 5 spatial layers include mapped distributions of harvested area, yield and production at 5 arc-minute resolution for 26 major crops/crop groups, separately in rain-fed and irrigated cropland. Country totals are based on FAO statistics for the years 2009-2011. Also included are estimates of the spatial distribution of total crop production value and the production values of major crop groups (cereals, root crops, oil crops), all valued at year 2000 international prices, separately for rain-fed and irrigated cropland. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "067900d2",
   "metadata": {},
   "source": [
    "### Creating the dataset for calculating the biomass energy potential from the production values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a371b63",
   "metadata": {},
   "source": [
    " To be able to calculate the Biomass potential we will need the crop Residue-to-Product Ratio, Surplus Availability Factor/\n",
    " Availability Factor and Lower Heating Value. The following table summarizes these values for the crops we have from the \n",
    " actual yields and production data:\n",
    "\n",
    "| Crop         | Residue Type  | RPR   | SAF   | LHV (MJ/kg) | Sources         |\n",
    "|--------------|---------------|-------|-------|-------------|-----------------|\n",
    "| Maize        | Stalk         | 2     | 0.8   | 16.3        | a, b, j         |\n",
    "|              | Cob           | 0.273 | 1     | 16.63       |                 |\n",
    "|              | Husk          | 0.2   | 1     | 15.56       |                 |\n",
    "| Rice         | Straw         | 1.757 | 0.684 | 8.83        | b, d, e         |\n",
    "|              | Husk          | 0.23  | 0.83  | 12.9        | c, e, f         |\n",
    "| Sorghum      | Straw         | 1.25  | 0.8   | 12.38       | a, b            |\n",
    "|              | Husk          | 1.4   | 1     | 13          | c, j            |\n",
    "| Millet       | Straw         | 1.4   | 1     | 13          | c, j            |\n",
    "|              | Stalk         | 1.75  | 0.8   | 15.51       | a, b, f         |\n",
    "| Wheat        | Straw         | 1.2   | 0.29  | 15.6        | b, j            |\n",
    "|              | Husk          | 0.23  | 0.29  | 12.9        | b, f            |\n",
    "| Cassava      | Stalk         | 0.062 | 0.407 | 16.99       | a, d, e         |\n",
    "|              | Peelings      | 3     | 0.2   | 10.61       | a, i            |\n",
    "| Cocoyam      | Peelings      | 0.2   | 0.8   | 10.61       | i, j            |\n",
    "| Sweet potato | Peelings      | 0.6   | 0.8   | 10.61       | b, j            |\n",
    "| Yam          | Peelings      | 0.2   | 0.8   | 10.61       | i, j            |\n",
    "| Potatoes     | Peelings      | 0.75  | 0.8   | 10.61       | i, j            |\n",
    "| Groundnuts   | Shells/husks  | 0.477 | 1     | 15.56       | a, i, c         |\n",
    "|              | Straw         | 2.3   | 1     | 17.58       | a               |\n",
    "| Palm oil     | Fiber         | 0.147 | 1     | 19.94       | a, i            |\n",
    "|              | Shells        | 0.049 | 1     | 21.1        | a, i            |\n",
    "|              | Fronds        | 2.604 | 1     | 7.97        | i               |\n",
    "|              | Empty bunches | 0.428 | 1     | 19.41       | a, i            |\n",
    "|              | Male bunches  | 0.233 | 1     | 14.86       | i, j            |\n",
    "| Beans(In Pulses)| Straw      | 2.5   | 1     | 12.38       | j               |\n",
    "| Soybean      | Straw         | 2.66  | 0.8   | 18          | b, f            |\n",
    "|              | Pods          | 1     | 0.8   | 18          | a, b, f         |\n",
    "| Banana       | leaves        | 0.35  | 1     | 11.37       | g               |\n",
    "|              | stem          | 5.6   | 1     | 11.66       | a, j            |\n",
    "|              | peels         | 0.25  | 1     | 17          | h, j            |\n",
    "| Plantain(With Bananas)| leaves| 0.35 | 0.8   | 12.12       | g, i            |\n",
    "|              | stem          | 3.91  | 0.8   | 10.9        | g, i            |\n",
    "|              | peels         | 0.25  | 1     | 12.56       | a, h            |\n",
    "| Sugar Cane   | baggase       | 0.25  | 1     | 6.43        | b, c            |\n",
    "|              | tops/leaves   | 0.32  | 0.8   | 15.8        | b, c            |\n",
    "| Coffee       | husk          | 1     | 1     | 12.8        | b, c            |\n",
    "| Cocoa        | pods/husks    | 1     | 1     | 15.48       | j               |\n",
    "| Cotton       | stalk         | 2.1   | 1     | 15.9        | c, i            |\n",
    "| Barley       | straw         | 0.75  | 0.15  | 17.5        | k               |\n",
    "|              | stalk         | 1.60  | 0.60  | 18.5        | k               |\n",
    "| Tobacco      | stalk         | 1.20  | 0.60  | 16.1        | k               |\n",
    "| Sunflower    | stalk         | 2.50  | 0.60  | 14.2        | k               |\n",
    "| Sugarbeet    | residue       | 0.66  | 0.09  | 20.85       | p,q             |\n",
    "| Rapeseed     | straw         | 1.58  | 0.23  | 14.55       | l               |    \n",
    "| Olives       | cake          | 0.40  | 0.90  | 19.7        | k               |\n",
    "| Lettuce      | waste         |1.2    | 0.50  | 12.8        | l               |\n",
    "| Tomatoes     | stem          |0.3    | 0.50  | 13.7        | l               |\n",
    "| Tomatoes     | leaves        |0.3    | 0.50  | 13.7        | l               |\n",
    "| Green peppers| residues      |0.45   | 0.50  | 12.0        | l               |\n",
    "| Red Peppers  | residues      |0.45   | 0.50  | 12.0        | l               |\n",
    "| Other Cereals| straw         | 1.2   | 0.40  | 16.845      | n,o             |\n",
    "| Rest of Crops| residue       | 0     | 0     | 0           | Not found       |\n",
    "| Fodder Crops | straw         | 0.4   | 0     | 0           | p               |   \n",
    "| Tur          | stalk         | 2.5   | 0.38  |18.58        | m               |\n",
    "| Lentils      | stalk         | 1.8   | 0.38  |14.65        | m               |\n",
    "| Gaur         | stalk         | 1.0   | 0.38  |16.02        | m               |\n",
    "| Gram         | stalk         | 1.1   | 0.38  |16.02        | m               |\n",
    "| Fruits and Nuts|pruning      | 0     | 0     |0            | Not found       |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ed1999",
   "metadata": {},
   "source": [
    "So for calculating the biomass potential of residues we have the production values data from GAEZ v4 for the years 2000 and 2010 which we will be using. The data in these rasters gives us the production of the particular crop in 1000 tonnes or 1 mln GK\\\\$. The other unit present in the data is mln GK\\\\$ which is used by FAO for crop groups like Fodder Crops, Pulses, Vegetables etc.\n",
    "\n",
    "The documentation of the GAEZ v4 describes the yield as either tonnes/hectare or 1000 GK\\\\$/hectare. From this, we derive that 1 million GK\\\\$ = 1000 tonnes.\n",
    "\n",
    "The following group contains the following crops according the GAEZ v4 Documentation:\n",
    "- Fodder Crops: All commodities in FAOSTAT primary crop production domain ranging from forage and silage, maize to vegetables and roots fodders.\n",
    "- Pulses: Bambara beans; beans, dry; broad beans, dry; chick peas; cow peas, dry;lentils; peas, dry; pigeon peas; pulses, other\n",
    "- Other cereals: Buckwheat; canary seed; fonio; mixed grain; oats; pop corn; quinoa; rye; triticale;\n",
    "- Yams and other roots: Taro; yautia; yams; roots and tubers;\n",
    "- Other crops: Includes all other crops from FAOSTAT production domain not covered by 25 crop groups above and excluding coir, vegetable tallow, oil of stillinga, oil of citronella, essential oils and rubber, natural.\n",
    "\n",
    "\n",
    "Assumptions made on the basis of GAEZ V4 documentation, from which we are using actual yield data:\n",
    "- Stimulants in the GAEZ v4 includes Cocoa Beans, Coffee, Green Tea, Tea. We have till now considered only Cocoa and Coffee.\n",
    "- In Yams and other roots, till now we have only included Yams and Coco Yams(Taro)\n",
    "- In Vegetables we have taken Green and Red Peppers (Residues), Tomatoes (Stem and Leaves) and Lettuce (Waste) (We will take their mean value for RPR, SAF, and LHV since individual weightage in the crop yield is not available.)\n",
    "- In Pulses, we have taken Tur, Gaur, Gram, Beans and Lentils.(We will use their average based on same reason as above).\n",
    "- In Other Cereals, we have taken oats and rye.(We will use their average based on same reason as above).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa310c16",
   "metadata": {},
   "source": [
    "## References for the above data:\n",
    "RPR and LHV values given were obtained from already published studies conducted in other countries, such as Ghana, Uganda, Zambia and China.\n",
    "\n",
    "A source for most of these references was: https://www.aimspress.com/article/doi/10.3934/energy.2023002?viewType=HTML.\n",
    "The references are as follows:\n",
    "- a. Jekayinfa SO, Scholz V (2009) Potential availability of energetically usable crop residues in Nigeria. Energy Sources, Part A: Recovery, Util, Environ Effects 31: 687–697. https://doi.org/10.1080/15567030701750549 doi: 10.1080/15567030701750549.\n",
    "- b. Gabisa EW, Gheewala SH (2018) Potential of bio-energy production in Ethiopia based on available biomass residues. Biomass Bioenergy 111: 77–87. https://doi.org/10.1016/j.biombioe.2018.02.009 doi: 10.1016/j.biombioe.2018.02.009.\n",
    "- c. Okello C, Pindozzi S, Faugno S, et al. (2013) Bioenergy potential of agricultural and forest residues in Uganda. Biomass Bioenergy 56: 515–525. https://doi.org/10.1016/j.biombioe.2013.06.003 doi: 10.1016/j.biombioe.2013.06.003.\n",
    "- d. Koopmans A, Koppenjan J (1998) The Resource Base. Reg Consult Mod Appl Biomass Energy, 6–10.\n",
    "- e. San V, Ly D, Check NI (2013) Assessment of sustainable energy potential on non-plantation biomass resources in Sameakki Meanchey district in Kampong Chhnan pronice, Cambonia. Int J Environ Rural Dev 4: 173–178.\n",
    "- f. Yang J, Wang X, Ma H, et al. (2014) Potential usage, vertical value chain and challenge of biomass resource: Evidence from China's crop residues. Appl Energy 114: 717–723. https://doi.org/10.1016/j.apenergy.2013.10.019 doi: 10.1016/j.apenergy.2013.10.019.\n",
    "- g. Patiño FGB, Araque JA, Kafarov DV (2016) Assessment of the energy potential of agricultural residues in non-interconnected zones of Colombia: Case study of Chocó and Putumayo katherine Rodríguez cáceres. Chem Eng Trans 50: 349–354. https://doi.org/10.3303/CET1650059 doi: 10.3303/CET1650059.\n",
    "- h. Milbrandt A (2011) Assessment of biomass resources in Liberia. Liberia: Dev Resour, 117–166.\n",
    "- i.Kemausuor F, Kamp A, Thomsen ST, et al. (2014) Assessment of biomass residue availability and bioenergy yields in Ghana. Resou Conser Recycl 86: 28–37. https://doi.org/10.1016/j.resconrec.2014.01.007 doi: 10.1016/j.resconrec.2014.01.007.\n",
    "- j. Mboumboue E, Njomo D (2018) Biomass resources assessment and bioenergy generation for a clean and sustainable development in Cameroon. Biomass Bioenergy 118: 16–23. https://doi.org/10.1016/j.biombioe.2018.08.002 doi: 10.1016/j.biombioe.2018.08.002.\n",
    "- k. https://www.researchgate.net/publication/342000532_Agricultural_Residues_Potential_of_Hatay.\n",
    "- l. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9941997/.\n",
    "- m. https://www.saarcenergy.org/wp-content/uploads/2020/03/Final-Draft-SEC-report-on-crop-residue_14022020-1540-IM-1.pdf.\n",
    "- n. https://www.sciencedirect.com/science/article/pii/S0956053X10002436?via%3Dihub.\n",
    "- o. https://www.sciencedirect.com/science/article/pii/S0921344920305280.\n",
    "- p. https://www.diva-portal.org/smash/get/diva2:1208954/FULLTEXT01.pdf.\n",
    "- q.https://www.researchgate.net/publication/317490809_VALORIZATION_OF_SUGAR_BEET_PULP_RESIDUE_AS_A_SOLID_FUEL_VIA_TORREFACTION"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2abf356a",
   "metadata": {},
   "source": [
    "Now the final table sorted by crop names which will be converted into a pandas data frame for us to use will contain: \n",
    "- All the vegetables combined into Vegetables row.\n",
    "- Coffee & Cocoa combined under stimulants row.\n",
    "- All the pulses as mentioned above will be grouped under pulses row.\n",
    "\n",
    "<h2><center>Final Crop Table with RPR, SAF and LHV values</center></h2>\n",
    "\n",
    "| Crop         | Residue Type  | RPR   | SAF   | LHV (MJ/kg) |\n",
    "|--------------|---------------|-------|-------|-------------|\n",
    "| Banana       | leaves        | 0.35  | 0.9   | 11.745      |\n",
    "| Banana       | peels         | 0.25  | 1     | 14.78       |\n",
    "| Banana       | stem          | 4.90  | 0.9   | 11.66       |\n",
    "| Barley       | stalk         | 1.60  | 0.60  | 18.5        |\n",
    "| Barley       | straw         | 0.75  | 0.15  | 17.5        |\n",
    "| Cassava      | Peelings      | 3     | 0.2   | 10.61       |\n",
    "| Cassava      | Stalk         | 0.062 | 0.407 | 16.99       |\n",
    "| Cotton       | stalk         | 2.1   | 1     | 15.9        |\n",
    "| Fodder Crops | straw         | 0.4   | 0     | 0           |\n",
    "|Fruits and nuts| pruning      | 0     | 0     | 0           |\n",
    "| Groundnut    | Shells/husks  | 0.477 | 1     | 15.56       |\n",
    "| Groundnut    | Straw         | 2.3   | 1     | 17.58       |\n",
    "| Maize        | Cob           | 0.273 | 1     | 16.63       |\n",
    "| Maize        | Husk          | 0.2   | 1     | 15.56       |\n",
    "| Maize        | Stalk         | 2     | 0.8   | 16.3        |\n",
    "| Millet       | Stalk         | 1.75  | 0.8   | 15.51       |\n",
    "| Millet       | Straw         | 1.4   | 1     | 13          |\n",
    "| Other Cereals| straw         | 1.2   | 0.40  | 16.845      |\n",
    "| Oil palm     | Empty bunches | 0.428 | 1     | 19.41       |\n",
    "| Oil palm     | Fiber         | 0.147 | 1     | 19.94       |\n",
    "| Oil palm     | Fronds        | 2.604 | 1     | 7.97        |\n",
    "| Oil palm     | Male bunches  | 0.233 | 1     | 14.86       |\n",
    "| Oil palm     | Shells        | 0.049 | 1     | 21.1        |\n",
    "| Olive        | Cake          | 0.40  | 0.9   | 19.7        |\n",
    "| Potato and Sweet Potato | Peelings | 0.675 |0.8  |10.61|\n",
    "| Pulses       | stalk         | 1.78  | 0.504 |15.53        |\n",
    "| Rapeseed     | straw         | 1.58  | 0.23  | 14.55       |\n",
    "| Wetland rice | Husk          | 0.23  | 0.83  | 12.9        |\n",
    "| Wetland rice | Straw         | 1.757 | 0.684 | 8.83        |\n",
    "| Sorghum      | Husk          | 1.4   | 1     | 13          |\n",
    "| Sorghum      | Straw         | 1.25  | 0.8   | 12.38       |\n",
    "| Soybean      | Pods          | 1     | 0.8   | 18          |\n",
    "| Soybean      | Straw         | 2.66  | 0.8   | 18          |\n",
    "| Stimulants   | husks         | 1     | 1     | 14.14       |\n",
    "| Sugar Cane   | bagasse       | 0.25  | 1     | 6.43        |\n",
    "| Sugar Cane   | tops/leaves   | 0.32  | 0.8   | 15.8        |\n",
    "| Sugar beet   | residue       | 0.66  | 0.09  | 20.85       |\n",
    "| Sunflower    | stalk         | 2.50  | 0.60  | 14.2        |\n",
    "| Tobacco      | stalk         | 1.20  | 0.60  | 16.1        |\n",
    "| Vegetables   | residue       | 0.675 | 0.50  | 12.625      |\n",
    "| Wheat        | Husk          | 0.23  | 0.29  | 12.9        |\n",
    "| Wheat        | Straw         | 1.2   | 0.29  | 15.6        |\n",
    "| Yam and others| Peelings     | 0.2   | 0.8   | 10.61       |\n",
    "| Rest of crops| Residue       | 0.0   | 0.0   | 0           |\n",
    "|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfbda11b",
   "metadata": {},
   "source": [
    "### REMARK : Incase you have your own RPR , SAF and LHV values for your region, we request you to fork this repository and modify the values accordingly to obtain more region specific results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2d25437a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now importing the table in pandas format so that we can use it for geospatial analysis\n",
    "# Defining the table data\n",
    "data = [\n",
    "    ['Banana', 'leaves', 0.35, 0.9, 11.745],\n",
    "    ['Banana', 'peels', 0.25, 1, 14.78],\n",
    "    ['Banana', 'stem', 4.90, 0.9, 11.66],\n",
    "    ['Barley', 'stalk', 1.60, 0.60, 18.5],\n",
    "    ['Barley', 'straw', 0.75, 0.15, 17.5],\n",
    "    ['Cassava', 'Peelings', 3, 0.2, 10.61],\n",
    "    ['Cassava', 'Stalk', 0.062, 0.407, 16.99],\n",
    "    ['Cotton', 'stalk', 2.1, 1, 15.9],\n",
    "    ['Fodder crops', 'straw', 0.4, 0, 0],\n",
    "    ['Fruits and nuts', 'Pruning' , 0, 0, 0],\n",
    "    ['Groundnut', 'Shells/husks', 0.477, 1, 15.56],\n",
    "    ['Groundnut', 'Straw', 2.3, 1, 17.58],\n",
    "    ['Maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Millet', 'Stalk', 1.75, 0.8, 15.51],\n",
    "    ['Millet', 'Straw', 1.4, 1, 13],\n",
    "    ['Other cereals', 'straw', 1.2, 0.40, 16.845],\n",
    "    ['Oil palm', 'Empty bunches', 0.428, 1, 19.41],\n",
    "    ['Oil palm', 'Fiber', 0.147, 1, 19.94],\n",
    "    ['Oil palm', 'Fronds', 2.604, 1, 7.97],\n",
    "    ['Oil palm', 'Male bunches', 0.233, 1, 14.86],\n",
    "    ['Oil palm', 'Shells', 0.049, 1, 21.1],\n",
    "    ['Olive' , 'Cake', 0.4, 0.9, 19.7],\n",
    "    ['Potato and Sweet Potato', 'Peelings', 0.675, 0.8, 10.61],\n",
    "    ['Pulses', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Rapeseed', 'straw', 1.58, 0.23, 14.55],\n",
    "    ['Wetland rice', 'Husk', 0.23, 0.83, 12.9],\n",
    "    ['Wetland rice', 'Straw', 1.757, 0.684, 8.83],\n",
    "    ['Sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Soybean', 'Pods', 1, 0.8, 18],\n",
    "    ['Soybean', 'Straw', 2.66, 0.8, 18],\n",
    "    ['Stimulants', 'husks', 1, 1, 14.14],\n",
    "    ['Sugarcane', 'baggase', 0.25, 1, 6.43],\n",
    "    ['Sugarcane', 'tops/leaves', 0.32, 0.8, 15.8],\n",
    "    ['Sugarbeet', 'residue', 0.66, 0.09, 20.85],\n",
    "    ['Sunflower', 'stalk', 2.50, 0.60, 14.2],\n",
    "    ['Tobacco', 'stalk', 1.20, 0.60, 16.1],\n",
    "    ['Vegetables', 'residue', 0.675, 0.50, 12.625],\n",
    "    ['Wheat', 'Husk', 0.23, 0.29, 12.9],\n",
    "    ['Wheat', 'Straw', 1.2, 0.29, 15.6],\n",
    "    ['Yams and other roots', 'Peelings', 0.2, 0.8, 10.61],\n",
    "    ['Rest of crops', 'Residue' , 0, 0, 0]\n",
    "]\n",
    "\n",
    "# Defining the column names\n",
    "columns = ['Crop', 'Residue Type', 'RPR', 'SAF', 'LHV (MJ/kg)']\n",
    "\n",
    "# Create the DataFrame\n",
    "residue_values = pd.DataFrame(data, columns=columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff87c56",
   "metadata": {},
   "source": [
    "### Defining the function for calculating raw biomass energy potential in the past (years 2000 and 2010)\n",
    "This function will output an xarray containing all the crops and their corresponding biomass energy potential in each pixel and a final xarray called 'Combined' which gives the sum of all of these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3dc8d412",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The band dimension comes up in a lot of places and is not needed for our calculations\n",
    "def remove_band_dimension(array):\n",
    "    if array.ndim == 3:\n",
    "        return array[0]\n",
    "    return array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "458c2d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_lat_lon_from_transform(transform, shape):\n",
    "    ny, nx = shape\n",
    "    lons , lats = transform * np.mgrid[:nx, :ny]\n",
    "    return lats, lons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3e38bfed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def biomass_potential_past(shapefile, time_period, water_supply):\n",
    "\n",
    "    unique_crops_actual = production_values['Crop'].unique()\n",
    "\n",
    "    filtered_production_values = production_values[(production_values['Time Period'] == time_period) &\n",
    "                                               (production_values['Water Supply'] == water_supply)]\n",
    "    required_production_values = filtered_production_values[['Crop', 'Download URL']]\n",
    "\n",
    "    # For defining size of the array to be used \n",
    "    with rasterio.open(potential_yield.iloc[2, 14].strip()) as src:\n",
    "        clipped_shapefile_init, transform_init = mask(src, shapefile.geometry, crop=True)\n",
    "        clipped_shapefile_init = remove_band_dimension(clipped_shapefile_init)\n",
    "        lats_init, lons_init = get_lat_lon_from_transform(transform_init, clipped_shapefile_init.shape)\n",
    "\n",
    "    # Create xarray Dataset to store individual biomass potentials for each crop\n",
    "    individual_biomass_potentials = {}\n",
    "\n",
    "    # Variable to hold the net sum of all crops and their residues\n",
    "    net_sum = 0.0\n",
    "\n",
    "    # Initialize 'net_biomass_potential_array' before the loop\n",
    "    net_biomass_potential_array = xr.DataArray(data=np.zeros_like(clipped_shapefile_init,  dtype='float32'),\n",
    "                                               dims=('y', 'x'),\n",
    "            coords={'y': range(clipped_shapefile_init.shape[0]), 'x': range(clipped_shapefile_init.shape[1]),\n",
    "                                'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "                                               attrs={'units': 'PetaJoules',\n",
    "                                                      'sum production': 0.0})  # Add initial sum as an attribute\n",
    "\n",
    "    for crop in unique_crops_actual:\n",
    "        required_url = required_production_values[required_production_values['Crop'] == crop]['Download URL'].values[0].strip()\n",
    "\n",
    "        with rasterio.open(required_url) as src:\n",
    "            crs_crop = src.crs\n",
    "            clipped_shapefile, clipped_transform = mask(src, shapefile.geometry, crop=True)\n",
    "            clipped_shapefile = remove_band_dimension(clipped_shapefile)\n",
    "            sum_value_shapefile = np.nansum(clipped_shapefile)\n",
    "\n",
    "        # Get the residues for the current crop\n",
    "        residue_rows = residue_values.loc[residue_values['Crop'] == crop]\n",
    "\n",
    "        # Calculate the sum of residues for the current crop\n",
    "        crop_residue_sum = 0.0\n",
    "        crop_residue_shapefile= np.zeros_like(clipped_shapefile)\n",
    "        \n",
    "        for _, residue_row in residue_rows.iterrows():\n",
    "            LHV = residue_row['LHV (MJ/kg)']\n",
    "            SAF = residue_row['SAF']\n",
    "            RPR = residue_row['RPR']\n",
    "            crop_residue_sum += sum_value_shapefile * LHV * SAF * RPR # Unit conversion for MJ to J and \n",
    "            # 1000 tonnes to kilograms ; Then multiplying by 10**-12 for PetaJoules\n",
    "            crop_residue_shapefile += clipped_shapefile * LHV * SAF * RPR \n",
    "\n",
    "    # Calculate the net sum for all crops and their residues\n",
    "        net_sum += crop_residue_sum\n",
    "\n",
    "    # Create xarray DataArray to store individual biomass potential for the current crop\n",
    "        crop_biomass_potential_array = xr.DataArray(data=crop_residue_shapefile,\n",
    "                                                    dims=('y', 'x'),\n",
    "            coords={'y': range(clipped_shapefile_init.shape[0]), 'x': range(clipped_shapefile_init.shape[1]),\n",
    "                                'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "                                                    attrs={'units': 'PetaJoules',\n",
    "                                                           'sum_production': crop_residue_sum})  # Add sum as an attribute\n",
    "\n",
    "        # Sum the individual biomass potential with the net biomass potential\n",
    "        net_biomass_potential_array += crop_biomass_potential_array\n",
    "\n",
    "        # Store the individual biomass potential for the current crop in the dictionary\n",
    "        individual_biomass_potentials[crop] = crop_biomass_potential_array\n",
    "\n",
    "    # Create xarray Dataset to hold all the individual biomass potentials for each crop\n",
    "    biomass_potentials_dataset = xr.Dataset(individual_biomass_potentials)\n",
    "\n",
    "    # Add the net sum of all crops and their residues as an attribute to the Dataset\n",
    "    biomass_potentials_dataset.attrs['Net Potential in PetaJ'] = net_sum\n",
    "\n",
    "    # Add the net_biomass_potential_array as a new variable named 'combined' to the biomass_potentials_dataset\n",
    "    biomass_potentials_dataset['Combined'] = net_biomass_potential_array\n",
    "    \n",
    "    biomass_potentials_dataset['Combined'].attrs['sum_production'] = net_sum\n",
    "\n",
    "    return biomass_potentials_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c81be42",
   "metadata": {},
   "source": [
    "#### Additional functions:\n",
    "- The following functions help estimate the final numbers for the raw biomass energy potential for a region. It has two options.\n",
    "    - Consolidated crop agnostic estimate\n",
    "    - Crop-specific estimate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3855cfad",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_actual_data_biomass_potential_all(shapefile, time_period, water_supply):\n",
    "    \n",
    "    value = biomass_potential_past(shapefile, time_period, water_supply)\n",
    "    \n",
    "    answer = value.attrs['net_sum']\n",
    "    \n",
    "    return answer\n",
    "\n",
    "def get_actual_data_biomass_potential_crop(shapefile, time_period, water_supply, crop):\n",
    "    \n",
    "    value = biomass_potential_past(shapefile, time_period, water_supply)\n",
    "    \n",
    "    answer = value[crop].attrs['sum_production']\n",
    "    \n",
    "    return answer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "859cb397",
   "metadata": {},
   "source": [
    "So the above functions and code finishes our task of getting the Raw Biomass Energy Potential from the Cropland in the past. Next we will see the functions for calculating the Raw Biomass Energy Potential from the Cropland in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44ff61ae",
   "metadata": {},
   "source": [
    "## II. Raw Biomass Energy Potential from Agricultural Residues using Actual Yields and Production for Harvested Area and Agro-Climatic Potential Yield for future yields [GAEZ-V4 Theme 5 and 3 respectively]\n",
    "\n",
    "Theme 3 provides crop-wise information about: (1) Agro-Climatic Yield, (2) Constraint Factors, (3) Growth Cycle Attributes, and (4) Land Utilization Types (LUT) Selection.\n",
    "\n",
    "Theme 5 spatial layers include mapped distributions of harvested area, yield and production at 5 arc-minute resolution for 26 major crops/crop groups, separately in rain-fed and irrigated cropland. Country totals are based on FAO statistics for the years 2009-2011. Also included are estimates of the spatial distribution of total crop production value and the production values of major crop groups (cereals, root crops, oil crops), all valued at year 2000 international prices, separately for rain-fed and irrigated cropland. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f089db83",
   "metadata": {},
   "source": [
    "Now since the future cropland data area data is not available to us we will be making an assumption. The assumption is that in the future the area under cropland which is required for us to ensure food security will remain the same as is was in the year 2010. So the harvest area data that we had from actual yields and production will serve as the cropland area for all future calculations. But since, with time the yield will vary and so will the residue from each crop. \n",
    "\n",
    "Assuming that RPR, SAF and LHV values also remain same for the crops in the future we will get the harvested area from the 2010 data for any shapefile and then multiply it with the future yield to get the future production. This will further be multiplied by RPR, SAF and LHV giving us raw biomass energy potential from cropland in the future. These values can then be compared with the past values calculated in part I to give us an idea as to how different future conditions affect the energy poential of the cropland. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a7810985",
   "metadata": {},
   "outputs": [],
   "source": [
    "def future_potential_cropland(time_period, climate_model, rcp, water_supply_future, input_level, shapefile_path, water_supply_2010):\n",
    "    \n",
    "    merged_df = pd.merge(harvested_area, potential_yield, on='Crop', how='inner')\n",
    "    unique_crops = merged_df['Crop'].unique()\n",
    "    \n",
    "    filtered_harvested_area = harvested_area[(harvested_area['Time Period'] == 2010) &\n",
    "                                             (harvested_area['Water Supply'] == water_supply_2010)]\n",
    "    required_harvested_area = filtered_harvested_area[['Crop', 'Download URL']]\n",
    "\n",
    "    filtered_potential_yield = potential_yield[(potential_yield['Time Period'] == time_period) &\n",
    "                                               (potential_yield['Climate Model'] == climate_model) &\n",
    "                                               (potential_yield['RCP'] == rcp) &\n",
    "                                               (potential_yield['Water Supply'] == water_supply_future) &\n",
    "                                               (potential_yield['Input Level'] == input_level)]\n",
    "    required_potential_yields = filtered_potential_yield[['Crop', 'Download URL']]\n",
    "    \n",
    "    # For defining size of the xarray\n",
    "    with rasterio.open(potential_yield.iloc[2, 14].strip()) as src:\n",
    "        clipped_shapefile_init, transform_init = mask(src, shapefile_path.geometry, crop=True)\n",
    "        clipped_shapefile_init = remove_band_dimension(clipped_shapefile_init)\n",
    "        lats_init, lons_init = get_lat_lon_from_transform(transform_init, clipped_shapefile_init.shape)\n",
    "\n",
    "\n",
    "    # Create xarray DataArray to store the net biomass potential for each pixel\n",
    "    net_biomass_potential_array = xr.DataArray(data=0.0,\n",
    "                                               dims=('y', 'x'),\n",
    "            coords={'y': range(clipped_shapefile_init.shape[0]), 'x': range(clipped_shapefile_init.shape[1]),\n",
    "                                'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "                                               attrs={'units': 'PetaJoules'})\n",
    "    \n",
    "\n",
    "    # Create xarray Dataset to store individual biomass potentials for each crop\n",
    "    individual_biomass_potentials = {}\n",
    "\n",
    "    # Variable to store the net sum of sum products\n",
    "    net_sum = 0.0\n",
    "    net_multiple = 0\n",
    "\n",
    "    for crop in unique_crops:\n",
    "        harvested_raster_url = required_harvested_area[required_harvested_area['Crop'] == crop]['Download URL'].values[0].strip()\n",
    "        potential_yield_raster_url = required_potential_yields[required_potential_yields['Crop'] == crop]['Download URL'].values[0].strip()\n",
    "\n",
    "        with rasterio.open(harvested_raster_url) as src:\n",
    "            clipped, _ = mask(src, shapefile_path.geometry, crop=True)\n",
    "\n",
    "        with rasterio.open(potential_yield_raster_url.strip()) as src:\n",
    "            clipped_2, _ = mask(src, shapefile_path.geometry, crop=True)\n",
    "        \n",
    "        clipped_1 = np.nan_to_num(clipped)\n",
    "        clipped_3 = np.nan_to_num(clipped_2)\n",
    "        \n",
    "         # Remove the 'bands' dimension if it exists (since it's not needed)\n",
    "        clipped_1 = remove_band_dimension(clipped_1)\n",
    "        clipped_3 = remove_band_dimension(clipped_3)\n",
    "\n",
    "        product = np.multiply(clipped_1, clipped_3)\n",
    "        sum_product = np.nansum(product)\n",
    "\n",
    "        # Extract RPR, SAF, and LHV values for the crop\n",
    "        residue_rows = residue_values[residue_values['Crop'] == crop]\n",
    "        temp_sum = 0\n",
    "        net_product_array = np.zeros_like(product)\n",
    "\n",
    "        for _, residue_row in residue_rows.iterrows():\n",
    "            LHV = residue_row['LHV (MJ/kg)']\n",
    "            SAF = residue_row['SAF']\n",
    "            RPR = residue_row['RPR']\n",
    "\n",
    "            # Multiply the sum_product with RPR, SAF, and LHV\n",
    "            result = sum_product * RPR * SAF * LHV * (10 ** -3) # Unit conversion factor to PetaJoules\n",
    "            temp_sum += result\n",
    "            product_array = product * RPR * SAF * LHV * (10 ** -3)\n",
    "            net_product_array += product_array\n",
    "            \n",
    "            \n",
    "\n",
    "        net_sum += temp_sum\n",
    "        \n",
    "\n",
    "        # Create xarray DataArray to store individual biomass potential for the current crop\n",
    "        crop_biomass_potential_array = xr.DataArray(data= net_product_array,\n",
    "                                                    dims=('y', 'x'),\n",
    "            coords={'y': range(clipped_shapefile_init.shape[0]), 'x': range(clipped_shapefile_init.shape[1]),\n",
    "                                'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "                                                    attrs={'units': 'PetaJoules',\n",
    "                                                           'sum_production': temp_sum})  # Add sum as an attribute\n",
    "\n",
    "        # Sum the individual biomass potential with the net biomass potential\n",
    "        net_biomass_potential_array += crop_biomass_potential_array\n",
    "\n",
    "        # Store the individual biomass potential for the current crop in the dictionary\n",
    "        individual_biomass_potentials[crop] = crop_biomass_potential_array\n",
    "\n",
    "    # Create xarray Dataset to hold all the individual biomass potentials for each crop\n",
    "    biomass_potentials_dataset = xr.Dataset(individual_biomass_potentials)\n",
    "\n",
    "    # Add the net sum of all crops and their residues as an attribute to the Dataset\n",
    "    biomass_potentials_dataset.attrs['net_sum in PJ'] = net_sum\n",
    "\n",
    "    # Add the net_biomass_potential_array as a new variable named 'combined' to the biomass_potentials_dataset\n",
    "    biomass_potentials_dataset['Combined'] = net_biomass_potential_array\n",
    "    \n",
    "    biomass_potentials_dataset['Combined'].attrs['sum_production'] = net_sum\n",
    "\n",
    "    return biomass_potentials_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60e0e0a2",
   "metadata": {},
   "source": [
    "#### Additional functions:\n",
    "The following functions are available if you just need the final numbers for the biomass energy potential for the \n",
    "region. It has two options either it can give you the net or it can give you the values for a specific crop as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "611219c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function doing as described above\n",
    "\n",
    "def future_residues_all(time_period, climate_model, rcp, water_supply_future, input_level, shapefile_path, water_supply_2010):\n",
    "    \n",
    "    value = future_potential_cropland(time_period, climate_model, rcp, water_supply_future, input_level, shapefile_path, water_supply_2010)\n",
    "    \n",
    "    answer = value.attrs['net_sum']\n",
    "    \n",
    "    return answer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ba050404",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We also wanted to create a function that does this for a single crop as well.\n",
    "\n",
    "def future_residues_crop(crop, time_period, climate_model, rcp, water_supply_future, input_level, shapefile_path, water_supply_2010):\n",
    "    \n",
    "    value = future_potential_cropland(time_period, climate_model, rcp, water_supply_future, input_level, shapefile_path, water_supply_2010)\n",
    "    \n",
    "    answer = value[crop].attrs['sum_production']\n",
    "    \n",
    "    return answer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7e150c6",
   "metadata": {},
   "source": [
    "So this completes all cropland residue calculations for us. Next we will move to residue and biomass energy potential from the marginal land but before we get into this we will have to generate the marginal land. So we need to extract certain pixels from certain rasters which will be masked later to account for removal of deserts, water, glaciers etc. as described in the paper"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f24c0bf9",
   "metadata": {},
   "source": [
    "## III. Raw Biomass Energy Potential from Agricultural Residues using Agro-Climatic Potential Yield for future yields data [GAEZ- V4 Theme: 3]\n",
    "\n",
    "Theme 3 provides crop-wise information about: (1) Agro-Climatic Yield, (2) Constraint Factors, (3) Growth Cycle Attributes, and (4) Land Utilization Types (LUT) Selection."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3428266",
   "metadata": {},
   "source": [
    "Before moving on to describe the functions, we need to obtain the RPR, SAF and LHV values for the crops available in the potential yield theme in GAEZ. The following is the table followed by the assumptions made and references:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d440e3",
   "metadata": {},
   "source": [
    "<h2><center>Final Potential Yield Crop Table with RPR, SAF and LHV values</center></h2>\n",
    "\n",
    "| Crop                  | Residue Type    | RPR   | SAF   | LHV (MJ/kg) |\n",
    "|:----------------------:|:-----------------:|:-------:|:-------:|:-------------:|\n",
    "| Alfalfa               | residue         | 0.25  | 0.0   | 0.0         |\n",
    "| Banana                | leaves          | 0.35  | 0.9   | 11.745      |\n",
    "| Banana                | peels           | 0.25  | 1     | 14.78       |\n",
    "| Banana                | stem            | 4.90  | 0.9   | 11.66       |\n",
    "| Barley                | stalk           | 1.60  | 0.60  | 18.5        |\n",
    "| Barley                | straw           | 0.75  | 0.15  | 17.5        |\n",
    "| Biomass highland sorghum |Husk            | 1.4   | 1     | 13          |\n",
    "| Biomass highland sorghum | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Biomass lowland sorghum  |Husk            | 1.4   | 1     | 13          |\n",
    "|Biomass lowland sorghum  | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Biomass sorghum       |Husk            | 1.4   | 1     | 13          |\n",
    "|Biomass sorghum        | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Biomass temperate sorghum | Husk            | 1.4   | 1     | 13          |\n",
    "| Biomass temperate sorghum | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Buckwheat             | straw         | 1.2   | 0.40  | 16.845      |\n",
    "| Cabbage               | residue       | 0.675 | 0.50  | 12.625      |\n",
    "| Carrot                | residue       | 0.675 | 0.50  | 12.625      |\n",
    "| Cassava               | Peelings        | 3     | 0.2   | 10.61       |\n",
    "| Cassava               | Stalk           | 0.062 | 0.407 | 16.99       |\n",
    "| Chickpea              | stalk         | 1.78  | 0.504 |15.53        |\n",
    "| Citrus                | pruning        | 0.29     | 0.8   | 17.85       |\n",
    "| Cocoa                 |  pods/husks    | 1     | 1     | 15.48       |\n",
    "| Cocoa cumoun          |  pods/husks    | 1     | 1     | 15.48       |\n",
    "| Cocoa hybrid          |  pods/husks    | 1     | 1     | 15.48       |\n",
    "| Coconut               | Husk           | 1.03  | 1     | 18.6        |\n",
    "| Coconut               | Coir dust      | 0.62  | 1     | 13.4        |\n",
    "| Cocoyam               |  Peelings      | 0.2   | 0.8   | 10.61       |\n",
    "| Coffee                | husk          | 1     | 1     | 12.8         |\n",
    "| Coffee arabica        |husk          | 1     | 1     | 12.8          |\n",
    "| Coffee robusta        |husk          | 1     | 1     | 12.8          |\n",
    "| Cotton                | stalk           | 2.1   | 1     | 15.9       |\n",
    "| Cowpea                | stalk         | 1.78  | 0.504 |15.53         |\n",
    "| Dry pea               |stalk         | 1.78  | 0.504 |15.53          |\n",
    "| Dryland rice          | Husk            | 0.23  | 0.83  | 12.9       |\n",
    "| Dryland rice          | Straw           | 1.757 | 0.684 | 8.83       |\n",
    "| Flax                  |  stalk           | 2.50  | 0.60  | 14.2      |\n",
    "| Foxtail millet        |  Stalk           | 1.75  | 0.8   | 15.51     |\n",
    "| Foxtail millet                | Straw           | 1.4   | 1     | 13 |\n",
    "| Gram                  |stalk         | 1.1   | 0.38  |16.02          |\n",
    "| Grass                 |  straw         | 0.4   | 0     | 0           |\n",
    "| Greater yam           |Peelings      | 0.2   | 0.8   | 10.61         |\n",
    "| Groundnut             | Shells/husks    | 0.477 | 1     | 15.56      |\n",
    "| Groundnut             | Straw           | 2.3   | 1     | 17.58      |\n",
    "| Highland maize        | Cob             | 0.273 | 1     | 16.63      |\n",
    "| Highland maize                 | Husk            | 0.2   | 1     | 15.56|\n",
    "| Highland maize                 | Stalk           | 2     | 0.8   | 16.3 |\n",
    "| Highland sorghum      | Husk            | 1.4   | 1     | 13         |\n",
    "|Highland sorghum       | Straw           | 1.25  | 0.8   | 12.38      |\n",
    "| Jatropha | Woody Biomass | 0.25 | 0.8 | 15.5 |\n",
    "| Jatropha | Leaves | 0.25 | 0.8 | 12 |\n",
    "| Jatropha | Pericarp | 0.20 | 0.8 | 10 |\n",
    "| Jatropha | Tegument | 0 | 10.8 | 16.9 |\n",
    "| Jatropha | Endosperm Cake | 0.2 | 0.8 | 13.6 |\n",
    "| Lowland maize         | Cob             | 0.273 | 1     | 16.63       |\n",
    "| Lowland maize                 | Husk            | 0.2   | 1     | 15.56       |\n",
    "| Lowland maize                 | Stalk           | 2     | 0.8   | 16.3        |\n",
    "| Lowland sorghum       |Husk            | 1.4   | 1     | 13          |\n",
    "|Lowland sorghum        | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Maize                 | Cob             | 0.273 | 1     | 16.63       |\n",
    "| Maize                 | Husk            | 0.2   | 1     | 15.56       |\n",
    "| Maize                 | Stalk           | 2     | 0.8   | 16.3        |\n",
    "| Millet                | Stalk           | 1.75  | 0.8   | 15.51       |\n",
    "| Millet                | Straw           | 1.4   | 1     | 13          |\n",
    "| Miscanthus            | residue        | 0.42   | 1     | 17.44       |\n",
    "| Napier grass          |  straw         | 0.4   | 0     | 0            |\n",
    "| Oat                   | straw  |    1.15 |  0.40 | 18.45|\n",
    "| Oil palm              | Empty bunches   | 0.428 | 1     | 19.41       |\n",
    "| Oil palm              | Fiber           | 0.147 | 1     | 19.94       |\n",
    "| Oil palm              | Fronds          | 2.604 | 1     | 7.97        |\n",
    "| Oil palm              | Male bunches    | 0.233 | 1     | 14.86       |\n",
    "| Oil palm              | Shells          | 0.049 | 1     | 21.1        |\n",
    "| Olive                 | cake          | 0.40  | 0.90  | 19.7        |\n",
    "| Onion                 | residue       | 0.675 | 0.50  | 12.625      |\n",
    "| Para rubber           | residue               | 0     | 0     | 0           |\n",
    "| Pasture legumes       | stalk         | 1.78  | 0.504 |15.53        |\n",
    "| Pearl millet          |  Stalk           | 1.75  | 0.8   | 15.51       |\n",
    "| Pearl millet                | Straw           | 1.4   | 1     | 13          |\n",
    "| Phaseolus bean        |stalk         | 1.78  | 0.504 |15.53        |\n",
    "| Pigeonpea             |stalk         | 1.78  | 0.504 |15.53        |\n",
    "| Rapeseed              | straw           | 1.58  | 0.23  | 14.55       |\n",
    "| Reed canary grass     |  straw         | 0.4   | 0     | 0           |\n",
    "| Rye                   |   straw| 1.25|   0.40  | 15.24|\n",
    "| Silage maize          | Cob             | 0.273 | 1     | 16.63       |\n",
    "| Silage maize                 | Husk            | 0.2   | 1     | 15.56       |\n",
    "| Silage maize                 | Stalk           | 2     | 0.8   | 16.3        |\n",
    "| Sorghum               | Husk            | 1.4   | 1     | 13          |\n",
    "| Sorghum               | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Soybean               | Pods            | 1     | 0.8   | 18          |\n",
    "| Soybean               | Straw           | 2.66  | 0.8   | 18          |\n",
    "| Spring barley         | stalk           | 1.60  | 0.60  | 18.5        |\n",
    "| Spring barley                | straw           | 0.75  | 0.15  | 17.5        |\n",
    "| Spring rye            |  straw| 1.25|   0.40  | 15.24|\n",
    "| Spring wheat          |  Husk            | 0.23  | 0.29  | 12.9        |\n",
    "| Spring wheat                 | Straw           | 1.2   | 0.29  | 15.6        |\n",
    "| Sugarbeet             | residue         | 0.66  | 0.09  | 0           |\n",
    "| Sugar Cane   | baggase       | 0.25  | 1     | 6.43        |\n",
    "| Sugar Cane   | tops/leaves   | 0.32  | 0.8   | 15.8        |\n",
    "| Sunflower             | stalk           | 2.50  | 0.60  | 14.2        |\n",
    "| Sweet potato          |  Peelings      | 0.6   | 0.8   | 10.61       |\n",
    "| Switchgrass           | straw         | 0.4   | 0     | 0           |\n",
    "| Tea                   | husks         | 1     | 1     | 14.14       |\n",
    "| Temperate maize       | Cob             | 0.273 | 1     | 16.63       |\n",
    "| Temperate maize                 | Husk            | 0.2   | 1     | 15.56       |\n",
    "| Temperate maize                 | Stalk           | 2     | 0.8   | 16.3        |\n",
    "| Temperate sorghum     |Husk            | 1.4   | 1     | 13          |\n",
    "| Temperate sorghum     | Straw           | 1.25  | 0.8   | 12.38       |\n",
    "| Tobacco               | stalk           | 1.20  | 0.60  | 16.1        |\n",
    "| Tomato                | stem          |0.3    | 0.50  | 13.7        | l               |\n",
    "| Tomato                | leaves        |0.3    | 0.50  | 13.7        | l               |\n",
    "| Wetland rice          | Husk            | 0.23  | 0.83  | 12.9        |\n",
    "| Wetland rice          | Straw           | 1.757 | 0.684 | 8.83        |\n",
    "| Wheat                 | Husk            | 0.23  | 0.29  | 12.9        |\n",
    "| Wheat                 | Straw           | 1.2   | 0.29  | 15.6        |\n",
    "| White potato          | Peelings      | 0.75  | 0.8   | 10.61       |\n",
    "| White yam             | Peelings      | 0.2   | 0.8   | 10.61       |\n",
    "| Winter barley         |  stalk           | 1.60  | 0.60  | 18.5        |\n",
    "| Winter barley                | straw           | 0.75  | 0.15  | 17.5        |\n",
    "| Winter rye            |   straw| 1.25|   0.40  | 15.24|\n",
    "| Winter wheat          |  Husk            | 0.23  | 0.29  | 12.9        |\n",
    "| Winter wheat                 | Straw           | 1.2   | 0.29  | 15.6        |\n",
    "| Yam                   |Peelings      | 0.2   | 0.8   | 10.61       |\n",
    "| Yellow yam            | Peelings      | 0.2   | 0.8   | 10.61       |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36455f0d",
   "metadata": {},
   "source": [
    "Since some extra crops have been added to this table in comparison to the previous table here are the references and few assumptions that have been made in order to make this table complete:\n",
    "\n",
    "For a few of the crops I have assigned them value based on the following assumptions:\n",
    "\n",
    "- Biomass highland sorghum, Biomass lowland sorghum, Biomass temperate sorghum, Highland sorghum, Lowland sorghum, Temperate sorghum will all be assigned just the sorghum value.\n",
    "- Buckwheat will take the value of other cereals\n",
    "- Onion, Cabbage, Carrot will take the value for vegetables\n",
    "- Cocoa, Cocoa cumoun, Cocoa hybrid will take the value of Cocoa.\n",
    "- Coffee, Coffee arabica, Coffee robusta will take the value of Coffee.\n",
    "- Dryland Rice will take the value of Rice/Wetland Rice.\n",
    "- Pearl millet, Foxtail millet, Millet will take the value of Millet.\n",
    "- Greater yam, White yam, Yam, Yellow yam will all be assigned the value of Yam.\n",
    "- Highland maize, Lowland maize, Silage maize, Temperate maize and Maize will be assigned the value of Maize.\n",
    "- Spring barley, Winter barley will take the value of Barley \n",
    "- Spring rye, Winter rye will take the value of rye\n",
    "- Spring wheat, Winter wheat will take the value of wheat\n",
    "- White potato will take the value of Potato\n",
    "- Tea will take the value of stimulants.\n",
    "- ChickPea, Cowpea, Dry Pea, Pasture Legumes, Phaseolus bean, Pigeonpea will take the value of pulses.\n",
    "- Alfalfa, Grass, Napier grass, Para, Rubber, Reed canary grass will take the value of fodder crops.(Alfalfa takes green fodder values which will have different RPR(0.25) but SAF will be 0.)\n",
    "- Flax, being similar to sunflower, takes the values of sunflower since they belong to similar crop category.\n",
    "- Citrus here includes the average of Oranges and Lemons\n",
    "\n",
    "For the crops included above and not in the previous table here are the references:\n",
    "\n",
    "- Coconut : Reference: https://www.aimspress.com/article/doi/10.3934/energy.2023002?viewType=HTML\n",
    "- Citrus :  Reference:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9941997/\n",
    "- Jatropha : Reference : doi: 10.1016/j.rser.2015.10.009.\n",
    "- Miscanthus RPR/SAF & LHV: https://www.sciencedirect.com/science/article/pii/S1161030101001022. & https://www.researchgate.net/publication/338950136_Calorific_values_of_Miscanthus_x_giganteus_biomass_cultivated_under_suboptimal_conditions_in_marginal_soils            \n",
    "- Para Rubber : None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d155b5b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The final crop data being added to a pandas dataframe:\n",
    "\n",
    "data = [\n",
    "    ['Alfalfa', 'residue', 0.25, 0.0, 0.0],\n",
    "    ['Banana', 'leaves', 0.35, 0.9, 11.745],\n",
    "    ['Banana', 'peels', 0.25, 1, 14.78],\n",
    "    ['Banana', 'stem', 4.90, 0.9, 11.66],\n",
    "    ['Barley', 'stalk', 1.60, 0.60, 18.5],\n",
    "    ['Barley', 'straw', 0.75, 0.15, 17.5],\n",
    "    ['Biomass highland sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Biomass highland sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Biomass lowland sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Biomass lowland sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Biomass sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Biomass sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Biomass temperate sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Biomass temperate sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Buckwheat', 'straw', 1.2, 0.40, 16.845],\n",
    "    ['Cabbage', 'residue', 0.675, 0.50, 12.625],\n",
    "    ['Carrot', 'residue', 0.675, 0.50, 12.625],\n",
    "    ['Cassava', 'Peelings', 3, 0.2, 10.61],\n",
    "    ['Cassava', 'Stalk', 0.062, 0.407, 16.99],\n",
    "    ['Chickpea', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Citrus', 'prunings', 0.29 , 0.80 , 17.85],\n",
    "    ['Cocoa', 'pods/husks', 1, 1, 15.48],\n",
    "    ['Cocoa cumoun', 'pods/husks', 1, 1, 15.48],\n",
    "    ['Cocoa hybrid', 'pods/husks', 1, 1, 15.48],\n",
    "    ['Coconut', 'husk', 1.03, 1, 18.6],\n",
    "    ['Coconut','coir dust', 0.62, 1, 13.4],\n",
    "    ['Cocoyam', 'Peelings', 0.2, 0.8, 10.61],\n",
    "    ['Coffee', 'husk', 1, 1, 12.8],\n",
    "    ['Coffee arabica', 'husk', 1, 1, 12.8],\n",
    "    ['Coffee robusta', 'husk', 1, 1, 12.8],\n",
    "    ['Cotton', 'stalk', 2.1, 1, 15.9],\n",
    "    ['Cowpea', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Dry pea', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Dryland rice', 'Husk', 0.23, 0.83, 12.9],\n",
    "    ['Dryland rice', 'Straw', 1.757, 0.684, 8.83],\n",
    "    ['Flax', 'stalk', 2.5, 0.6, 14.2],\n",
    "    ['Foxtail millet', 'Stalk', 1.75, 0.8, 15.51],\n",
    "    ['Foxtail millet', 'Straw', 1.4, 1, 13],\n",
    "    ['Gram', 'stalk', 1.1, 0.38, 16.02],\n",
    "    ['Grass', 'straw', 0.4, 0, 0],\n",
    "    ['Greater yam', 'Peelings', 0.2, 0.8, 10.61],\n",
    "    ['Groundnut', 'Shells/husks', 0.477, 1, 15.56],\n",
    "    ['Groundnut', 'Straw', 2.3, 1, 17.58],\n",
    "    ['Highland maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Highland maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Highland maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Highland sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Highland sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Jatropha', 'Woody Biomass', 0.25, 0.8, 15.5],\n",
    "    ['Jatropha', 'Leaves', 0.25, 0.8, 12],\n",
    "    ['Jatropha', 'Pericarp', 0.20, 0.8, 10],\n",
    "    ['Jatropha', 'Tegument', 0, 10.8, 16.9],\n",
    "    ['Jatropha', 'Endosperm Cake', 0.2, 0.8, 13.6],\n",
    "    ['Lowland maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Lowland maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Lowland maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Lowland sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Lowland sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Millet', 'Stalk', 1.75, 0.8, 15.51],\n",
    "    ['Millet', 'Straw', 1.4, 1, 13],\n",
    "    ['Miscanthus', 'residue', 0.42, 1, 17.44],\n",
    "    ['Napier grass', 'straw', 0.4, 0, 0],\n",
    "    ['Oat', 'straw', 1.15, 0.4, 18.45],\n",
    "    ['Oil palm', 'Empty bunches', 0.428, 1, 19.41],\n",
    "    ['Oil palm', 'Fiber', 0.147, 1, 19.94],\n",
    "    ['Oil palm', 'Fronds', 2.604, 1, 7.97],\n",
    "    ['Oil palm', 'Male bunches', 0.233, 1, 14.86],\n",
    "    ['Oil palm', 'Shells', 0.049, 1, 21.1],\n",
    "    ['Olive', 'cake', 0.4, 0.9, 19.7],\n",
    "    ['Onion', 'residue', 0.675, 0.5, 12.625],\n",
    "    ['Para rubber', 'residue', 0, 0, 0],\n",
    "    ['Pasture legumes', 'stalk', 0, 0, 0],\n",
    "    ['Pearl millet', 'Stalk', 1.75, 0.8, 15.51],\n",
    "    ['Pearl millet', 'Straw', 1.4, 1, 13],\n",
    "    ['Phaseolus bean', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Pigeonpea', 'stalk', 1.78, 0.504, 15.53],\n",
    "    ['Rapeseed', 'straw', 1.58, 0.23, 14.55],\n",
    "    ['Reed canary grass', 'straw', 0.4, 0, 0],\n",
    "    ['Rye', 'straw', 1.25, 0.4, 15.24],\n",
    "    ['Silage maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Silage maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Silage maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Soybean', 'Pods', 1, 0.8, 18],\n",
    "    ['Soybean', 'Straw', 2.66, 0.8, 18],\n",
    "    ['Spring barley', 'stalk', 1.6, 0.6, 18.5],\n",
    "    ['Spring barley', 'straw', 0.75, 0.15, 17.5],\n",
    "    ['Spring rye', 'straw', 1.25, 0.4, 15.24],\n",
    "    ['Spring wheat', 'Husk', 0.23, 0.29, 12.9],\n",
    "    ['Spring wheat', 'Straw', 1.2, 0.29, 15.6],\n",
    "    ['Sugarbeet', 'residue', 0.66, 0.09, 0],\n",
    "    ['Sugar Cane', 'baggase', 0.25, 1, 6.43],\n",
    "    ['Sugar Cane', 'tops/leaves', 0.32, 0.8, 15.8],\n",
    "    ['Sunflower', 'stalk', 2.5, 0.6, 14.2],\n",
    "    ['Sweet potato', 'Peelings', 0.6, 0.8, 10.61],\n",
    "    ['Switchgrass', 'straw', 0.4, 0, 0],\n",
    "    ['Tea', 'husks', 1, 1, 14.14],\n",
    "    ['Temperate maize', 'Cob', 0.273, 1, 16.63],\n",
    "    ['Temperate maize', 'Husk', 0.2, 1, 15.56],\n",
    "    ['Temperate maize', 'Stalk', 2, 0.8, 16.3],\n",
    "    ['Temperate sorghum', 'Husk', 1.4, 1, 13],\n",
    "    ['Temperate sorghum', 'Straw', 1.25, 0.8, 12.38],\n",
    "    ['Tobacco', 'stalk', 1.2, 0.6, 16.1],\n",
    "    ['Tomato', 'stem', 0.3, 0.5, 13.7],\n",
    "    ['Tomato', 'leaves', 0.3, 0.5, 13.7],\n",
    "    ['Wetland rice', 'Husk', 0.23, 0.83, 12.9],\n",
    "    ['Wetland rice', 'Straw', 1.757, 0.684, 8.83],\n",
    "    ['Wheat', 'Husk', 0.23, 0.29, 12.9],\n",
    "    ['Wheat', 'Straw', 1.2, 0.29, 15.6],\n",
    "    ['White potato', 'Peelings', 0.75, 0.8, 10.61],\n",
    "    ['White yam', 'Peelings', 0.2, 0.8, 10.61],\n",
    "    ['Winter barley', 'stalk', 1.6, 0.6, 18.5],\n",
    "    ['Winter barley', 'straw', 0.75, 0.15, 17.5],\n",
    "    ['Winter rye', 'straw', 1.25, 0.4, 15.24],\n",
    "    ['Winter wheat', 'Husk', 0.23, 0.29, 12.9],\n",
    "    ['Winter wheat', 'Straw', 1.2, 0.29, 15.6],\n",
    "    ['Yam', 'Peelings', 0.2, 0.8, 10.61],\n",
    "    ['Yellow yam', 'Peelings', 0.2, 0.8, 10.61]\n",
    "]\n",
    "\n",
    "# Defining the column names\n",
    "columns = ['Crop', 'Residue Type', 'RPR', 'SAF', 'LHV (MJ/kg)']\n",
    "\n",
    "# Create the DataFrame\n",
    "all_residue_values = pd.DataFrame(data, columns=columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87d3dcc6",
   "metadata": {},
   "source": [
    "### Helper function for clipping a raster according to the selected region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "c25bae78",
   "metadata": {},
   "outputs": [],
   "source": [
    "def maskingwithshapefile(shapefile, raster_path):\n",
    "    \n",
    "    with rasterio.open(raster_path) as src:\n",
    "        crs= src.crs\n",
    "        shapefile.crs=crs\n",
    "        clipped, transform = mask(src, shapefile.geometry, crop=True )\n",
    "        # This returns a numpy array on which we will conduct operations.\n",
    "    return clipped"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7145512e",
   "metadata": {},
   "source": [
    "### Helper functions for converting numpy arrays to raster format again to give us a clipped raster ( We have used the memory file data type in rasterio which allows us to create rasters in the active memory without the need to download these to the computer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "fe028136",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For clipped arrays :\n",
    "\n",
    "# The clipped one has extra parameters which ensure that the clipped raster has the correct latitudes and \n",
    "# longitudes acc. to the chosen CRS. This will be mostly useful when we will try to extract rasters with the\n",
    "# correct longitude and latitude values.\n",
    "\n",
    "\n",
    "# Create a function to convert a NumPy array to an in-memory raster\n",
    "def array_to_inmemory_raster_for_clipped(array, transform, crs, shapefile):\n",
    "    \n",
    "    if len(array.shape) == 2:\n",
    "        array = array.reshape((1, array.shape[0], array.shape[1]))  # Add a singleton dimension\n",
    "\n",
    "    _, height, width = array.shape  # Extract the height and width from the shape\n",
    "\n",
    "    # Reshape the array to 2 dimensions\n",
    "    array = array[0]  # Extract the first dimension to remove the extra dimension\n",
    "    \n",
    "    # Correcting the transform for correct axes value representation\n",
    "    \n",
    "    shapefile.geometry.iloc[0]\n",
    "    xmin = shapefile.bounds.iloc[0,3]\n",
    "    ymax = shapefile.bounds.iloc[0,0]\n",
    "    topleft_corner = (xmin,ymax)\n",
    "\n",
    "\n",
    "    a = transform.a\n",
    "    b = transform.b\n",
    "    d = transform.d\n",
    "    e = transform.e\n",
    "\n",
    "    transform = Affine(a, b, ymax, d, e, xmin)\n",
    "    \n",
    "    # Update the transform with corrected values to make sure only the region of shapefile is represented\n",
    "    transform = Affine(transform.a, transform.b, ymax, transform.d, transform.e, xmin)\n",
    "    \n",
    "    # Define the raster metadata\n",
    "    meta = {\n",
    "        'count': 1,\n",
    "        'dtype': array.dtype,\n",
    "        'width': width,\n",
    "        'height': height,\n",
    "        'crs': crs,\n",
    "        'transform': transform\n",
    "    }\n",
    "\n",
    "    # Create an in-memory raster file\n",
    "    memory_file = rasterio.MemoryFile()\n",
    "    with memory_file.open(driver='GTiff', **meta) as dst:\n",
    "        dst.write(array, 1)\n",
    "\n",
    "    return memory_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7c50c445",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For Non-Clipped :\n",
    "\n",
    "# This function might not be as used like the above one and will mostly be used only when converting the full \n",
    "# rasters from one resolution to the other along with resampling of data.\n",
    "\n",
    "def array_to_inmemory_raster_for_non_clipped(array, transform, crs):\n",
    "    if len(array.shape) == 2:\n",
    "        array = array.reshape((1, array.shape[0], array.shape[1]))  # Add a singleton dimension\n",
    "\n",
    "    _, height, width = array.shape  # Extract the height and width from the shape\n",
    "\n",
    "    # Reshape the array to 2 dimensions\n",
    "    array = array[0]  # Extract the first dimension to remove the extra dimension\n",
    "\n",
    "    # Define the raster metadata\n",
    "    meta = {\n",
    "        'count': 1,\n",
    "        'dtype': array.dtype,\n",
    "        'width': width,\n",
    "        'height': height,\n",
    "        'crs': crs,\n",
    "        'transform': transform\n",
    "    }\n",
    "\n",
    "    # Create an in-memory raster file\n",
    "    memory_file = rasterio.MemoryFile()\n",
    "    with memory_file.open(driver='GTiff', **meta) as dst:\n",
    "        dst.write(array, 1)\n",
    "\n",
    "    return memory_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8391bd5",
   "metadata": {},
   "source": [
    "### Helper Functions for extracting select pixels from selected rasters and then consolidating them in a data frame along with their coordinates so that they can be removed/masked later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "459d722b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once we have the rasters we are ready to extract specific pixels from it. The following functions\n",
    "# will help us do exactly that. One of them will extract pixels with specific values and the other will \n",
    "# extract pixels having value above a certain threshold.\n",
    "\n",
    "def coordinates_and_values(raster_path, pixel_values):\n",
    "    \n",
    "    with rasterio.open(raster_path) as src:\n",
    "        df = pd.DataFrame()  \n",
    "        raster_band = src.read(1) \n",
    "        transform_coordinate_conversion = src.transform\n",
    "\n",
    "        for i in pixel_values:\n",
    "            rows, cols = np.where(raster_band == i)  # extract row and column numbers for each pixel\n",
    "            coords = transform_coordinate_conversion * (cols, rows)\n",
    "            lon, lat = coords[0], coords[1]  \n",
    "\n",
    "            values = raster_band[rows,cols]\n",
    "            df_i = pd.DataFrame({'lon': lon, 'lat': lat, 'pixel': values, 'row': rows, 'col': cols}) \n",
    "            df = pd.concat([df, df_i], ignore_index=True)  # append to get DataFrame of lon, lat, row, col, and \n",
    "            # pixel\n",
    "\n",
    "    return df\n",
    "\n",
    "def coordinates_and_threshold(raster_path, threshold):\n",
    "    \n",
    "    with rasterio.open(raster_path) as src:\n",
    "        df = pd.DataFrame()  \n",
    "        raster_band = src.read(1)  \n",
    "        transform_coordinate_conversion = src.transform\n",
    "\n",
    "        rows, cols = np.where(raster_band > threshold)  # find row and column indices for pixel values above \n",
    "        # threshold\n",
    "        coords = transform_coordinate_conversion * (cols, rows)\n",
    "        lon, lat = coords[0], coords[1] \n",
    "\n",
    "        values = raster_band[rows, cols]  # extract pixel values above threshold\n",
    "        df = pd.DataFrame({'lon': lon, 'lat': lat, 'row': rows, 'col': cols, 'pixel': values})  \n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b1c307c",
   "metadata": {},
   "source": [
    "### Helper functions for converting the clipping to shape file and conversion back to raster in a single function . It also converts the data frames generated above to a GeoDataFrame (gdf) so that the stored locations can be used to remove/mask these pixels later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4c446da4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clipper(shapefile, raster_path):\n",
    "    \n",
    "    with rasterio.open(raster_path) as src:\n",
    "        transform_store = src.transform\n",
    "        crs_store = src.crs\n",
    "        \n",
    "    masked = maskingwithshapefile(shapefile , raster_path)\n",
    "    output_raster = array_to_inmemory_raster_for_clipped(masked, transform_store, crs_store, shapefile)\n",
    "    \n",
    "    return output_raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d0f80538",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating this function since both the value and threshold functions give us dataframes and we will need\n",
    "# geodataframes containing the geometry column with the Point values for each pixel.\n",
    "\n",
    "def convert_df_to_gdf(dataframe):\n",
    "    \n",
    "    # Convert the DataFrame of coordinates to a GeoDataFrame\n",
    "    geometry = [Point(lon, lat) for lon, lat in zip(dataframe['lon'], dataframe['lat'])]\n",
    "    coordinates_gdf = gpd.GeoDataFrame(dataframe, geometry=geometry)\n",
    "    \n",
    "    return coordinates_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec7f7b81",
   "metadata": {},
   "source": [
    "### Helper function for conversion of rasters from higher to lower resolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b5708a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Also since the crop data is in a lower resolution than the AEZ classsification and other data so we will also \n",
    "# create a fn. for the conversion of Resolution from a higher to lower resolution along with resampling of the \n",
    "# pixel values. The current method of resampling used is 'Mode' method.\n",
    "\n",
    "def resolution_converter_mode(raster_path , resampling_method):\n",
    "    \n",
    "    downscale_factor = 10 \n",
    "    with rasterio.open(raster_path) as dataset:\n",
    "\n",
    "        # Compute the downsampled shape and final resolution based on the downscale factor\n",
    "        downsampled_height = int(dataset.height / downscale_factor)\n",
    "        downsampled_width = int(dataset.width / downscale_factor)\n",
    "        downsampled_shape = (dataset.count, downsampled_height, downsampled_width)\n",
    "        initial_resolution = dataset.res[0]\n",
    "        final_resolution = initial_resolution * downscale_factor\n",
    "        crs_final = dataset.crs\n",
    "\n",
    "        # Resample the data to the target shape\n",
    "        data = dataset.read(\n",
    "            out_shape=downsampled_shape,\n",
    "            resampling=resampling_method\n",
    "        )\n",
    "\n",
    "     # Compute the scale factors for the image transform\n",
    "    scale_x =1/(initial_resolution / final_resolution)\n",
    "    scale_y = 1/(-initial_resolution / -final_resolution)\n",
    "    \n",
    "    initial_transform = Affine(initial_resolution, 0.0, -180.0, 0.0, -initial_resolution, 90.0 )\n",
    "\n",
    "    # Scale the image transform\n",
    "    transform = initial_transform * initial_transform.scale(scale_x, scale_y)\n",
    "\n",
    "    # Reshape the data array to have two dimensions (height x width)\n",
    "    reshaped_data = data[0, :, :]\n",
    "\n",
    "    return array_to_inmemory_raster_for_non_clipped(reshaped_data, transform, crs_final)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc774e3e",
   "metadata": {},
   "source": [
    "### Helper function for removal of accumulated pixels in the GeoDataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "d172dfe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To explain it a bit more, what we do is get those pixels and assign the pixels to be removed a nodata value.\n",
    "\n",
    "def remove_pixels(raster_path, shapefile, geodataframe):\n",
    "    with rasterio.open(raster_path) as src:\n",
    "        # Mask the raster using the shapefile boundary\n",
    "        clipped, transform = mask(src, shapefile.geometry, crop=True)\n",
    "        crs = src.crs\n",
    "\n",
    "    raster_used = array_to_inmemory_raster_for_clipped(clipped, transform, crs, shapefile)\n",
    "\n",
    "    dataset = rasterio.open(raster_used)\n",
    "    \n",
    "    # Load the GeoDataFrame\n",
    "    gdf = geodataframe\n",
    "    \n",
    "    # Convert the row and column numbers to pixel coordinates\n",
    "    pixel_coords = zip(gdf['row'], gdf['col'])\n",
    "    \n",
    "    # Create a mask with the same dimensions as the raster\n",
    "    masked = np.zeros((dataset.height, dataset.width), dtype=np.uint8)\n",
    "    \n",
    "    # Set the pixels at the specified coordinates to a nodata value\n",
    "    nodata_value = 255\n",
    "    for row, col in pixel_coords:\n",
    "        masked[row, col] = nodata_value\n",
    "    \n",
    "    # Read the raster data\n",
    "    data = dataset.read(1)\n",
    "    \n",
    "    # Apply the mask\n",
    "    data = np.where(masked == nodata_value, np.nan, data)\n",
    "    \n",
    "    return data\n",
    "\n",
    "\n",
    "# So what this function does, is that it takes a sample raster path as input first clips it. Then it creates a\n",
    "# mask layer using the sample raster while assigning nodata values to the points we want removed. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf2cb724",
   "metadata": {},
   "source": [
    "### Helper function for going over the selected region and identifying which crop to grow to maximize energy extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "1816bf6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_max_for_each_pixel(time_period, climate_model, rcp, water_supply_future, input_level,\n",
    "                            shapefile, geodataframe):\n",
    "    # Initialize dictionaries to store the data for each crop and their residues\n",
    "    crop_data = {}\n",
    "    max_values = None\n",
    "    max_crops = None\n",
    "    crop_residue_sum_dict = {}  # Dictionary to store crop_residue_sum_array for each crop\n",
    "    \n",
    "    filtered_potential_yield = potential_yield[(potential_yield['Time Period'] == time_period) &\n",
    "                                               (potential_yield['Climate Model'] == climate_model) &\n",
    "                                               (potential_yield['RCP'] == rcp) &\n",
    "                                               (potential_yield['Water Supply'] == water_supply_future) &\n",
    "                                               (potential_yield['Input Level'] == input_level)]\n",
    "    required_potential_yields = filtered_potential_yield[['Crop', 'Download URL']]\n",
    "    \n",
    "    # For defining size of the xarray\n",
    "    with rasterio.open(potential_yield.iloc[2, 14].strip()) as src:\n",
    "        clipped_shapefile_init, transform_init = mask(src, shapefile.geometry, crop=True)\n",
    "        clipped_shapefile_init = remove_band_dimension(clipped_shapefile_init)\n",
    "        lats_init, lons_init = get_lat_lon_from_transform(transform_init, clipped_shapefile_init.shape)\n",
    "    \n",
    "    # Iterate over the global rasters\n",
    "    for crop in potential_yield['Crop'].unique():\n",
    "        # Find the correct raster path for the rasters you want to access\n",
    "        raster_path = required_potential_yields[required_potential_yields['Crop'] == crop]['Download URL'].values[0].strip()\n",
    "        \n",
    "        # Remove pixels from the raster\n",
    "        data = remove_pixels(raster_path, shapefile, geodataframe)\n",
    "        \n",
    "        # Multiply the values with the corresponding RPR, SAF, and LHV\n",
    "        residue_rows = all_residue_values.loc[all_residue_values['Crop'] == crop]\n",
    "        crop_residue_sum_array = np.zeros_like(data)\n",
    "        for _, residue_row in residue_rows.iterrows():\n",
    "            LHV = residue_row['LHV (MJ/kg)']\n",
    "            SAF = residue_row['SAF']\n",
    "            RPR = residue_row['RPR']\n",
    "            data_final = data * LHV * SAF * RPR\n",
    "            \n",
    "            # Handle NaN values by converting them to zero\n",
    "            data_final = np.nan_to_num(data_final, nan=0)\n",
    "            positive_data_mask = data_final >= 0\n",
    "            data_final = data_final * positive_data_mask\n",
    "            \n",
    "            crop_residue_sum_array += data_final \n",
    "        \n",
    "        # Store the data for the current crop and their residues in dictionaries\n",
    "        crop_residue_sum_dict[crop] = crop_residue_sum_array\n",
    "        \n",
    "        # Check if max_values array is None (first iteration)\n",
    "        if max_values is None:\n",
    "            max_values = crop_residue_sum_array\n",
    "            max_crops = np.empty_like(crop_residue_sum_array).astype(str)\n",
    "            max_crops[crop_residue_sum_array == None] = None\n",
    "            max_crops[crop_residue_sum_array != None] = crop\n",
    "        else:\n",
    "            # Find the maximum value between the current raster and previous maximum values\n",
    "            max_mask = crop_residue_sum_array > max_values\n",
    "            max_values[max_mask] = crop_residue_sum_array[max_mask]\n",
    "            max_crops[max_mask] = crop\n",
    "    \n",
    "    # Replace instances of \"Alfalfa\" with None/NaN in max_crops\n",
    "    max_crops[max_crops == \"Alfalfa\"] = None\n",
    "\n",
    "    # Create xarray Dataset to hold all the individual biomass potentials for each crop\n",
    "    biomass_potentials_dataset = xr.Dataset()\n",
    "\n",
    "    # Add the crop_residue_sum_array as a variable for each crop\n",
    "    for crop, crop_residue_sum_array in crop_residue_sum_dict.items():\n",
    "        # Remove spaces and special characters from crop name to make it a valid variable name\n",
    "        variable_name = crop.replace(\" \", \"_\").replace(\"-\", \"_\").replace(\"(\", \"\").replace(\")\", \"\").replace(\",\", \"\")\n",
    "    \n",
    "        # Create a DataArray for the crop_residue_sum_array\n",
    "        data_array = xr.DataArray(\n",
    "            crop_residue_sum_array,\n",
    "            dims=('y', 'x'),\n",
    "            coords={'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "            attrs={'units': 'PetaJoules'}\n",
    "        )\n",
    "    \n",
    "        # Add the DataArray as a variable to the dataset\n",
    "        biomass_potentials_dataset[variable_name] = data_array\n",
    "    \n",
    "        # Add the yield sum as an attribute for the variable\n",
    "        biomass_potentials_dataset[variable_name].attrs['sum'] = np.nansum(crop_residue_sum_array)\n",
    "\n",
    "    # Add the max_values as a variable to the dataset\n",
    "    biomass_potentials_dataset['max_values'] = xr.DataArray(\n",
    "        max_values,\n",
    "        dims=('y', 'x'),\n",
    "        coords={'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)},\n",
    "        attrs={'units': 'PetaJoules'}\n",
    "    )\n",
    "    biomass_potentials_dataset['max_values'].attrs['sum'] = np.nansum(max_values)\n",
    "    \n",
    "    # Add the max_crops as a variable to the dataset\n",
    "    biomass_potentials_dataset['max_crops'] = xr.DataArray(\n",
    "    max_crops,\n",
    "    dims=('y', 'x'),\n",
    "    coords={'latitude': (('x', 'y'), lats_init), 'longitude': (('x', 'y'), lons_init)}\n",
    "        ).astype(str)\n",
    "    \n",
    "    # Add attributes for the sum of max_values and units\n",
    "    biomass_potentials_dataset.attrs['net_sum'] = np.nansum(max_values)\n",
    "    biomass_potentials_dataset.attrs['units'] = 'PetaJoules'\n",
    "\n",
    "    return biomass_potentials_dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "378abbed",
   "metadata": {},
   "source": [
    "### Helper function for obtaining harvested area and the net pixel area for each pixel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "af3e72a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To be able to find the harvested area per pixel we need to sum up the harvested area for each crop in each pixel\n",
    "# and store it in the form of an numpy array. \n",
    "\n",
    "def get_net_harvested_area(shapefile, geodataframe):\n",
    "    \n",
    "    net_harvested_area_obtained = None  # Initialize the net harvested area array\n",
    "    \n",
    "    filtered_harvested_area = harvested_area[(harvested_area['Time Period'] == 2010) &\n",
    "                                                   (harvested_area['Water Supply'] == 'Total')]\n",
    "    required_harvested_area = filtered_harvested_area['Download URL']\n",
    "    \n",
    "    for required_url in required_harvested_area: # This will loop pover all the available crops in cropland.\n",
    "        required_url = required_url.strip()\n",
    "        \n",
    "        with rasterio.open(required_url) as src:\n",
    "            # Assuming shapefile is used for clipping\n",
    "            # Adjust accordingly if using a different method for clipping\n",
    "            \n",
    "            # Clip the raster using the shapefile\n",
    "            data = remove_pixels(required_url, shapefile , geodataframe)\n",
    "            \n",
    "            # Get the harvested area values from the clipped data\n",
    "            harvested_area_obtained = data  # Adjust this line if necessary\n",
    "            \n",
    "            # Update the net harvested area array\n",
    "            if net_harvested_area_obtained is None:\n",
    "                net_harvested_area_obtained = harvested_area_obtained\n",
    "            else:\n",
    "                net_harvested_area_obtained += harvested_area_obtained\n",
    "    \n",
    "    return net_harvested_area_obtained*1000 # To make the unit as hectares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "190cd571",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simply using a simple cosine of latitude approximation to account for the curvature of the Earth. The\n",
    "# consequences of making this simple approximation have been explored in detail in the paper.\n",
    "\n",
    "def extract_pixel_area(raster_path, shapefile):\n",
    "    with rasterio.open(raster_path) as src:\n",
    "    # Clip the raster using the shapefile boundaries\n",
    "        clipped_data, clipped_transform = mask(src, shapefile.geometry, crop=True)\n",
    "\n",
    "    # Get the pixel dimensions of the clipped raster\n",
    "    pixel_width = clipped_transform[0]\n",
    "    pixel_height = clipped_transform[4]\n",
    "\n",
    "    # Calculate the area of each pixel based on latitude\n",
    "    pixel_area = np.zeros_like(clipped_data, dtype=np.float32)\n",
    "\n",
    "    for band in range(clipped_data.shape[0]):\n",
    "        for row in range(clipped_data.shape[1]):\n",
    "            for col in range(clipped_data.shape[2]):\n",
    "                if clipped_data[band, row, col] != src.nodata:\n",
    "                    lat = rasterio.transform.xy(clipped_transform, row, col, offset='center')[1]\n",
    "                    lat_rad = np.radians(lat)\n",
    "                    pixel_area[band, row, col] = (np.abs(pixel_width * pixel_height * (111319.9)**2)*(np.cos(lat_rad)))/(10000) # conversion to hectares\n",
    "    return pixel_area"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c6f754",
   "metadata": {},
   "source": [
    "#### Final output function for marginal land data. \n",
    "- It outputs an xarray containing the energy yields of each crop for the selected geography (These are not multiplied by the area since the area is the same and this saves compute). \n",
    "- It also outputs the net biomass energy potential from the marginal land using the the crop energy yields which give the max energy extractable from each pixel. \n",
    "- Finally it outputs the sum of this final array described in the second point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "dce54f4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_biomass_potential_for_marginal(shapefile,time_period, climate_model, rcp, water_supply_future,\n",
    "                                       input_level):\n",
    "    \n",
    "    # Selecting and correcting the resolution of the selected files\n",
    "    \n",
    "    initial_aez_chosen = aez_classification[(aez_classification['Time Period']== time_period)&\n",
    "                                         (aez_classification['RCP']== rcp)]\n",
    "    initial_aez_raster = initial_aez_chosen['Download URL'].values[0].strip()\n",
    "    \n",
    "    final_aez_raster = resolution_converter_mode(initial_aez_raster , Resampling.mode)\n",
    "    \n",
    "    initial_exclusion = exclusion_areas.iloc[0,9].strip()\n",
    "    \n",
    "    final_exclusion_raster = resolution_converter_mode(initial_exclusion , Resampling.mode)\n",
    "    \n",
    "    initial_tree = tree_cover_share.iloc[0,9].strip()\n",
    "    \n",
    "    final_tree_cover_raster = resolution_converter_mode(initial_tree , Resampling.average)\n",
    "    \n",
    "    pasture_raster = \"./dataset/pasture_cubic_reproject.tif\" # the reproject here is just to make the pixel size 0.083333 instead of 0.08328\n",
    "    \n",
    "    clipped_aez = clipper(shapefile , final_aez_raster)\n",
    "    clipped_exclusion = clipper(shapefile, final_exclusion_raster)# Now remove 2 to 7 value range.\n",
    "    clipped_tree = clipper(shapefile, final_tree_cover_raster)# Threshold == 50\n",
    "    clipped_pasture = clipper(shapefile, pasture_raster)\n",
    "    \n",
    "    coordinates_aez= coordinates_and_values(clipped_aez, [49,50,52,53,55,56,57])\n",
    "    coordinates_exclusion = coordinates_and_values(clipped_exclusion, [2,3,4,5,6,7])\n",
    "    coordinates_tree = coordinates_and_threshold(clipped_tree, 50) \n",
    "    coordinates_pasture = coordinates_and_threshold(clipped_pasture, 0.5)\n",
    "    \n",
    "    coordinates_all = pd.concat([coordinates_aez, coordinates_exclusion, coordinates_tree,coordinates_pasture],\n",
    "                                ignore_index=True)  # append to get DataFrame of lon, lat, row, col, and pixel\n",
    "    \n",
    "    gdf_final = convert_df_to_gdf(coordinates_all)\n",
    "    \n",
    "    biomass_potential_xarray = find_max_for_each_pixel(time_period, climate_model, rcp,\n",
    "                                                       water_supply_future, input_level,\n",
    "                                                       shapefile, gdf_final)\n",
    "    \n",
    "    harvested_area_from_shapefile = get_net_harvested_area(shapefile, gdf_final)\n",
    "    \n",
    "    net_area = extract_pixel_area(potential_yield.iloc[2,14].strip(), shapefile) # Just a reference raster so doesn't matter\n",
    "    \n",
    "    # Using sample raster potential_yield.iloc[2,14] to extract the base transform and crs used for all\n",
    "    with rasterio.open(potential_yield.iloc[2,14].strip()) as src:\n",
    "        standard_transform = src.transform\n",
    "        standard_crs = src.crs\n",
    "    \n",
    "    remaining_area = np.subtract(net_area,harvested_area_from_shapefile)\n",
    "    final_potential = np.multiply(biomass_potential_xarray['max_values'].values, remaining_area)*(10**-5)# Unit conversion MJ to Joules and 10 Kg to Kg and \n",
    "    # then 10**-12 for peta joules\n",
    "    total_biomass_marginal_potential = np.nansum(final_potential)\n",
    "\n",
    "    \n",
    "    return total_biomass_marginal_potential, final_potential, biomass_potential_xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c98f5704",
   "metadata": {},
   "source": [
    "#### Total Raw Biomass Energy Potential : \n",
    "- This is the final total raw biomass energy potential function which gives us the sum of energy potential from each pixel for cropland and marginal land, which is the final intended output of our tool. \n",
    "- It also produces a final xarray which contains everything under their respective headings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "158e69c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_total_biomass_potential(shapefile, time_period, climate_model, rcp, water_supply_future, input_level, water_supply_2010):\n",
    "    # Calculate biomass potential for cropland\n",
    "    cropland_dataset = future_potential_cropland(time_period, climate_model, rcp, water_supply_future, input_level, shapefile, water_supply_2010)\n",
    "\n",
    "    # Calculate biomass potential for marginal land\n",
    "    marginal_land_potential, marginal_land_array, marginal_land_dataset = get_biomass_potential_for_marginal(shapefile, time_period, climate_model, rcp, water_supply_future, input_level)\n",
    "\n",
    "    # Create an empty dataset to hold the outputs\n",
    "    total_biomass_dataset = xr.Dataset(coords=cropland_dataset.coords)\n",
    "\n",
    "    # Add attributes for the total sum of attributes for each pixel\n",
    "    total_sum_cropland = cropland_dataset.attrs['net_sum in PJ']\n",
    "    total_sum_marginal_land = marginal_land_potential\n",
    "    total_biomass_dataset.attrs['net_energy_potential'] = total_sum_cropland + total_sum_marginal_land\n",
    "    total_biomass_dataset.attrs['final_energy_array'] = marginal_land_array + cropland_dataset['Combined'].values\n",
    "\n",
    "    return total_biomass_dataset , cropland_dataset , marginal_land_dataset , marginal_land_array"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83f46dd9",
   "metadata": {},
   "source": [
    "### Functions to visualize the raster with pixel values shown & to display the crop selected int eh case of marginal lands"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "dfe02c88",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assuming the rest of your code remains the same\n",
    "def bokeh_plot(shapefile, array ):\n",
    "    \n",
    "    with rasterio.open(potential_yield.iloc[2,14].strip()) as src:\n",
    "            standard_transform = src.transform \n",
    "            standard_crs= src.crs\n",
    "    # Convert the GeoDataFrame to GeoJSONDataSource\n",
    "    \n",
    "    test = array_to_inmemory_raster_for_clipped(array, standard_transform,standard_crs,shapefile)\n",
    "    \n",
    "    geojson = gpd.GeoSeries(shapefile.geometry).to_json()\n",
    "    geojson_data = GeoJSONDataSource(geojson=geojson)\n",
    "\n",
    "    # Read the raster data using rasterio.open\n",
    "    raster_path = test\n",
    "    src = rasterio.open(raster_path)\n",
    "\n",
    "    # Read the raster data into a NumPy array\n",
    "    data = src.read(1)\n",
    "\n",
    "    # Flip the raster data along the y-axis to align with Bokeh's coordinate system\n",
    "    data = np.flipud(data)\n",
    "\n",
    "    # Calculate the minimum and maximum x and y coordinates\n",
    "    min_x, min_y, max_x, max_y = src.bounds\n",
    "\n",
    "    # Calculate the minimum and maximum values for the color mapper using the 'combined' array\n",
    "    combined_min = array.min().item()\n",
    "    combined_max = array.max().item()\n",
    "\n",
    "    # Create a new color mapper for the raster data\n",
    "    mapper = LinearColorMapper(palette=bp.Blues8, low=combined_min, high=combined_max)\n",
    "\n",
    "\n",
    "    # Create a Bokeh figure\n",
    "    p = figure(width=600, height=600, x_axis_label='Longitude', y_axis_label='Latitude', title='Interactive Raster')\n",
    "\n",
    "    # Add the raster image to the figure using the 'image' glyph\n",
    "    p.image(image=[data], x=min_x, y=min_y, dw=(max_x - min_x), dh= (max_y - min_y), color_mapper=mapper)\n",
    "\n",
    "    # Add color bar to the plot\n",
    "    color_bar = ColorBar(color_mapper=mapper, location=(0, 0))\n",
    "    p.add_layout(color_bar, 'right')\n",
    "\n",
    "# Add the boundary of Spain to the plot as a polygon\n",
    "    p.patches('xs', 'ys', source=geojson_data, line_color='black', fill_alpha=0)\n",
    "\n",
    "    # Add hover tool to display the pixel value when hovering over the raster\n",
    "    hover = HoverTool(tooltips=[('Value', '@image')], mode='mouse')\n",
    "    p.add_tools(hover)\n",
    "\n",
    "    # Show the plot\n",
    "    bpl.show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2611cb9e-03d8-4bb6-aed2-704c74d6851c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bokeh_max_min_plot(shapefile, array):\n",
    "    with rasterio.open(potential_yield.iloc[2,14].strip()) as src:\n",
    "            standard_transform = src.transform \n",
    "            standard_crs= src.crs\n",
    "\n",
    "    test = array_to_inmemory_raster_for_clipped(array, standard_transform, standard_crs, shapefile)\n",
    "    geojson = gpd.GeoSeries(shapefile.geometry).to_json()\n",
    "    geojson_data = GeoJSONDataSource(geojson=geojson)\n",
    "\n",
    "    src = rasterio.open(test)\n",
    "    data = src.read(1)\n",
    "    data = np.flipud(data)\n",
    "\n",
    "    min_x, min_y, max_x, max_y = src.bounds\n",
    "\n",
    "    # Ignore 0 as the minimum value and take the next higher value\n",
    "    # Filter the array to include only non-zero values, then find the minimum\n",
    "    non_zero_array = array[array > 0]\n",
    "    combined_min = non_zero_array.min().item() if non_zero_array.size > 0 else array.min().item()\n",
    "    combined_max = array.max().item()\n",
    "\n",
    "    palette = plasma(256)  # Using the Magma palette, or choose another\n",
    "\n",
    "    mapper = LinearColorMapper(palette=palette, low=combined_min, high=combined_max)\n",
    "\n",
    "    p = figure(width=600, height=600, x_axis_label='Longitude', y_axis_label='Latitude', title='Interactive Raster')\n",
    "    p.image(image=[data], x=min_x, y=min_y, dw=(max_x - min_x), dh=(max_y - min_y), color_mapper=mapper)\n",
    "\n",
    "    color_bar = ColorBar(color_mapper=mapper, location=(0, 0))\n",
    "    p.add_layout(color_bar, 'right')\n",
    "\n",
    "    p.patches('xs', 'ys', source=geojson_data, line_color='black', fill_alpha=0)\n",
    "\n",
    "    hover = HoverTool(tooltips=[('Value', '@image')], mode='mouse')\n",
    "    p.add_tools(hover)\n",
    "\n",
    "    bpl.show(p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "d65342be",
   "metadata": {},
   "outputs": [],
   "source": [
    "def crop_show(crop_array,shapefile):\n",
    "    \n",
    "    crop_names = np.unique(crop_array)\n",
    "    crop_indices = {crop: i for i, crop in enumerate(crop_names)}\n",
    "    color_index = np.vectorize(crop_indices.get)(crop_array)\n",
    "     # Define colormap with an additional color for 'None'\n",
    "    cmap = plt.get_cmap('tab20b', len(crop_indices))\n",
    "    \n",
    "    with rasterio.open(potential_yield.iloc[2, 14].strip()) as src:\n",
    "        standard_transform = src.transform \n",
    "        standard_crs= src.crs\n",
    "    \n",
    "    # Get the pixel size from the standard_transform\n",
    "    pixel_size_x = standard_transform.a\n",
    "    pixel_size_y = standard_transform.e\n",
    "    \n",
    "    graph_bounds = shapefile.bounds\n",
    "    \n",
    "    # Get the boundary values from the shapefile\n",
    "    xmin = graph_bounds['minx'].min()\n",
    "    xmax = graph_bounds['maxx'].max()\n",
    "    ymin = graph_bounds['miny'].min()\n",
    "    ymax = graph_bounds['maxy'].max()\n",
    "    \n",
    "    # Calculate the number of rows and columns in the raster\n",
    "    rows, cols = color_index.shape\n",
    "    \n",
    "    # Create the transformation matrix for the raster\n",
    "    transform = Affine(pixel_size_x, 0, xmin,\n",
    "                       0, pixel_size_y, ymax)\n",
    "    \n",
    "    # Create a memory file to store the raster\n",
    "    with MemoryFile() as memfile:\n",
    "        # Create a new raster dataset\n",
    "        with memfile.open(driver='GTiff', height=rows, width=cols, count=1,\n",
    "                          dtype=color_index.dtype, crs=standard_crs,\n",
    "                          transform=transform) as dataset:\n",
    "            # Write the color_index array to the raster dataset\n",
    "            dataset.write(color_index, 1)\n",
    "    \n",
    "            # Plot setup\n",
    "            fig, ax = plt.subplots()\n",
    "    \n",
    "            # Plot the raster with shapefile boundaries\n",
    "            show((dataset, 1), ax=ax, cmap=cmap, vmin=0, vmax=len(crop_indices)-1)\n",
    "    \n",
    "            # Plot the shapefile boundary\n",
    "            shapefile.plot(ax=ax, facecolor='none', edgecolor='black')\n",
    "            \n",
    "            legend_patches = [mpatches.Patch(color=cmap(i), label=crop) for crop, i in crop_indices.items()]\n",
    "            plt.legend(handles=legend_patches, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)\n",
    "    \n",
    "            plt.title(\"Marginal Crop Distribution Across The Selected Region\")\n",
    "            plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e4ae3df",
   "metadata": {},
   "source": [
    "### These are the final visualisation functions which output the net raw biomass energy potential from the marginal and the cropland respectively and show them with an interactive plotly graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8b039618",
   "metadata": {},
   "outputs": [],
   "source": [
    "def graph_plotter_cropland(shapefile, climate_model, water_supply_future, input_level):\n",
    "    time_periods = ['2000', '2010', '2011-2040', '2041-2070', '2071-2100']\n",
    "    RCPs = ['RCP2.6', 'RCP4.5', 'RCP6.0', 'RCP8.5']\n",
    "    \n",
    "    xarrays = {}  # Dictionary to store the xarray for each RCP and time period\n",
    "    \n",
    "    # Colors from the Google logo: Blue, Red, Yellow, Blue-Green\n",
    "    colors = ['#4285F4', '#DB4437', '#F4B400', '#0F9D58']\n",
    "\n",
    "    fig = go.Figure()\n",
    "    \n",
    "    value_1 = biomass_potential_past(shapefile, 2000, 'Total')\n",
    "    value_2 = biomass_potential_past(shapefile, 2010, 'Total')\n",
    "    \n",
    "    xarrays[('2000')] = value_1\n",
    "    xarrays[('2010')] = value_2\n",
    "\n",
    "    initial_potential_1 = value_1.attrs['Net Potential in PetaJ']\n",
    "    initial_potential_2 = value_2.attrs['Net Potential in PetaJ']\n",
    "\n",
    "    # Add traces for 2000 and 2010 with different colors\n",
    "    fig.add_trace(go.Bar(x=['2000'], y=[initial_potential_1], name='2000', marker_color='#000000'))\n",
    "    fig.add_trace(go.Bar(x=['2010'], y=[initial_potential_2], name='2010', marker_color='#808080'))\n",
    "\n",
    "    x_values = ['2000', '2010']  # Initialize x-axis values\n",
    "\n",
    "    for i, RCP in enumerate(RCPs):\n",
    "        biomass_potentials = []  # Initialize biomass_potentials for each time period\n",
    "\n",
    "        for j, time_period in enumerate(time_periods[2:]):\n",
    "            \n",
    "            value_calculator = future_potential_cropland(time_period, climate_model, RCP, water_supply_future, input_level, shapefile, 'Total')\n",
    "            \n",
    "            potential_value = value_calculator.attrs['net_sum in PJ']\n",
    "            biomass_potentials.append(potential_value)\n",
    "            \n",
    "            # Store the xarray in the xarrays dictionary with a tuple key (RCP, time_period)\n",
    "            xarrays[(RCP, time_period)] = value_calculator\n",
    "\n",
    "        # Assign the corresponding color from the Google logo to each time period\n",
    "        color = colors[i % len(colors)]\n",
    "\n",
    "        fig.add_trace(go.Bar(x=time_periods[2:], y=biomass_potentials, name=RCP, marker_color=color))\n",
    "\n",
    "    fig.update_layout(\n",
    "        barmode='group',\n",
    "        xaxis_title='Years',\n",
    "        yaxis_title='Biomass Energy Potential from Cropland Land <br> in PetaJoules',\n",
    "        title='Cropland Biomass Energy Potential from different RCPs'\n",
    "    )\n",
    "    \n",
    "    fig.show()\n",
    "\n",
    "    return fig, xarrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32aa378c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def graph_plotter_marginal(shapefile, climate_model, water_supply_future, input_level):\n",
    "    time_periods = ['2011-2040', '2041-2070', '2071-2100']\n",
    "    RCPs = ['RCP2.6', 'RCP4.5', 'RCP6.0', 'RCP8.5']\n",
    "\n",
    "    fig = go.Figure()\n",
    "    xarrays = {}  # Dictionary to store the xarray for each RCP and time period\n",
    "    final_potentials = {}  # Dictionary to store the final potential array for each RCP\n",
    "\n",
    "    for i, RCP in enumerate(RCPs):\n",
    "        biomass_potentials = []  # Initialize biomass_potentials for each time period\n",
    "\n",
    "        for time_period in time_periods:\n",
    "            total_biomass_marginal_potential, final_potential, xarray = get_biomass_potential_for_marginal(\n",
    "                shapefile, time_period, climate_model, RCP, water_supply_future, input_level\n",
    "            )\n",
    "            biomass_potentials.append(total_biomass_marginal_potential)\n",
    "            xarrays[(RCP, time_period)] = xarray\n",
    "            final_potentials[(RCP, time_period)] = final_potential\n",
    "\n",
    "        fig.add_trace(go.Bar(x=time_periods, y=biomass_potentials, name=RCP))\n",
    "\n",
    "    fig.update_layout(\n",
    "        barmode='group',\n",
    "        xaxis_title='Time Periods',\n",
    "        yaxis_title='Biomass Energy Potential from Marginal Land <br> in PetaJoules',\n",
    "        title='Marginal Land Biomass Energy Potential from different RCPs in Time Periods'\n",
    "    )\n",
    "    \n",
    "    fig.show()\n",
    "\n",
    "    return fig, xarrays, final_potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6e606128",
   "metadata": {},
   "outputs": [],
   "source": [
    "def graph_plotter_all(shapefile, climate_model, water_supply_future, input_level, water_supply_2010):\n",
    "    time_periods = ['2011-2040','2041-2070', '2071-2100']\n",
    "    RCPs = ['RCP2.6', 'RCP4.5', 'RCP6.0', 'RCP8.5']\n",
    "    colors = ['#4285F4', '#DB4437', '#F4B400', '#0F9D58']\n",
    "    \n",
    "    fig_crop = go.Figure()\n",
    "    fig_marg = go.Figure()\n",
    "    fig_total = go.Figure()\n",
    "    \n",
    "    \n",
    "    # For cropland \n",
    "    \n",
    "    xarrays_crop = {}  \n",
    "    value_1 = biomass_potential_past(shapefile, 2000, 'Total')\n",
    "    value_2 = biomass_potential_past(shapefile, 2010, 'Total')\n",
    "    \n",
    "    xarrays_crop[('2000')] = value_1\n",
    "    xarrays_crop[('2010')] = value_2\n",
    "\n",
    "    initial_potential_1 = value_1.attrs['Net Potential in PetaJ']\n",
    "    initial_potential_2 = value_2.attrs['Net Potential in PetaJ']\n",
    "\n",
    "    # Add traces for 2000 and 2010 with different colors\n",
    "    fig_crop.add_trace(go.Bar(x=['2000'], y=[initial_potential_1], name='2000', marker_color='#000000'))\n",
    "    fig_crop.add_trace(go.Bar(x=['2010'], y=[initial_potential_2], name='2010', marker_color='#808080'))\n",
    "\n",
    "    x_values = ['2000', '2010']  # Initialize x-axis values\n",
    "    \n",
    "    # For marginal land\n",
    "    \n",
    "    xarrays_marg = {}  # Dictionary to store the xarray for each RCP and time period\n",
    "    final_potentials_marg = {}  # Dictionary to store the final potential array for each RCP\n",
    "    \n",
    "    # For total land \n",
    "    \n",
    "    arrays_total = {}\n",
    "\n",
    "    for i, RCP in enumerate(RCPs):\n",
    "        biomass_potentials_crop = []  # Initialize biomass_potentials for each time period\n",
    "        biomass_potentials_marg = []\n",
    "        biomass_potentials_total = []\n",
    "\n",
    "        for j, time_period in enumerate(time_periods):\n",
    "            \n",
    "            total_energy_potential, cropland_dataset, marginal_dataset, marginal_array   = get_total_biomass_potential(\n",
    "                shapefile, time_period, climate_model, RCP, water_supply_future, input_level, water_supply_2010\n",
    "            )\n",
    "            \n",
    "            # Total Land\n",
    "            biomass_potentials_total.append(total_energy_potential.attrs['net_energy_potential'])\n",
    "            arrays_total[(RCP, time_period)] = total_energy_potential.attrs['final_energy_array']\n",
    "            \n",
    "            # Marginal Land\n",
    "            biomass_potentials_marg.append(np.nansum(marginal_array))\n",
    "            xarrays_marg[(RCP, time_period)] = marginal_dataset\n",
    "            final_potentials_marg[(RCP, time_period)] = marginal_array\n",
    "            \n",
    "            # Cropland\n",
    "            potential_value = cropland_dataset.attrs['net_sum in PJ']\n",
    "            biomass_potentials_crop.append(potential_value)\n",
    "            xarrays_crop[(RCP, time_period)] = cropland_dataset\n",
    "            \n",
    "\n",
    "        # Assign the corresponding color from the Google logo to each time period\n",
    "        color = colors[i % len(colors)]\n",
    "        fig_crop.add_trace(go.Bar(x=time_periods, y=biomass_potentials_crop, name=RCP, marker_color=color))\n",
    "        fig_marg.add_trace(go.Bar(x=time_periods, y=biomass_potentials_marg, name=RCP))\n",
    "        fig_total.add_trace(go.Bar(x=time_periods, y=biomass_potentials_total, name=RCP))\n",
    "\n",
    "    fig_total.update_layout(\n",
    "        barmode='group',\n",
    "        xaxis_title='Time Periods',\n",
    "        yaxis_title='Biomass Energy Potential from Total Land <br> in PetaJoules',\n",
    "        title='Total Biomass Energy Potential from different RCPs in Time Periods'\n",
    "    )\n",
    "\n",
    "    fig_crop.update_layout(\n",
    "        barmode='group',\n",
    "        xaxis_title='Years',\n",
    "        yaxis_title='Biomass Energy Potential from Cropland Land <br> in PetaJoules',\n",
    "        title='Cropland Biomass Energy Potential from different RCPs'\n",
    "    )\n",
    "    \n",
    "    fig_marg.update_layout(\n",
    "        barmode='group',\n",
    "        xaxis_title='Time Periods',\n",
    "        yaxis_title='Biomass Energy Potential from Marginal Land <br> in PetaJoules',\n",
    "        title='Marginal Land Biomass Energy Potential from different RCPs in Time Periods'\n",
    "    )\n",
    "    \n",
    "    fig_crop.show()\n",
    "    fig_marg.show()\n",
    "    fig_total.show()\n",
    "\n",
    "    \n",
    "    return fig_crop, fig_marg , fig_total , xarrays_marg , xarrays_crop , final_potentials_marg , arrays_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da73ec16",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo_env",
   "language": "python",
   "name": "geo_env"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}