crop_modeling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0fe56655",
   "metadata": {},
   "source": [
    "# Modeling Crop Yield\n",
    "## Python modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b4f19f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "import time\n",
    "import os\n",
    "\n",
    "import dask\n",
    "from dask.distributed import Client\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "\n",
    "import geopandas\n",
    "\n",
    "import pyarrow\n",
    "from sklearn.linear_model import RidgeCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.impute import SimpleImputer\n",
    "from scipy.stats import spearmanr\n",
    "from scipy.linalg import LinAlgWarning\n",
    "from scipy.stats import pearsonr\n",
    "\n",
    "import math\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd2f0e8f",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "#### Choose a satellite.\n",
    "\n",
    "For a description of the Landsat 8 mission, see the US Geological metadata [here.]()\n",
    "\n",
    "For a description of the Sentinel 2 mission, see the US Geological metadata [here.]()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb13001c",
   "metadata": {},
   "outputs": [],
   "source": [
    "satellite = \"landsat-8-c2-l2\"\n",
    "# satellite = \"sentinel-2-l2a\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed1ab70e",
   "metadata": {},
   "source": [
    "#### Choose band combination.\n",
    "\n",
    "For a description of **Landsat** bands, see the [US Geological Survey documentation here.](https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites)\n",
    "\n",
    "For a description of **Sentinel bands**, see the [US Geological Survey documentation here.](https://www.usgs.gov/centers/eros/science/usgs-eros-archive-sentinel-2#:~:text=4%20bands%20at%2010%20meter,%2Dinfrared%20(842%20nm)\n",
    "\n",
    "According to our results, bands **(insert band selection here)** result in the best model performance for Landsat, and **(insert band selection here)** result in the best model performance for Sentinel for the task of predicting maize yields in Zambia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6155e23e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#bands = \"2-3-4\"\n",
    "# bands = \"2-3-4-8\"\n",
    "bands = \"1-2-3-4-5-6-7\"\n",
    "# bands = \"2-3-4-5-6-7-8-11-12\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a52ecb9d",
   "metadata": {},
   "source": [
    "#### Choose the number of points that were featurized.\n",
    "\n",
    "Each value in the following chunk represents the amount of thousands of points that were featurized in each respective feature file. These points represent a uniform subset of the spatial grid of Zambia. Points are spaced at uniform intervals for each selection, measured in kilometers in the longitudinal direction for each set of features. The kilometer distance interval differs for each selection below; 42,000 points results in the smallest uniform distance between points, and 4,000 points results in the greatest uniform distance between points. Selecting a greater quantity of points results in a denser spatial sample, which increases computational cost and time, but increases the spatial resolution of the model. Regardless of the quantity of points selected, each point is buffered by the same distance, resulting in a 1km^2 cell around each point.\n",
    "\n",
    "These specific options point quantities is a result of uniformly increasing the distance between points in units of kilometers prior to matching satellite images to each point. These options represent the number of points that fall within the borders of Zambia, and the numbers have been rounded to the nearest thousandth for consistency in naming files. See the [CropMOSAIKS Featurization repository](https://github.com/cropmosaiks/Featurization) for more information regarding how these distances we calculated. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a903025",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = 15\n",
    "# points = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b135301",
   "metadata": {},
   "source": [
    "#### Choose to keep only areas with crops (`True`) or to keep all points (`False`)\n",
    "\n",
    "Selecting `True` applies a \"cropland mask\" to the spatial grid of Zambia. This retains only the regions of the country in which maize is  grown, according to the **(insert source here)**. As a result, the spatial extent of the features that are fed into the model are highly subset for the specific task at hand: modeling maize yields. According to our results, selecting `True` **(insert increases or decreases here)** model performance.\n",
    "\n",
    "Selecting `False` results in modeling with the maximum spatial extent of the features, with more generalized features as a result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "def9380e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# crop_mask = True\n",
    "crop_mask = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fede5fbb",
   "metadata": {},
   "source": [
    "Choose a weighted average (`True`) or a simple mean (`False`) to use when collapsing features to administrative boundary level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fccb6992",
   "metadata": {},
   "outputs": [],
   "source": [
    "#weighted_avg = True\n",
    "weighted_avg = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aca0aad3",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Choose which months to use in the model.\n",
    "\n",
    "Note that months 10, 11, and 12 get pushed to the next year because the growing season (November - May) spans the calendar year. Maize is planted in November, starts to change color with maturity in May, and is harvested in June - August. According to our results, subsetting the months to **(insert month selection here)** increases model performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e1a5587",
   "metadata": {},
   "outputs": [],
   "source": [
    "# month_range = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]\n",
    "# month_range = [      3, 4, 5, 6, 7, 8, 9            ]\n",
    "month_range = [         4, 5, 6, 7, 8, 9            ]\n",
    "# month_range = [            5, 6, 7, 8, 9            ]\n",
    "# month_range = [         4, 5, 6, 7, 8               ]\n",
    "# month_range = [            5, 6, 7, 8               ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c15b7c5",
   "metadata": {},
   "source": [
    "#### Impute NA values by descending group levels (True) or `scikit learn`'s simple imputer (False)\n",
    "\n",
    "Imputing \"manually\" by descending group levels imputes NA values in multiple \"cascading\" steps, decreasing the proportion of inoutated values with each step. First, the NA values are imputed at by both `year` and `district`, which should yield imputed values that most closely match the feature values that would be present in the data if there was no clouds obscuring the satellite images. Next, the remaining NA values that could not be imputed by both `year` and `district` are imputed by only `district`. Lastly, the remaining NA vlaues that could not be imputed by both `year` and `district` or by just `district` are imputed by `year` only. This option gives the user more control and transparency over how the imputation is executed.\n",
    "\n",
    "Imputing using `scikit learn`'s simple imputer executes standard imputation, the details of which can be found in the `scikitlearn` documentation [here.](https://scikit-learn.org/stable/modules/generated/sklearn.impute.SimpleImputer.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9662e99",
   "metadata": {},
   "outputs": [],
   "source": [
    "impute_manual = True\n",
    "# impute_manual = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "404b7ee8-3289-467c-9dcd-84d50deb9a92",
   "metadata": {},
   "source": [
    "### Select month start\n",
    "\n",
    "If using Landsat 8 satellite imagery and want to include 2013 imagery and you have data that overlaps with the 2013 month range, select `True`. This is necessary because the Landsat 8 mission data starts in month 4. Ensure month range selected is 4 - 9."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a114312-ae97-4575-81a9-a715217a7f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "include_2013 = True\n",
    "# include_2013 = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a06291c4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Unchanging parameters\n",
    "\n",
    "The parameters in the following chunk are set for the country of Zambia for with 1000 features, regardless of the satellite selected. The start years for each satellite reflect the respective years that Landsat 8 and Sentinel 2A missions began.\n",
    "\n",
    "The number of features is set to 1000 to serve as a staple parameter among the several other parameters varied during the model optimization process. Changing this parameter in the following code chunk will result in an error because featurizing landsat imagery for a different number of features was outside the scope of this project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "026f29eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_code = \"ZMB\"\n",
    "num_features = 1000\n",
    "\n",
    "if satellite == \"landsat-8-c2-l2\":\n",
    "    year_start = 2013 # Landsat\n",
    "else:\n",
    "    year_start = 2015 # Sentinel\n",
    "year_end = 2021\n",
    "year_end_crops = 2019\n",
    "\n",
    "data_dir = \"/capstone/cropmosaiks/data\"  \n",
    "# data_dir = \"data\"  \n",
    "feature_file_name = (f'{satellite}_bands-{bands}_{country_code}_{points}k-points_{num_features}-features')\n",
    "weight_file_name = (f'{country_code}_crop_weights_{points}k-points')\n",
    "\n",
    "if points == \"4\":\n",
    "    marker_sz = 60\n",
    "elif points == \"15\":\n",
    "    marker_sz = 15\n",
    "elif points == \"24\":\n",
    "    marker_sz = 10\n",
    "else:\n",
    "    marker_sz = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aef6bb69",
   "metadata": {},
   "source": [
    "## Administrative boundaries \n",
    "\n",
    "Administrative boundaries reflect the **(insert number of districts in dataset)** district boundaries within the country of Zambia. A district can be likened to a state within the larger U.S.A. We subset the spatial grid to district level becuase the crop yield data is at the district level of specificity. The features are originally produced at higher spatial resolution, then summarized to the district level in order to train the model with ground-truth crop data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a28fd026",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_shp = geopandas.read_file(f'{data_dir}/boundaries/gadm36_{country_code}_2.shp')\n",
    "country_shp = country_shp.rename(columns = {'NAME_2': 'district'})[['district', 'geometry']]\n",
    "country_shp.district = country_shp.district.replace(\"MPongwe\", 'Mpongwe', regex=True)\n",
    "country_districts = country_shp.district.sort_values().unique().tolist()\n",
    "country_shp = country_shp.set_index('district')\n",
    "country_shp.shape\n",
    "country_shp.plot(figsize = (12,10), linewidth = 1, edgecolor = 'black' )\n",
    "# country_shp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "567ba4ec",
   "metadata": {},
   "source": [
    "## Crop yield\n",
    "\n",
    "Zambian maize yield data reflects the predicted annual maize yield provided by farmers in the month of May, when the maize matures and changes colors prior to harvest, which allows the farmers to estimate what their yield will be in the following months. These predictions are in units of metric tons per hectare and provide valuable insight to the Zambian government as they plan for the quanitites of food to import into the country in the future. For more metadata, see the websites for the [Central Statistics Office of Zambia (CSO)](https://www.zamstats.gov.zm/) and the [Summary statistics from CSO.](https://www.zamstats.gov.zm/agriculture-environment-statistics/)\n",
    "\n",
    "In order to standardize the names of all districts shared between the geoboundaries and the crop yield data, we correct for spelling, dashes, and apostrophes. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1cc9106",
   "metadata": {},
   "outputs": [],
   "source": [
    "# crop_df_full = pd.read_csv(data_dir+'/crops/cfs_maize_districts_zambia_2009_2018.csv')\n",
    "crop_df_full = pd.read_csv(data_dir+'/crops/cfs_maize_districts_zambia_2009_2022.csv')\n",
    "crop_df_full = crop_df_full[crop_df_full.year <= year_end_crops]\n",
    "crop_df_full.district = crop_df_full.district.replace(\n",
    "    {\"Itezhi-tezhi\": 'Itezhi-Tezhi',\n",
    "     \"Kapiri-Mposhi\": 'Kapiri Mposhi',\n",
    "     \"Kapiri mposhi\": 'Kapiri Mposhi',\n",
    "     \"Shang'ombo\": 'Shangombo',\n",
    "     \"Chienge\": 'Chiengi'\n",
    "    }, regex=True)\n",
    "crop_districts = crop_df_full.district.sort_values().unique().tolist()\n",
    "crop_df = crop_df_full[['district', 'year', 'yield_mt']]\n",
    "ln = len(crop_df[crop_df.year == 2016].district)\n",
    "crop_df = crop_df.set_index('district')\n",
    "ln\n",
    "# crop_df_full\n",
    "# crop_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f293d3b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(set(crop_districts) - set(country_districts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5c55cbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "list(set(country_districts) - set(crop_districts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bdc3c62",
   "metadata": {},
   "outputs": [],
   "source": [
    "country_crop = geopandas.GeoDataFrame(crop_df.join(country_shp), crs = country_shp.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ace3e2f",
   "metadata": {},
   "source": [
    "## Crop land"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "575c05e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = pd.read_feather(f\"{data_dir}/weights/{weight_file_name}.feather\")\n",
    "# weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee9d980",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights_gdf = geopandas.GeoDataFrame(\n",
    "    weights, \n",
    "    geometry = geopandas.points_from_xy(x = weights.lon, y = weights.lat), \n",
    "    crs='EPSG:4326'\n",
    ")\n",
    "weights_gdf.plot(figsize = (10,10),\n",
    "                 cmap = 'inferno',\n",
    "                 markersize = marker_sz,\n",
    "                 alpha = .9,\n",
    "                 column = 'crop_perc')\n",
    "# plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9aa0cfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights.crop_perc = weights.crop_perc.fillna(0)\n",
    "# #weights.crop_perc = weights.crop_perc + 0.0001"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1caa5a1f",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Features\n",
    "\n",
    "Append annual features files together into one file: `features_raw`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55409ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = geopandas.GeoDataFrame()\n",
    "\n",
    "for yr in range(year_start, year_end + 1):\n",
    "    print(f\"Opening: {feature_file_name}_{yr}.feather\")\n",
    "    features_x = pd.read_feather(f\"{data_dir}/features/{satellite}/{feature_file_name}_{yr}.feather\")\n",
    "\n",
    "    # concatenate the feather files together, axis = 0 specifies to stack rows (rather than adding columns)\n",
    "    features = pd.concat([features, features_x], axis=0)\n",
    "    \n",
    "    print(\"feature.shape\", features.shape)\n",
    "    print(\"Appending:\", yr)\n",
    "    print(\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "470e9d09-4c14-4adc-af1f-5291c6160f6f",
   "metadata": {},
   "source": [
    "The code below filters the `features` dataset for the first year based on `month_start`, which is dependent on satellite. These three conditions need to be met. See "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8335a372-43e8-48ca-a2dd-bfb5df5d73da",
   "metadata": {},
   "outputs": [],
   "source": [
    "if include_2013 & (satellite == \"landsat-8-c2-l2\") & (month_range == [4, 5, 6, 7, 8, 9]):\n",
    "    month_start = 4\n",
    "else:\n",
    "    month_start = 10\n",
    "    \n",
    "month_start\n",
    "\n",
    "keep = np.where(\n",
    "    ((features.year == year_start) & (features.month >= month_start)) | (features.year > year_start),\n",
    "    True, False\n",
    ")\n",
    "features = features[keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a8b4918",
   "metadata": {},
   "outputs": [],
   "source": [
    "# carry months October, November, and December over to the following year's data\n",
    "# these months represent the start of the growing season for the following year's maize yield\n",
    "features['year'] = np.where(\n",
    "    features['month'].isin([10, 11, 12]),\n",
    "    features['year'] + 1, \n",
    "    features['year'])\n",
    "\n",
    "features = features[features['year'] <= year_end]\n",
    "\n",
    "features.sort_values(['year', 'month'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8539a4a",
   "metadata": {},
   "source": [
    "### Filter month range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4b2d613",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset the features to only the month range selected at the top of the notebook\n",
    "features = features[features.month.isin(month_range)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "554ff055",
   "metadata": {},
   "source": [
    "### Pivot wider\n",
    "Here we pivot the data from long format to wide by indexing on 'lon', 'lat', 'year', 'month' and using the unstack function. We then map column names based on the month index and the associated features so month '01' is appended to each feature for that month making 0_01, 1_01 etc. This results in a Tidy data structure, with each row representing an image, and each column representing a feature for a certain month."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cac4697",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = features.set_index(['lon','lat', \"year\", 'month']).unstack()\n",
    "features.columns = features.columns.map(lambda x: '{}_{}'.format(*x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35e9b16e",
   "metadata": {},
   "source": [
    "### Replace \"inf\" values with `NaN`\n",
    "\n",
    "Infinity values are the result of **(insert reason here)**. We replace them with `NaN` because **(insert reason here)**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1845e697",
   "metadata": {},
   "outputs": [],
   "source": [
    "features.replace([np.inf, -np.inf], np.nan, inplace=True)\n",
    "features = features.reset_index()\n",
    "# features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a234a9c",
   "metadata": {},
   "source": [
    "### Attach crop weights\n",
    "Attach weight to each point (% area cropped of surrounding 1 km^2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e42e177",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = features.join(weights.set_index(['lon', 'lat']), on = ['lon', 'lat'])\n",
    "features = features.drop([\"geometry\"], axis = 1)\n",
    "# features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2be3147c",
   "metadata": {},
   "source": [
    "### Mask croppped regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9d65ec1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# any 1 km^2 cell with a crop percentage > 0 will be retained\n",
    "# the mask will not be applied if crop_mask is set to False at the top of this notebook\n",
    "if crop_mask:\n",
    "    features = features[features.crop_perc > 0]\n",
    "else:\n",
    "    pass\n",
    "# features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04c3acd3",
   "metadata": {},
   "source": [
    "### Make \"features\" a `GeoDataFrame`\n",
    "\n",
    "The coordinate reference system is set to EPSG 4326 - WGS 84, the latitude/longitude coordinate system based on the Earth's center of mass, used by the Global Positioning System."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dee1f8b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = geopandas.GeoDataFrame(\n",
    "    features, \n",
    "    geometry = geopandas.points_from_xy(x = features.lon, y = features.lat), \n",
    "    crs='EPSG:4326'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc66e85d",
   "metadata": {},
   "source": [
    "### Drop 'lat' and 'lon' columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94b47a12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop the redundant independent lon and lat columns because now that they are in a separate geometry column, there is no need for these columns\n",
    "features = features.drop(['lon', 'lat'], axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd4c25e0",
   "metadata": {},
   "source": [
    "### Join features to country geometry\n",
    "\n",
    "Join the partially processed feature data to Zambia geometry to join to district level (the highest resolution at which we have maize yield data). After imputation and some minor data porcessing steps, the features will be summarized to each district level. Although this lowers spatial resolution of the features, it is a necessary step to train the model because the district-level crop yields need to be paired with the district-level features in order to execute supervised machine learning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0fba002",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = features.sjoin(country_shp, how = 'left', predicate = 'within')\n",
    "features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3230ad6",
   "metadata": {},
   "source": [
    "### Correct column names and drop geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84705c9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop NA's from the district column (called index_right) then rename the column index_right to district\n",
    "features = (\n",
    "    features\n",
    "    # drop NA values in the `index_right` column, then rename that column to `district` for clarity\n",
    "    .dropna(subset=['index_right'])\n",
    "    .rename(columns = {\"index_right\": \"district\",})\n",
    "    .reset_index(drop = True)\n",
    ")\n",
    "# make a copy of the features, assigned to an object so we can plot the points as they are in this state\n",
    "# later, after imputation, we also copy the entire features df and save that as an object in case we want to plot or work with them in that stage \n",
    "points = features.copy()\n",
    "# save the geometries as an object to join them later to the rows of crop yield predictions\n",
    "points = features[['geometry']]\n",
    "# save the years as an object to join them later to the rows\n",
    "year = features[['year']]\n",
    "# drop geometry column for 20/21 features\n",
    "features = features.drop(['geometry'], axis = 1)\n",
    "features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2deeefc9-f501-4843-a5e3-66d02988aeea",
   "metadata": {},
   "source": [
    "### Intermediate Processing Check: Plot the Features\n",
    "\n",
    "So far, the features have been masked to just the crop regions of Zambia. Therefore, we expect the plotted points to roughly show the border of Zambia and be present throughout the country, with certain chunks of points missing in areas where crops are not grown. This serves as a \"sanity check\" that we are on-track with our spatial processing. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5753fa30-32e2-44ce-972e-e33d5d567f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "points.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca644884",
   "metadata": {},
   "source": [
    "### Impute missing values\n",
    "\n",
    "Imputing \"manually\" by descending group levels imputes NA values in multiple \"cascading\" steps, decreasing the proportion of inoutated values with each step. First, the NA values are imputed at by both `year` and `district`, which should yield imputed values that most closely match the feature values that would be present in the data if there was no clouds obscuring the satellite images. Next, the remaining NA values that could not be imputed by both `year` and `district` are imputed by only `district`. Lastly, the remaining NA vlaues that could not be imputed by both `year` and `district` or by just `district` are imputed by `year` only. This option gives the user more control and transparency over how the imputation is executed.\n",
    "\n",
    "Imputing using `scikit learn`'s simple imputer executes standard imputation, the details of which can be found in the `scikitlearn` documentation [here.](https://scikit-learn.org/stable/modules/generated/sklearn.impute.SimpleImputer.html)\n",
    "\n",
    "The imputation approach depends on the selection made at the top of this notebook for `impute_manual`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1c9cfbb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# compute the number of cells in the features dataframe, based on the amount of rows (images), months, and feature columns\n",
    "num_cells = len(features) * len(month_range) * num_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cca7eda-15bb-4d0f-8706-f8665424aaeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class bcolors:\n",
    "    BL = '\\x1b[1;34m' #GREEN\n",
    "    GR = '\\x1b[1;36m' #GREEN\n",
    "    YL = '\\x1b[1;33m' #YELLOW\n",
    "    RD = '\\x1b[1;31m' #RED\n",
    "    RESET = '\\033[0m' #RESET COLOR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1be2c7c-37e1-4bf9-80fc-d0554a14f2d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if impute_manual:\n",
    "    ln_ft = len(features)\n",
    "    ln_na = len(features.dropna())\n",
    "    print(f'Starting total row count: {bcolors.BL}{ln_ft}{bcolors.RESET}',\n",
    "          f'\\nPre-Impute NaN row count: {bcolors.RD}{ln_ft - ln_na}{bcolors.RESET}',\n",
    "          f'\\nPre-Impute NaN row %: {bcolors.RD}{((ln_ft - ln_na) / ln_ft)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\nPre-Impute NaN cell %: {bcolors.RD}{(features.isna().sum().sum() / num_cells)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\n\\nStep 1: Filling NaN values by month, year, and district group average')\n",
    "    features = (\n",
    "        features\n",
    "        .fillna(features\n",
    "                .groupby(['year', 'district'], as_index=False)\n",
    "                .transform('mean')\n",
    "               )\n",
    "    )\n",
    "    ln_ft = len(features)\n",
    "    ln_na = len(features.dropna())\n",
    "    print(f'Post step 1 NaN row count: {bcolors.YL}{ln_ft - ln_na}{bcolors.RESET}',\n",
    "          f'\\nPost step 1 NaN row %: {bcolors.YL}{((ln_ft - ln_na) / ln_ft)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\nPost step 1 NaN cell %: {bcolors.YL}{(features.isna().sum().sum() / num_cells)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\n\\nStep 2: Filling NaN values by month and district across group average')\n",
    "    features = (\n",
    "        features\n",
    "        .fillna(features\n",
    "                .groupby(['district'], as_index=False)\n",
    "                .transform('mean')\n",
    "               )\n",
    "    )\n",
    "    ln_ft = len(features)\n",
    "    ln_na = len(features.dropna())\n",
    "    print(f'Post step 2 NaN row count: {bcolors.GR}{ln_ft - ln_na}{bcolors.RESET}',\n",
    "          f'\\nPost step 2 NaN row %: {bcolors.GR}{((ln_ft - ln_na) / ln_ft)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\nPost step 2 NaN cell %: {bcolors.GR}{(features.isna().sum().sum() / num_cells)*100:.02f}{bcolors.RESET}',\n",
    "          f'\\n\\nStep 3: Drop remaining NaN values\\n')\n",
    "    features = features.dropna(axis=0)\n",
    "    print(f'Ending total row count: {bcolors.BL}{len(features)}{bcolors.RESET}')\n",
    "else:\n",
    "    features = features.set_index(['year', 'district'])\n",
    "    imputer = SimpleImputer(missing_values=np.nan, strategy='mean')\n",
    "    imputer.fit_transform(features)\n",
    "    features[:] = imputer.transform(features)\n",
    "    features = features.reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c51dd396",
   "metadata": {},
   "source": [
    "### Save copy of processed features before sumarizing training features to district level\n",
    "\n",
    "Duplicate the features dataframe at this stage so we can retain a copy of features at point resolution for all years available, which is `2013/2014/2016-2021`.\n",
    "\n",
    "    - The start year is `2016` if the satellite selected is Sentinel 2 (due to the fact that Sentinel 2 launched in June of `2015`)\n",
    "    - The start year is `2013` if the satellite selected is Landsat 8 and the month range selected was anything besides all months (due to the fact that Landsat 8 launched in February of `2013`)\n",
    "    - The start year is `2014` if the satellite selected is Landsat 8 and the month range selected was all months\n",
    "    \n",
    "This duplicated dataframe we create in the following code is called `features_all_years`. The purpose for this dataframe comes into play after the model is trained; we will be able to plug in point-resolution features from _any and all_ years from this dataframe into the trained model and observe how the model predicts crop years across space and time. It would be interesting to plot these features for each year sequentially to show how the crop prediction landscape changes by year. These point-resolution features increase the spatial resolution of the ground-truth crop data we have for the years through 2018, because our ground-truth crop data is at a the coarser  district-resolution. Furthermore, these point-resolution features are the _only_ crop data we have for the years 2020-2021. The reason we lack data from 2019 is because the Zanbia Sattistics Agency has not yet released their Crop Forecast Survey data for that year. The reason we do not have data for 2020-2021 is because Covid-19 prevented any Crop Forecast Surveys from being conducted. \n",
    "\n",
    "\n",
    "After we create the dataframe `features_all_years`, we are free to further process the original features dataframe, `features`, in order to train the model with these features and their paired ground-truth crop yields. Processing this dataframe further requires us to subset the years to the start year through the years for which we have crop data: `2013/2014/2016-2019`. This dataframe is called `features_through_2019`. The reason we subset this dataframe is because we are training the model using _supervised_ machine learning, which means we are feeding it only features that have ground-truth crop data accosicated with them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d421fce",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_all_years = features.copy()\n",
    "\n",
    "# assign the geometry column to features_2014_2021 so it can serve 2 purposes:\n",
    "# 1. plotting features sequentially by year\n",
    "# 2. the entire dataframe can be fed into the model after the model is trained on only the summarized features for 2014-2018 and the associated crop data\n",
    "features_all_years['geometry'] = points.geometry\n",
    "# moving forward in the immediate sections, summarize the `features` dataframe to district level"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf565a2-5486-43a8-9eec-9a75ba19d21a",
   "metadata": {},
   "source": [
    "### Subset the `features` dataframe to the years `2013/2014/2016-2019`\n",
    "\n",
    "As discussed above, all feature data required the previous processing steps. However, the following code is not executed on _all_ features.\n",
    "1. The `2013/2014-2019` (Landsat 8) or `2016-2019` (Sentinel 2) set of features is named `features_through_2019` so the name is agnostic to the selected satellite and month range.\n",
    "2. The features in each row (point) in `features_through_2019` will be summarized to administrative boundary level (all rows (points) have already been assigned to the district in which they reside).\n",
    "2. These summarized features in `features_through_2019` will then be joined with the district-resolution crop data. Next, we split these features into train and test sets, and then use them to train the model using ridge regression.\n",
    "3. Lastly, apply the trained model to the larger, comprehensive feature set for `2013/2014-2021` (Landsat 8) or `2016-2021` (Sentinel 2), called `features_all_years`, which includes all years regardless if those years have ground truth crop data. \n",
    "\n",
    "After we train the model, recall that we are _most interested_ in generating crop predictions for the three years `2020-2021` because the model was not trained on these years at all, and it is an \"out-of-bag sample\" in machine learning terms. The following code splits the data correctly no matter which satellite was selected at the top of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9023df42-a53c-4e92-8738-118ce91096eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# subset the features dataframe to just the years 2013/2014/2016-2019\n",
    "# this will allow us to join these years to the crop data over space and time in the next section\n",
    "features_through_2019 = features[features.year <= year_end_crops]\n",
    "\n",
    "# check that subsetting years worked\n",
    "print(features_through_2019[\"year\"].unique())\n",
    "# output should be 2016, 2017, 2018, 2019 for Sentinel 2\n",
    "# output should be 2013/2014, 2015, 2016, 2017, 2018, 2019 for Landsat 8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d7398ff",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Summarise to administrative boundary level\n",
    "Weighted by cropped area, or simple mean, depending on the selection at the top of this notebook for `weighted_avg`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "664f80ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the order of the columns in the dataframe that will be summarized and then fed into the ridge regression in order to train the model\n",
    "# we care about the order of columns specifically because in the following steps we assign only the feature columns to an object, so we need to know which 3 columns to omit by indexing\n",
    "features_through_2019.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5bb4b3-eae0-4e0e-98ee-84d61b8358a1",
   "metadata": {},
   "source": [
    "The output above show that the 3 columns that are _not_ features are the first 2 columns `year` and `district`, and the last column, `crop_perc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccdbc894-8004-4fa3-b02c-432e1acc84c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the shape to the dataframe as a sanity check\n",
    "features_through_2019.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f9f4620-f31d-4e85-afc2-5a30b942754a",
   "metadata": {},
   "source": [
    "The output above shows the number of rows and columns in the dataframe, respectively. Recall that the number of rows represents the number of points for which we have features, and the number of columns is all features for all months selected plus the columns `year`, `district`, and `crop_perc`. There are 13866 rows, meaning that is the amount of training points we have to feed into the model _before they are summarized to district level_, so this number will shrink after we summarize to district level. There are 12003 columns, which will not change after we summarize the features to district level. The number of columns that we include in the features object in the next chunk will be this number minus the 3 non-feature columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "310666d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create object that contains only feature columns, rather than all columns that would include `district`, `year`, and `crop_perc`\n",
    "# python index starts at 0, so here we specify to retain columns starting at 3 through every column besides the last column\n",
    "# the columns we omit stay in the dataframe, because we assign the selected columns to an object, but the omitted columns are not included in the calculation in the next chunk\n",
    "var_cols = features_through_2019.columns[2:-1].values.tolist()\n",
    "\n",
    "# call the object `var_cols` to check that it only includes feature columns, but do not view it in list format because it is more readable not as a list \n",
    "features_through_2019.columns[2:-1]\n",
    "# these are all the feature columns that will be fed into the `weighted_avg` calculation in the next chunk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2fcd08",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# execute the weighted average on the cropped area\n",
    "# if the object `weighted_avg` was set to True at the top of the notebook, this calculation gives more weight (importance) to area in Zambia that has a high proportion of cropland\n",
    "# this differs from the masking step that was executed earlier, because the masking step completely removed all land that was 0% cropland, and this step processes the land that is >1% cropland\n",
    "if weighted_avg:\n",
    "    features_summary = (\n",
    "        features_through_2019\n",
    "        .groupby(['year', 'district'], as_index=False)\n",
    "        .apply(lambda x: pd.Series([sum(x[v] * x.crop_perc) / sum(x.crop_perc) for v in var_cols]))\n",
    "    )\n",
    "# if weighted_avg was set to FALSE at the beginning of the notebook, the weighted avergae is not executed\n",
    "# instead, the features are simply summarized by district and year, which also occurred above if weighted_avg was set to TRUE\n",
    "else:\n",
    "    features_summary = features_through_2019.groupby(['district',\"year\"], as_index = False).mean()\n",
    "# now that we have renamed our features object to `features_summary`, note that we use this object moving forward"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcf5f1ce-57d8-46e4-bec3-c0bafe0d305e",
   "metadata": {},
   "source": [
    "Now that the features have been summarized to district and year, there are fewer rows. The dataframe we were working with before this step,  `features_through_2018`, had 13866 rows that represented points. Now we have 216 rows, as shown by the following output. Notice we still have all 12003 columns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de6e828-dd90-4ca8-b3f6-1c09096711e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_summary.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f1dc498",
   "metadata": {},
   "source": [
    "### Join crop data\n",
    "\n",
    "Combine the crop data and features into 1 dataframe in order to train the model on this object in the next section. The crop data is at district-resolution, as is the feature data in the dataframe `features_summary`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f91a28ef-1cac-4ece-84b3-03c251df1f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remind yourself of the object `crop_df` we created much earlier in the notebook when we read in the crop data in the section \"Crop yield\"\n",
    "crop_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95679bec-1884-4c2d-951e-887bff2480a5",
   "metadata": {},
   "source": [
    "This dataframe needs to undergo a little processing before we can join it to `features_summary`.\n",
    "1. Firstly, we need to subset the years to just the years for which we have features. \n",
    "2. Next, we remove 2 districts. We do this because the district shapefile that we used to assign districts to the feature data does not include the districts Mafinga or Ikelenge. It is not uncommon for district boundaries to change slightly over time in sub-Saharan Africa. This was determined by examining the unique district values in both dataframes before executing the follwoing chunk.\n",
    "3. Lastly, we reset the index so all 3 variables are columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4ad4b05-d1b3-4d9a-a4de-58de95d6c91a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take subset of crop dataframe so that it only includes the years for which we have feature data, through 2018\n",
    "# this subsets the dataframe correctly no matter which satellite was selected\n",
    "# COME BACK TO THIS PART TO DETEMRINE HOW TO HANDLE +1 \n",
    "# delete: crop_through_2018 = crop_df[crop_df.year >= year_start]\n",
    "# delete: crop_through_2018 = crop_df[crop_df.year >= year_start + 1]\n",
    "# suggested change:\n",
    "crop_through_2019 = crop_df[crop_df.year >= min(features_summary.year)]\n",
    "\n",
    "# now the crop data starts in 2013, 2014 or 2016, depending on which satellite was selected and which month range was selected\n",
    "# remove 2 districts because they are not present in one shapefile \n",
    "crop_through_2019 = crop_through_2019[~crop_through_2019.index.isin(['Mafinga', 'Ikelenge'])]\n",
    "# reset index so all 3 variables are columns\n",
    "crop_through_2019.reset_index(inplace = True)\n",
    "# take a look at the crop dataframe as a sanity check that the processing steps were correct\n",
    "# the `yield_mt` data is in units of metric tonnes per hectare of cropland planted\n",
    "crop_through_2019\n",
    "# if interested, please see more details about the crop data in our technical documentation for this project"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aaf0c93-9fc6-4c06-91a4-4e003ad3d10a",
   "metadata": {},
   "source": [
    "Now remind yourself of the columns of the features dataframe so we know which column(s) we can use to join the dataframes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a8ddc6d-7127-4e5e-987f-1fb62a306d81",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_summary.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7516d69c-9a05-4479-81d6-83829932b2c6",
   "metadata": {},
   "source": [
    "The columns we can join by are `district` and `year`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f065da58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# join the crop data with the summarized feature data through 2019\n",
    "features_summary = (\n",
    "    features_summary\n",
    "    .set_index([\"district\", \"year\"])\n",
    "    .join(other = crop_through_2019.set_index([\"district\", \"year\"]))\n",
    "    .reset_index())\n",
    "\n",
    "features_summary.columns\n",
    "# now that the data is joined, the column `yield_mt` is present in `features_summary`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab60bf3-9eae-4867-8e53-cd41dd3c8a51",
   "metadata": {},
   "source": [
    "Now that the summarized features are joined to the crop data, each row represents summarized features for a district for one year. For example, the top rows of this dataframe show the district Chadiza for several years, then the district switches to Chama, with a Chama row for each year, and so on all the way through the districts until the last one, Zambezi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2b570ee-635e-417a-b29b-1131a9b94d56",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_summary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85e3288b",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ea11764",
   "metadata": {},
   "source": [
    "### Define `x`'s and `y`'s that will be a part of training the model\n",
    "\n",
    "Since our independent variable is the features, these are the `x`'s. Our dependent variable is the crop yield in metric tonnes per hectare planted, so that will be the `y`'s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c44f3d06-6248-4c6b-bcc1-62a12b5ed0d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# this if-else statement is determined by the selection of True or False made at the top of the notebook for the object `weighted_avg`\n",
    "if weighted_avg:\n",
    "    drop_cols = ['district', 'year', 'yield_mt']\n",
    "else:\n",
    "# if the selection was False, we can drop the `crop_perc` columns because we will not use it moving forward\n",
    "    drop_cols = ['district', 'year', 'yield_mt', \"crop_perc\"]\n",
    "# note that this chunk does not actually execute the column dropping, it just defines the object that will determine which are dropped in the next chunk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcc5752e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define ALL x's (the features_summary data that will be fed into the model to train) by dropping the columns specified in the previous chunk\n",
    "# if the `axis` argument was set to 0, this code would incorrectly drop rows instead of columns\n",
    "x_all = features_summary.drop(drop_cols, axis = 1)\n",
    "\n",
    "# log transform the `yield_mt` column because that was recommended by the client\n",
    "y_all = np.log10(features_summary.yield_mt.to_numpy() + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ccecfe-5e5c-4895-8fdb-d0d165364900",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at the first 10 values of the log transformed `yield_mt` array to get an idea of the scale of the crop yields we use to train the model\n",
    "print(y_all[:10])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca312a5f",
   "metadata": {},
   "source": [
    "### Split into train and test sets\n",
    "\n",
    "This step is executed right before training the model so we can train on 80% of the data and preserve 20% for testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79eb18ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train, x_test, y_train, y_test = train_test_split(\n",
    "    x_all, y_all, test_size = 0.2, random_state = 0\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5ea5f2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of total points: \", len(x_all), \"\\n\", \n",
    "      \"Number of training points: \", len(x_train), \"\\n\",\n",
    "      \"Number of testing points: \", len(x_test), sep = \"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be2d84d0",
   "metadata": {},
   "source": [
    "### Train model using cross-validated ridge regression\n",
    "\n",
    "Please see the documentation for the function that executes this regression [here.](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RidgeCV.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf329e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ridge_cv_random = RidgeCV(cv = 5, alphas = np.logspace(-8, 8, base = 10, num = 17))\n",
    "ridge_cv_random.fit(x_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90e9d48b-28dd-47ec-9098-970958bc9b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Estimated regularization parameter: {ridge_cv_random.alpha_}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c58e251",
   "metadata": {},
   "source": [
    "### Validation set $R^2$ performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "536045c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Validation R2 performance: {ridge_cv_random.best_score_:0.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "680ac1d1",
   "metadata": {},
   "source": [
    "### Train set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9fcd436",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = np.maximum(ridge_cv_random.predict(x_train), 0)\n",
    "r2_train = r2_score(y_train, y_pred)\n",
    "\n",
    "fig, ax = plt.subplots(ncols=1)\n",
    "plt.scatter(y_pred, y_train, alpha=1, s=4)\n",
    "plt.xlabel(\"Predicted\", fontsize=15, x = .3)\n",
    "plt.ylabel(\"Ground Truth\", fontsize=15)\n",
    "plt.suptitle(r\"$\\log_{10}(1 + Crop Yield)$\", fontsize=20, y=1.02)\n",
    "plt.title((f\"Model applied to train data n = {len(x_train)}, R$^2$ = {r2_train:0.2f}\"),\n",
    "          fontsize=12, y=1.01)\n",
    "\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "\n",
    "ax.axline([0, 0], [1, 1], c = \"k\")\n",
    "\n",
    "plt.gca().spines.right.set_visible(False)\n",
    "plt.gca().spines.top.set_visible(False)\n",
    "\n",
    "\n",
    "# plt.savefig(f'images/{feature_file_name}_train_data.jpg', dpi=300)\n",
    "plt.show()\n",
    "plt.close()\n",
    "# the model is plotted with a black 45 degree line that serves as a reference of what a perfect correlation would look like\n",
    "# deviation of the line indicates that there is not a perfect correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbbc0486",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Training R^2 = {r2_train:0.2f}\\nPearsons r = {pearsonr(y_pred, y_train)[0]:0.2f}\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2550c544-4a28-4d34-841a-837223fa0bd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pearson r^2\n",
    "pearsonr(y_pred, y_train)[0] ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16c01413-8e64-4ba8-b61e-5fd8c9d10c36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# alternative way to calculate Training R^2\n",
    "ridge_cv_random.score(x_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dff83102",
   "metadata": {},
   "source": [
    "### Test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcb42c16",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "y_pred = np.maximum(ridge_cv_random.predict(x_test), 0)\n",
    "r2_test = r2_score(y_test, y_pred)\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(y_pred, y_test, alpha=1, s=4)\n",
    "plt.xlabel(\"Predicted\", fontsize=15)\n",
    "plt.ylabel(\"Ground Truth\", fontsize=15)\n",
    "plt.suptitle(r\"$\\log_{10}(1 + Crop Yield)$\", fontsize=20, y=1.02)\n",
    "plt.title(f\"Model applied to test data n = {len(x_test)}, R$^2$ = {r2_test:0.2f}\",\n",
    "          fontsize=12, y=1)\n",
    "\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "\n",
    "ax.axline([0, 0], [.75, .75], c = \"k\")\n",
    "\n",
    "plt.gca().spines.right.set_visible(False)\n",
    "plt.gca().spines.top.set_visible(False)\n",
    "\n",
    "# plt.savefig(f'images/{feature_file_name}_test_data.jpg', dpi=300)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4f6268f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Testing set R^2 = {r2_test:0.2f}\")\n",
    "print(f\"Testing set pearsons R = {pearsonr(y_pred, y_test)[0]:0.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aef4215-081f-4537-aade-bf7c3e5d3d93",
   "metadata": {},
   "source": [
    "Summary of both train and test data sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78675318-1a34-4c3d-90f8-e738abd6c2f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = np.maximum(ridge_cv_random.predict(x_all), 0)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 7))\n",
    "ax.axline([0, 0], [.75, .75], c = \"k\")\n",
    "plt.scatter(y_pred, y_all, alpha=.9, s=15)\n",
    "plt.xlabel(\"Predicted\", fontsize=15)\n",
    "plt.ylabel(\"Observed\", fontsize=15)\n",
    "plt.text(\n",
    "    0, .8, fontsize=15, fontweight=\"bold\",\n",
    "    s=f\"R$^2$={r2_train:0.2f} - Train set\",\n",
    ")\n",
    "plt.text(\n",
    "    0, .75, fontsize=15, fontweight=\"bold\",\n",
    "    s=f\"R$^2$={ridge_cv_random.best_score_:0.2f} - Validation set\",\n",
    ")\n",
    "plt.text(\n",
    "    0, .7, fontsize=15, fontweight=\"bold\",\n",
    "    s=f\"R$^2$={r2_test:0.2f} - Test set\",\n",
    ")\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "\n",
    "plt.gca().spines.right.set_visible(False)\n",
    "plt.gca().spines.top.set_visible(False)\n",
    "\n",
    "# plt.savefig(f'images/{feature_file_name}_all_data.jpg', dpi=300)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a29cb5b",
   "metadata": {},
   "source": [
    "### Use the trained model to predict crop yields over all years from 1km grid-cell resolution features \n",
    "\n",
    "Recall that after we executed imputation on all feature years in the dataframe `features`, we copied the dataframe and named it `features_all_years`. Now we can plug that into the model to visualize how our model performs over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e360ddd6-1c76-4aa0-90cf-3b6ff35b1b94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# recall the object we created earlier, before we split the features by year into those that would train the model \n",
    "# and those that would be fed into the trained model to predict crop yields\n",
    "# in years for which we do not have crop data\n",
    "features_all_years.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd87a870-2e2b-4c03-b65b-d7431fa82f73",
   "metadata": {},
   "source": [
    "In the following chunk, we drop certain columns from `features_all_years` because we only need to feed the feature data into the model to generate predictions. Using the argument `axis = 1`, we specify that we are dropping columns rather than rows. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0ad5295",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_all = features_all_years.drop([\n",
    "    'year', \n",
    "    'geometry',\n",
    "    'district',\n",
    "    'crop_perc'\n",
    "], axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd1331c-9e70-4deb-83a1-ae87cfc7e46b",
   "metadata": {},
   "source": [
    "In the following chunk, we execute the model on the features from the dataframe `features_all_years`. The crop yield predictions for each row populate a new column in the dataframe.\n",
    "\n",
    "The model is run inside the `np.maximum()` function because if we run it without being wrapped inside function, some crop predictions are negative values, but we need them all to be positive because conceptually crop yields cannot be negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bad2e993-6581-4a67-ab60-508479ad53bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_all_years['yield_prediction'] = np.maximum(ridge_cv_random.predict(x_all), 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02a4758b-e84a-47c5-a5a5-561d6fbe6823",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check out the dataframe with the new column of predictions\n",
    "features_all_years.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f250da4-4a71-483d-bca9-248b8ecd6901",
   "metadata": {},
   "source": [
    "The dataframe is already a geodataframe, so we do not have to convert it to one before mapping predictions. However, we do need to replace all the zero value crop percentage areas with `NA`. We do this by applying the `mask()` function. This function is similar to an if-else statement. If the value of the `crop_perc` is equal to 0, that value is replaced by the value of the second argument, which is `NA`. If the value of `crop_prec` is _not_ equal to zero, we retain the current value. The argument `inplace = True` executes this replacement in the same cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fe0ee44-575a-45b8-ae14-b44ce3290ac4",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_all_years['yield_prediction'].mask(features_all_years['crop_perc']==0, np.nan, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82712032-c523-4569-9c2a-26f9dea91177",
   "metadata": {},
   "source": [
    "Recall that this dataframe has a geometry column, with latitude and longitude together. In order to map the predicted features, we separate this geometry column into separate `lon` and `lat` columns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "980b20d1-e9f8-43b8-888f-32df0f261592",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the longitude and latitude from the geometry column, and make then into independent columns\n",
    "features_all_years['lon'], features_all_years['lat'] = features_all_years.geometry.x, features_all_years.geometry.y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82eb0d10-ec00-4063-a106-1b73991b9df8",
   "metadata": {},
   "source": [
    "Plot the predicted features for each year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f4364cf-5eaf-4cb4-8c77-6bd38de288ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter(x, y, c, **kwargs):\n",
    "    plt.scatter(x, y, c=c, s = 1.25)\n",
    "sns.color_palette(\"viridis\", as_cmap=True)\n",
    "g = sns.FacetGrid(\n",
    "    features_all_years, \n",
    "    col=\"year\", \n",
    "    col_wrap = 4, \n",
    "    height=5, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(scatter, \"lon\", \"lat\", \"yield_prediction\")\n",
    "g.set_axis_labels(r\"Yield Prediction\")\n",
    "# save the figure and name the file so that it represents the model parameters that created the predictions\n",
    "# plt.savefig(f'images/{feature_file_name}_all_predictions.jpg', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adf3917b-9fea-429f-9416-23851b20229e",
   "metadata": {},
   "source": [
    "Plot the model's predicted features summarized to district level. In this visualization, we choose a specific year to examine rather than visualizing all years in one figure. Visualizing the the features summarized to district level is interesting because the crop data resolution provided by Zambia Statistics Agency is at the district level, and therefore it is easier to compare our model results to those ground-truth values when they are summarized to district level as well. Furthermore, our model's crop predictions for the years 2020 and 2021 might be more valuable when summarized to district level if Zambian governments, policy-makers, farmers, and researchers wish to use this data to determine crop imports, exports, and storage according to district summaries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dceaaf7-84d7-49b4-a2da-e0086c67c314",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_all_years_summary = (\n",
    "    features_all_years\n",
    "    .groupby(['district',\"year\"], as_index = False)['yield_prediction']\n",
    "    .mean()\n",
    "    .set_index('district')\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f59be6-5dd0-49f2-bfcc-eb2a2ebea5cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# join Zambia's shapefile to the summarized features to map the districts\n",
    "# reset the index so it is a properly formatted dataframe\n",
    "features_all_years_summary = features_all_years_summary.join(country_shp).reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d88e5b61-80c4-4b22-a0fa-5a2dbd415517",
   "metadata": {},
   "source": [
    "Now that the geometries have been converted to districts from points, the geomatries are now polygons. There is still a row for each district for each year."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60f98fce-e33f-49dc-a17c-c9d52215523a",
   "metadata": {},
   "source": [
    "In order to change the year visualized, simply change the year in the following code and re-run the chunk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7084cc6b-b854-46f4-b2e2-02ba71482002",
   "metadata": {},
   "outputs": [],
   "source": [
    "features_all_years_summary[features_all_years_summary.year == 2020].plot(column = \"yield_prediction\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45acaeba-bf52-49c3-89e9-00dd307b57e6",
   "metadata": {},
   "source": [
    "Plot a boxplot for each year to visualize the range and quantile distribution of each year's crop predictions, summarized to district level. This enables us to identify years with exceptional disparities between the predicted yields by district. It also allows us to identify years that have many outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9beb5eff-a4c8-4a54-85f2-52a1fcf8e37c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "sns.boxplot(x=\"year\", y=\"yield_prediction\", data = features_all_years_summary)\n",
    "plt.xlabel(\"Year\", fontsize=15)\n",
    "plt.ylabel(\"Predicted Yield\", fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d910fe1b-bb79-497a-952a-0c7d3a3e95ee",
   "metadata": {},
   "source": [
    "Visualize the total crop yield predictions by year. This bar chart shows the sum of all the district crop yields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "427592b4-686f-4ab1-b781-5f0194dbe081",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "sns.barplot(x=\"year\", y=\"yield_prediction\", data = features_all_years_summary, estimator = sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb543c8f-3e58-4075-b3c2-95804d3c6e7a",
   "metadata": {},
   "source": [
    "## Yield and Residual Plots\n",
    "\n",
    "Create a dataframe of residuals called `residuals_df` from the `features_summary` dataframe. Note that we are _not_ using the predicted crop yields for _all_ years for these residuals, but rather the ground-truth crop yields for just the years through 2018.\n",
    "\n",
    "The residuals give us an idea of the amount of uncertianty that is present in our model. By demeaning the residuals over space, we are able to remove the uncertainty over space and better determine our model performance over time and our uncertainty over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b4918a9-c0c4-4a0d-bb32-3ff240200571",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_all = features_summary.drop(drop_cols, axis = 1)\n",
    "\n",
    "# create empty dataframe to then populate with columns\n",
    "residual_df = pd.DataFrame()\n",
    "\n",
    "residual_df[\"yield_mt\"] = features_summary.yield_mt.to_numpy()\n",
    "residual_df[\"log_yield\"] = np.log10(features_summary.yield_mt.to_numpy() + 1)\n",
    "residual_df[\"prediction\"] = np.maximum(ridge_cv_random.predict(x_all), 0)\n",
    "residual_df[\"residual\"] = residual_df[\"log_yield\"] - residual_df[\"prediction\"]\n",
    "residual_df[\"year\"] = features_summary.year\n",
    "residual_df[\"district\"] = features_summary.district\n",
    "# join the district geometries\n",
    "residual_df = residual_df.join(country_shp, how = \"left\", on = \"district\")\n",
    "\n",
    "# demean by location so we can analyze the data over time\n",
    "residual_df[\"district_yield_mean\"] = residual_df.groupby('district')['log_yield'].transform('mean')\n",
    "residual_df[\"district_prediction_mean\"] = residual_df.groupby('district')['prediction'].transform('mean')\n",
    "residual_df[\"demean_yield\"] = residual_df[\"log_yield\"] - residual_df[\"district_yield_mean\"]\n",
    "residual_df[\"demean_prediction\"] = residual_df[\"prediction\"] - residual_df[\"district_prediction_mean\"]\n",
    "residual_gdf = geopandas.GeoDataFrame(residual_df)\n",
    "\n",
    "residual_gdf.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f30a2560-4eca-4061-897f-4bb7e2f7c708",
   "metadata": {},
   "source": [
    "Visualize the residuals for the ground truth crop yields through 2018 with a boxplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1676a675-bd63-4202-b647-85f0a4fb6152",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 5))\n",
    "sns.boxplot(x=\"year\", y=\"log_yield\", data=residual_df)\n",
    "plt.xlabel(\"Year\", fontsize=15)\n",
    "plt.ylabel(\"Log Yield\", fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc5edb9d-5840-4847-8f04-80f3de8ef0e1",
   "metadata": {},
   "source": [
    "Visualize the residuals as a sum by year with a bar plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f1ffa4d-6c2d-4f16-b10d-7ba68bad0b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 5))\n",
    "sns.barplot(x=\"year\", y=\"log_yield\", data=residual_df, estimator = sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd0e82a-7842-4ffe-aa54-0998a766d162",
   "metadata": {},
   "source": [
    "Visualize the crop yield residuals by year as a histogram to determine how they are distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd79b57f-1fcc-4364-ba15-db3eed9e7e3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.histplot, \"yield_mt\", bins = 20)\n",
    "g.set_axis_labels(\"Yield (MT)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45dd3e09-071a-4bb3-b706-9843ef4e58cc",
   "metadata": {},
   "source": [
    "Visualize the log-transformed crop yield residuals by year as a histogram to compare how they are distributed after the transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51204a02-7ac4-49bc-b4a6-d4d89f18b573",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.histplot, \"log_yield\", bins = 20)\n",
    "g.set_axis_labels(r\"$\\log_{10}(1 + Crop Yield)$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8b1529f-abd6-4dde-9af2-7128cfd908ce",
   "metadata": {},
   "source": [
    "#### Crop prediction histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f611311c-5df0-4cbf-b3ac-37931e00acd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.histplot, \"prediction\", bins = 20)\n",
    "g.set_axis_labels(r\"Crop yield predictions\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1da7448e-ad5f-488d-bde5-dab420c74859",
   "metadata": {},
   "source": [
    "#### Residual histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e0628f0-fce2-47ff-b96a-8a4deec322ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.histplot, \"residual\", bins = 20)\n",
    "g.set_axis_labels(r\"Residuals\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37f5d92d-0165-409b-aa92-e8224fce2553",
   "metadata": {},
   "outputs": [],
   "source": [
    "residual_gdf.residual.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88730cab-d02a-4df4-9d05-bbbf6b42eebd",
   "metadata": {},
   "outputs": [],
   "source": [
    "residual_gdf.residual.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cbac5aa-4029-43e2-9d94-dbaa930aa87b",
   "metadata": {},
   "source": [
    "#### Log crop yield vs residuals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32ec42e2-81d5-43f4-ab32-926f5f0ffc1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.scatterplot, \"log_yield\", \"residual\")\n",
    "g.set_axis_labels(r\"$\\log_{10}(1 + Crop Yield)$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "063416fe-51cb-4505-a5c1-59864bf277b3",
   "metadata": {},
   "source": [
    "#### District residuals "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a0f7002-7619-4fc4-bcb8-98f1317cd94f",
   "metadata": {},
   "outputs": [],
   "source": [
    "if satellite == 'landsat-8-c2-l2':\n",
    "    fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 5))\n",
    "    ax1 = (residual_gdf[residual_gdf.year == 2014]\n",
    "           .plot(ax = ax1, column = \"residual\", legend = True, norm=colors.Normalize(vmin= -0.4, vmax=0.4), cmap = \"BrBG\")\n",
    "           .set_title(\"2014 Residuals\"))\n",
    "    ax2 = (residual_gdf[residual_gdf.year == 2015]\n",
    "           .plot(ax = ax2, column = \"residual\", legend = True, norm=colors.Normalize(vmin= -0.4, vmax=0.4), cmap = \"BrBG\")\n",
    "           .set_title(\"2015 Residuals\"))\n",
    "else:\n",
    "    pass\n",
    "fig, (ax1,ax2,ax3) = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))\n",
    "ax1 = (residual_gdf[residual_gdf.year == 2016]\n",
    "       .plot(ax = ax1, column = \"residual\", legend = True, norm=colors.Normalize(vmin= -0.4, vmax=0.4), cmap = \"BrBG\")\n",
    "       .set_title(\"2016 Residuals\"))\n",
    "ax2 = (residual_gdf[residual_gdf.year == 2017]\n",
    "       .plot(ax = ax2, column = \"residual\", legend = True, norm=colors.Normalize(vmin= -0.4, vmax=0.4), cmap = \"BrBG\")\n",
    "       .set_title(\"2017 Residuals\"))\n",
    "ax3 = (residual_gdf[residual_gdf.year == 2018]\n",
    "       .plot(ax = ax3, column = \"residual\", legend = True, norm=colors.Normalize(vmin= -0.4, vmax=0.4), cmap = \"BrBG\")\n",
    "       .set_title(\"2018 Residuals\"))\n",
    "\n",
    "caption = \"A positive value is an underestimated prediction (the prediction is lower than the actual yield), a negative value is an over estimated prediction\"\n",
    "plt.figtext(0.5, 0.01, caption, wrap=True, horizontalalignment='center', fontsize=12)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c06eb6d-6b28-4493-982f-c2bdc6e18517",
   "metadata": {},
   "source": [
    "#### Difference from the mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f37e53cc-050b-4718-8882-410d02c5aeb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = sns.FacetGrid(\n",
    "    residual_gdf, \n",
    "    col=\"year\", \n",
    "#     col_wrap = 3, \n",
    "    height=4, \n",
    "    aspect=1\n",
    ")\n",
    "g.map(sns.scatterplot, \"demean_yield\", \"demean_prediction\")\n",
    "g.set_axis_labels('Difference from Yield Mean', 'Difference from Prediction Mean')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c993e12-982a-42e9-8c35-28770ab2b420",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize= (6, 5))\n",
    "ax.axline([-.2, -.2], [.2, .2], c = \"k\")\n",
    "plt.scatter(residual_gdf.demean_yield, residual_gdf.demean_prediction)\n",
    "plt.title(\"Demeaned truth and predictions by district\")\n",
    "plt.xlabel('Difference from Yield Mean')\n",
    "plt.ylabel('Difference from Predictions Mean')\n",
    "r_squared = r2_score(residual_gdf[\"demean_yield\"], residual_gdf[\"demean_prediction\"])\n",
    "plt.text(\n",
    "    -0.2,\n",
    "    .18,\n",
    "    s=f\"Demeaned R$^2$ = {r_squared:0.2f}\",\n",
    "    fontsize=15,\n",
    "    fontweight=\"bold\",\n",
    ")\n",
    "plt.savefig(f'images/{feature_file_name}_demean.jpg', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "167f9228-34b2-4f04-b4b0-d47e78542e26",
   "metadata": {},
   "outputs": [],
   "source": [
    "for yr in range(year_start+1, 2018):\n",
    "    r_squared = r2_score(residual_gdf[residual_gdf.year == yr][\"demean_yield\"], residual_gdf[residual_gdf.year == yr][\"demean_prediction\"])\n",
    "    pearson_r = pearsonr(residual_gdf[residual_gdf.year == yr][\"demean_yield\"], residual_gdf[residual_gdf.year == yr][\"demean_prediction\"])\n",
    "    \n",
    "    print(yr, f\"    R^2: {r_squared:.2f}\\n\",\n",
    "          f\"Pearson's r: {pearson_r[0]:.2f}\\n\", \n",
    "          sep = \"\")\n",
    "    \n",
    "r_squared = r2_score(residual_gdf[\"demean_yield\"], residual_gdf[\"demean_prediction\"])\n",
    "pearson_r = pearsonr(residual_gdf[\"demean_yield\"], residual_gdf[\"demean_prediction\"])\n",
    "print(f\"All     R^2: {r_squared:.2f}\\n\",\n",
    "      f\"Pearson's r: {pearson_r[0]:.2f}\", sep = \"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "034bd1d1-4490-4801-a729-331c00e0a347",
   "metadata": {},
   "outputs": [],
   "source": [
    "r2 = round(pearson_r[0] ** 2, 2)\n",
    "r2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e907866-1c8b-41dc-a346-01c0a53ee3fb",
   "metadata": {},
   "source": [
    "#### Join residuals to the features for _all_ years to visualize the residuals of the features before they were summarized to district level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2f0e606-93a8-477f-b15d-b9a98d2a84ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "complete_df = (\n",
    "    features_all_years_summary\n",
    "    .set_index(['district', 'year'])\n",
    "    .join(residual_df\n",
    "          .drop('geometry', axis = 1)\n",
    "          .set_index(['district', 'year'])\n",
    "         )\n",
    "    .reset_index()\n",
    ")\n",
    "\n",
    "complete_df.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f045cfe-8479-41c7-bd04-9247843e4094",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(figsize=(10, 5))\n",
    "tidy = complete_df.melt(id_vars='year').rename(columns=str.title)\n",
    "tidy = tidy[tidy.Variable.isin(['yield_prediction', 'log_yield'])]\n",
    "sns.barplot(x='Year', y='Value', hue='Variable', data=tidy, ax=ax1, ci = None)\n",
    "sns.despine(fig)\n",
    "\n",
    "h, l = ax1.get_legend_handles_labels()\n",
    "ax1.legend(h, ['Predicted Yield', 'Observed Yield'],loc='lower left')\n",
    "\n",
    "plt.savefig(f'images/{feature_file_name}_yield_pred.jpg', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4149c18b-e563-4575-9c13-db636db64af7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "sns.barplot(x=\"year\", y=\"yield_prediction\", data=complete_df, estimator = sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38476b4e-d3b5-4f79-8d2d-651d6416900a",
   "metadata": {},
   "source": [
    "### Congratulations on completing this analysis!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}