notebooks/error_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86f58e58",
   "metadata": {},
   "outputs": [],
   "source": [
    "from branca.element import Figure\n",
    "from datetime import datetime, date\n",
    "from dotenv import load_dotenv\n",
    "from matplotlib import pyplot as plt\n",
    "from shapely.geometry import Point\n",
    "from sklearn.metrics import (\n",
    "    f1_score, accuracy_score, confusion_matrix, ConfusionMatrixDisplay, precision_score, recall_score,\n",
    "    roc_auc_score\n",
    ")\n",
    "from tqdm.notebook import tqdm\n",
    "from typing import List, Tuple\n",
    "from rasterio.plot import reshape_as_image\n",
    "\n",
    "import folium\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.patches as patches\n",
    "import os\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import seaborn as sns\n",
    "import sys\n",
    "import xarray as xr\n",
    "\n",
    "sys.path.append(\"..\")\n",
    "from src.utils import get_dvc_dir\n",
    "from src.utils import load_tif, tifs_dir, raw_dir\n",
    "from src.ETL.constants import LAT, LON\n",
    "\n",
    "from notebook_utils import get_validation_df\n",
    "\n",
    "load_dotenv();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1595429",
   "metadata": {},
   "source": [
    "# 1. Loading model and data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e29068cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = get_validation_df(\"Malawi_2021\")\n",
    "df[\"start_date\"] = pd.to_datetime(df[\"start_date\"])\n",
    "len(df[df[\"errors\"]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9301e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def metrics_from_df(df):\n",
    "    y_true = df[\"y_true\"]\n",
    "    y_pred = df[\"y_pred\"]\n",
    "    metric_functions = [f1_score, accuracy_score, precision_score, recall_score, roc_auc_score]\n",
    "    metrics = {m.__name__: round(m(y_true, y_pred), 4) for m in metric_functions}\n",
    "    return metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78bfe634",
   "metadata": {},
   "source": [
    "# 2. Plotting all predictions on map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8e570ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "map_box_token = os.environ.get(\"MAPBOX_TOKEN\")\n",
    "mapbox_kwargs = {\n",
    "    \"tiles\": 'https://api.tiles.mapbox.com/v4/mapbox.satellite/{z}/{x}/{y}.png?access_token='+map_box_token,\n",
    "    \"attr\": 'mapbox.com',\n",
    "    \"name\": \"Mapbox\"\n",
    "}\n",
    "\n",
    "def plot_map(\n",
    "    color_bboxes: List[Tuple[str, List[Tuple[float, float]]]] = None, \n",
    "    color_markers: List[Tuple[str, List[Tuple[float, float]]]] = None, \n",
    "    zoom_start: int = 15, \n",
    "    width: int = 500, \n",
    "    height: int = 200,\n",
    "    map_kwargs=mapbox_kwargs\n",
    ") -> Figure:\n",
    "    if color_bboxes:\n",
    "        # First bbox is used as center\n",
    "        first_points = color_bboxes[0][1]\n",
    "        \n",
    "    elif color_markers:\n",
    "        # First point is used as center\n",
    "        first_points = color_markers[0][1]\n",
    "        \n",
    "    else:\n",
    "        raise ValueError(\"At least one of color_bboxes or color_markers must be set\")\n",
    "        \n",
    "    center = list(np.array([p[i] for p in first_points]).mean() for i in [0,1])\n",
    "\n",
    "    m = folium.Map(location=center, zoom_start=zoom_start, **map_kwargs)\n",
    "    \n",
    "    if color_bboxes:\n",
    "        for color, bbox in color_bboxes:\n",
    "            folium.Rectangle(bounds=bbox, color=color).add_to(m)\n",
    "        \n",
    "    if color_markers:\n",
    "        feature_group = folium.FeatureGroup('markers')\n",
    "        for color, markers in color_markers:\n",
    "            for marker in markers:\n",
    "                folium.Marker(location=marker, icon=folium.Icon(color=color)).add_to(feature_group)\n",
    "        feature_group.add_to(m)\n",
    "\n",
    "    fig = Figure(width=width, height=height)\n",
    "    return fig.add_child(m)\n",
    "\n",
    "def bbox_from_corners(x_min: float, x_max: float, y_min: float, y_max: float) -> List[Tuple[float, float]]:\n",
    "    points = [(y_min, x_min), (y_min, x_max), (y_max, x_min), (y_max, x_max)]\n",
    "    return points\n",
    "\n",
    "def label_to_bbox(label_lon, label_lat, bbox_size=0.00005):\n",
    "    max_lon = label_lon + bbox_size\n",
    "    max_lat = label_lat + bbox_size\n",
    "    min_lon = label_lon - bbox_size\n",
    "    min_lat = label_lat - bbox_size\n",
    "    bbox = bbox_from_corners(min_lon, max_lon, min_lat, max_lat)\n",
    "    return bbox\n",
    "\n",
    "def instance_to_bbox(inst):\n",
    "    max_lon = inst.instance_lon + bbox_size\n",
    "    max_lat = inst.instance_lat + bbox_size\n",
    "    min_lon = inst.instance_lon - bbox_size\n",
    "    min_lat = inst.instance_lat - bbox_size\n",
    "    bbox = bbox_from_corners(min_lon, max_lon, min_lat, max_lat)\n",
    "    return bbox"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "824298d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"color\"] = df[\"errors\"].apply(lambda err: \"red\" if err else \"yellow\")\n",
    "df[\"bboxes\"] = df.apply(lambda row: label_to_bbox(row[LON], row[LAT]), axis=1)\n",
    "df[\"tif_bboxes\"] = df.apply(lambda row: label_to_bbox(row[\"instance_lon\"], row[\"instance_lat\"]), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0a116eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "color_bboxes = [(\"black\", bbox) for bbox in df[\"tif_bboxes\"]] + [(color, bbox) for color, bbox in zip(df[\"color\"], df[\"bboxes\"])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "624f8143",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_map(color_bboxes=color_bboxes, width=800, height=600, zoom_start=7, map_kwargs=mapbox_kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "353adb12",
   "metadata": {},
   "source": [
    "**Conclusion** No consistent pattern from distribution of errors. We get all the water ones right at least."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93bf6660",
   "metadata": {},
   "outputs": [],
   "source": [
    "west_df = df[df[LON] > df[LON].median()].copy()\n",
    "east_df = df[df[LON] <= df[LON].median()].copy()\n",
    "def legend_str(prefix, df): \n",
    "    crop_percent = round(df[\"y_true\"].value_counts(normalize=True)[1], 4)\n",
    "    return f\"{prefix} (Size: {len(df) }, Crop {crop_percent})\"\n",
    "\n",
    "all_dfs = {\n",
    "    legend_str(\"All\", df): df, \n",
    "    legend_str(\"East\", east_df): east_df, \n",
    "    legend_str(\"West\", west_df): west_df\n",
    "}\n",
    "metrics = {k: metrics_from_df(d) for k, d in all_dfs.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86a3c54d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(metrics).plot(kind=\"bar\", figsize=(10,10), title=f\"Validation, Size: {len(df)}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7687ca0d",
   "metadata": {},
   "source": [
    "# 4. Analyzing individal predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec76fe06",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_tif(tif_name, start_date, inst_lon, inst_lat, crop_prob, max_val=2000):\n",
    "    tif = load_tif(tifs_dir / tif_name, start_date, 30)\n",
    "    fig, axes = plt.subplots(4, 6, figsize=(30,20), sharey=\"row\")\n",
    "\n",
    "    # Create a Rectangle patch\n",
    "    x = int(np.where(tif.x==inst_lon)[0])\n",
    "    y = int(np.where(tif.y==inst_lat)[0])\n",
    "    band_values = tif.sel(x=inst_lon).sel(y=inst_lat).values\n",
    "    \n",
    "    images = []\n",
    "    if crop_prob > 0.5:\n",
    "        color = \"green\"\n",
    "    else:\n",
    "        color = \"red\"\n",
    "    for i in range(12):\n",
    "        col = i%6\n",
    "        row = (i//6)*2\n",
    "        \n",
    "        # Plot images\n",
    "        images.append(get_colors(tif.values[i], max_val))\n",
    "        ax = axes[row+1][col]\n",
    "        ax.imshow(images[i])\n",
    "        # Add the patch to the Axes\n",
    "        rect = patches.Rectangle(xy=(x, y), width=1, height=1, linewidth=1, edgecolor=\"yellow\", facecolor='none')\n",
    "        ax.add_patch(rect)\n",
    "        ax.set_title(f\"Month {i+1}\")\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        for spine in ax.spines.values():\n",
    "            spine.set_edgecolor(color)\n",
    "            spine.set_linewidth(3)\n",
    "        \n",
    "        # Plot bar charts \n",
    "        ax = axes[row][col]\n",
    "        bars = ax.bar(list(range(band_values.shape[1])), band_values[i], color=\"lightgray\")\n",
    "        bars[1].set_color(\"blue\")\n",
    "        bars[2].set_color(\"green\")\n",
    "        bars[3].set_color(\"red\")\n",
    "        ax.set_xticks([])\n",
    "        ax.set_yticks([])\n",
    "        for spine in ax.spines.values():\n",
    "            spine.set_visible(False)\n",
    "        \n",
    "    fig.tight_layout()\n",
    "    return np.array(images)\n",
    "\n",
    "def get_colors(tif_values, max_val=2000):\n",
    "    colors = (tif_values[np.array([3,2,1])] * 10000).astype(np.float64)\n",
    "\n",
    "    min_val = 0\n",
    "\n",
    "    # Enforce maximum and minimum values\n",
    "    colors[colors[:, :, :] > max_val] = max_val\n",
    "    colors[colors[:, :, :] < min_val] = min_val\n",
    "\n",
    "    for b in range(colors.shape[0]):\n",
    "        colors[b, :, :] = colors[b, :, :] * 1 / (max_val - min_val)\n",
    "\n",
    "    # rasters are in the format [bands, rows, cols] whereas images are typically [rows, cols, bands]\n",
    "    # and so our array needs to be reshaped\n",
    "    colors_reshaped = reshape_as_image(colors)\n",
    "    return colors_reshaped\n",
    "\n",
    "def gmap_url(lon, lat):\n",
    "    return f\"http://maps.google.com/maps?z=12&t=m&q=loc:{lat}+{lon}\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d97e65b",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb8c8bcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_df = pd.read_csv(raw_dir / df[\"dataset\"].iloc[0] / df[df[\"errors\"]][\"source\"].iloc[0].split(\",\")[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f906f6c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "i += 1\n",
    "error = df[df[\"errors\"]].iloc[i]\n",
    "print(i)\n",
    "images = plot_tif(\n",
    "    error.source_file,\n",
    "    error.start_date,\n",
    "    error.instance_lon, \n",
    "    error.instance_lat, \n",
    "    error.crop_probability, \n",
    "    max_val=2000)\n",
    "\n",
    "# Prints google maps url\n",
    "url = gmap_url(error[LON], error[LAT])\n",
    "print(url)\n",
    "print(f\"Crop probability: {error.crop_probability}\")\n",
    "print(f\"Crop prediction: {error.y_pred_decimal}\")\n",
    "\n",
    "# Plots square on mapbox\n",
    "#plot_map(color_bboxes=bboxes[i], width=600, height=400, zoom_start=18, map_kwargs=mapbox_kwargs)\n",
    "\n",
    "# Red = non-crop, but crop was predicted\n",
    "# Green = crop, but non-crop was predicted\n",
    "# White = tif location\n",
    "\n",
    "raw_df[(np.isclose(raw_df[\"lon\"], error[LON])) & (np.isclose(raw_df[\"lat\"], error[LAT]))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8dcc9d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(df[\"y_true\"], df[\"y_pred\"])\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[\"Non-crop\", \"Crop\"])\n",
    "disp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e12565d",
   "metadata": {},
   "source": [
    "# 5. Individual error categories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17cd9cf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "non_crops_model_got_wrong = {\n",
    "    \"lot's of vegetation\": 21,\n",
    "    \"shrubs in a field\": 35,\n",
    "    \"field\": 5,\n",
    "    \"it should be crop\": 3,\n",
    "    \"swamp\": 1,\n",
    "    \"dirt\": 3,\n",
    "    \"on field boundary\": 2,\n",
    "    \"hard to tell\": 3\n",
    "}\n",
    "\n",
    "crops_model_got_wrong = {\n",
    "    \"fallow field?\": 5,\n",
    "    \"on field boundary\": 9,\n",
    "    \"clearly it's crop\": 9,\n",
    "    \"hard to tell\": 2\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8425a12",
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'size'   : 12}\n",
    "matplotlib.rc('font', **font)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(20, 10))\n",
    "ax.bar(\n",
    "    x=list(non_crops_model_got_wrong.keys()), \n",
    "    height=list(non_crops_model_got_wrong.values()), \n",
    "    width=0.5, \n",
    "    label=\"Non-crops model got wrong\")\n",
    "ax.tick_params(axis='x', rotation=45)\n",
    "ax.set_title(\"Model Errors\")\n",
    "\n",
    "ax.bar(\n",
    "    x=list(crops_model_got_wrong.keys()), \n",
    "    height=list(crops_model_got_wrong.values()), \n",
    "    width=0.5, \n",
    "    align=\"edge\", \n",
    "    label=\"Crops model got wrong\")\n",
    "ax.legend();\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c609acd4",
   "metadata": {},
   "source": [
    "**Conclusion**: Lot's of errors from shrubs on field.\n",
    "\n",
    "Are many of these in wildlife preserves and parks?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ec5db2e",
   "metadata": {},
   "source": [
    "# 6. Are many errors inside parks?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0304f6bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_dir = get_dvc_dir(\"raw\")\n",
    "dfs = []\n",
    "for filename in [\n",
    "    \"WDPA_WDOECM_Aug2021_Public_UGA_shp_0.zip\",\n",
    "    \"WDPA_WDOECM_Aug2021_Public_UGA_shp_1.zip\",\n",
    "    \"WDPA_WDOECM_Aug2021_Public_UGA_shp_2.zip\",\n",
    "]:\n",
    "    file_path = raw_dir / \"Uganda\" / filename\n",
    "    dfs.append(gpd.read_file(file_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e665750f",
   "metadata": {},
   "outputs": [],
   "source": [
    "parks = pd.concat(dfs)\n",
    "#df = df.set_crs(4326, allow_override=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed7960f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = pd.DataFrame([{\n",
    "    \"crop_probability\": inst.crop_probability,\n",
    "    \"lat\": inst.label_lat,\n",
    "    \"lon\": inst.label_lon,\n",
    "    \"error\": inst in errors,\n",
    "    \"in_park\": parks.geometry.contains(Point(inst.label_lon,inst.label_lat)).any(),\n",
    "    \"dist_to_park\": parks.geometry.boundary.distance(Point(inst.label_lon,inst.label_lat)).min(),\n",
    "    \"mean_ndvi\": inst.labelled_array.mean(axis=0)[13],\n",
    "    \"std_ndvi\": inst.labelled_array.std(axis=0)[13]\n",
    "} for inst in instances])\n",
    "points = gpd.GeoDataFrame(points, geometry=gpd.points_from_xy(points.lon, points.lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "352a9db0",
   "metadata": {},
   "outputs": [],
   "source": [
    "points[\"y_true\"] = y_true\n",
    "points[\"y_local\"] = y_local"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30242581",
   "metadata": {},
   "outputs": [],
   "source": [
    "def set_color(row):\n",
    "    condition = row[\"in_park\"] #or row[\"dist_to_park\"] < 0.01\n",
    "    \n",
    "    if row[\"error\"] and condition:\n",
    "        return \"darkred\"\n",
    "    elif row[\"error\"]:\n",
    "        return \"lightcoral\"\n",
    "    elif condition:\n",
    "        return \"forestgreen\"\n",
    "\n",
    "    return \"lawngreen\"\n",
    "\n",
    "points[\"color\"] = points.apply(set_color, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b8e8556",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(20,20))\n",
    "\n",
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "for i, ax in enumerate(axes):\n",
    "    world[world[\"name\"] == \"Uganda\"].plot(ax=ax[0], color=\"whitesmoke\", label=\"Uganda\")\n",
    "    parks.plot(ax=ax[0], color=\"lightgreen\", label=\"Parks\")\n",
    "    if i == 0:\n",
    "        points_for_ax = points[points[\"error\"]]\n",
    "        ax[0].set_title(\"Errors in Uganda\")\n",
    "        ax[1].set_title(\"Errors in Parks\")\n",
    "    else:\n",
    "        points_for_ax = points[~points[\"error\"]]\n",
    "        ax[0].set_title(\"Correct Guesses in Uganda\")\n",
    "        ax[1].set_title(\"Correct Guess in Parks\")\n",
    "    \n",
    "    points_for_ax.plot(ax=ax[0], color=points_for_ax[\"color\"])\n",
    "    \n",
    "    for side in [\"top\", \"bottom\", \"left\", \"right\"]:\n",
    "        ax[0].spines[side].set_visible(False)\n",
    "        if side in [\"top\", \"right\"]:\n",
    "            ax[1].spines[side].set_visible(False)\n",
    "    ax[0].set_yticks([])\n",
    "    ax[0].set_xticks([])\n",
    "    \n",
    "    colors = points_for_ax[\"color\"].value_counts()\n",
    "    ax[1].set_xticklabels([\"Out of parks\", \"In parks\"])\n",
    "    ax[1].bar(x=list(colors.keys()), height=colors, color=colors.keys())\n",
    "    ax[1].set_ylim([0,350])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbf07ac9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f8b8ac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "points[\"y_local_park_filter\"] = points[\"y_local\"]\n",
    "points.loc[points[\"in_park\"], \"y_local_park_filter\"] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2f2d6e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_local_park_filter = points[\"y_local_park_filter\"].to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55868750",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total predictions changed\n",
    "def compute_change(y_new, reason):\n",
    "    total_preds_changed = len(y_local[y_local != y_new])\n",
    "    correct_changes = len(y_local[(y_local != y_new) & (y_new == y_true)])\n",
    "    incorrect_changes = len(y_local[(y_local != y_new) & (y_new != y_true)])\n",
    "    print(f\"Total predictions changed from to non-crop because of {reason}: {total_preds_changed}\")\n",
    "    print(f\"Changes that resulted in correct prediction: {correct_changes}\")\n",
    "    print(f\"Changes that resulted in incorrect prediction: {incorrect_changes}\")\n",
    "    \n",
    "    print(\"\\nF1 Score change\")\n",
    "    print(f1_score(y_true, y_local))\n",
    "    print(f1_score(y_true, y_new))\n",
    "    \n",
    "    print(\"\\nAccuracy change\")\n",
    "    print(accuracy_score(y_true, y_local))\n",
    "    print(accuracy_score(y_true, y_new))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ae491ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_change(y_local_park_filter, \"parks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90df402d",
   "metadata": {},
   "outputs": [],
   "source": [
    "crops_in_parks = points[(y_local != y_local_park_filter) & (y_local_park_filter != y_true)]\n",
    "lat = crops_in_parks.iloc[1][\"lat\"]\n",
    "lon = crops_in_parks.iloc[1][\"lon\"]\n",
    "print(gmap_url(lon, lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0f2d5e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = confusion_matrix(y_true, y_local)\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[\"Non-crop\", \"Crop\"])\n",
    "disp.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00b7d9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filtering by parks helped reduce wrong crop guesses\n",
    "cm = confusion_matrix(y_true, y_local_park_filter)\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[\"Non-crop\", \"Crop\"])\n",
    "disp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f7ad33b",
   "metadata": {},
   "source": [
    "# 7. Can we threshold NDVI to get better results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5db017",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2,2, figsize=(20,10), sharex=\"col\", sharey=\"col\")\n",
    "\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        ax = axes[i][j]\n",
    "        col = \"mean_ndvi\" if j == 0 else \"std_ndvi\"\n",
    "        if col == \"mean_ndvi\":\n",
    "            ran = (-0.2, 0.7)\n",
    "        else:\n",
    "            ran = (0.06, 0.28)\n",
    "        \n",
    "        points[points[\"y_true\"] == i][col].plot.hist(ax=ax, range=ran, bins=40, label=col)\n",
    "        points[(points[\"y_true\"] == i) &  (points[\"error\"])][col].plot.hist(ax=ax, range=ran, bins=40, label=\"Errors\")\n",
    "        \n",
    "        title =  f\"{col} - \" + (\"Crop\" if i == 1 else \"Non-crop\")\n",
    "        ax.set_title(title)\n",
    "        ax.legend()\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3539ad37",
   "metadata": {},
   "source": [
    "- Generally crops have mean NDVI > 0.3, so mean NDVI < 0.3 => non-crop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ecddac5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "points[\"y_local_ndvi_filter\"] = points[\"y_local\"]\n",
    "points.loc[points[\"mean_ndvi\"] < 0.39, \"y_local_ndvi_filter\"] = 0\n",
    "y_local_ndvi_filter = points[\"y_local_ndvi_filter\"].to_list()\n",
    "\n",
    "# Total predictions changed\n",
    "compute_change(y_local_ndvi_filter, \"NDVI\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "713d9381",
   "metadata": {},
   "source": [
    "# 8. Combining improvements from parks + NDVI "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2c9bc8b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "points[\"y_local_parks_ndvi\"] = points[\"y_local_ndvi_filter\"]\n",
    "points.loc[points[\"y_local_park_filter\"] == 0, \"y_local_parks_ndvi\"] = 0\n",
    "y_local_parks_ndvi = points[\"y_local_parks_ndvi\"].to_list()\n",
    "compute_change(y_local_parks_ndvi, \"Parks + NDVI\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73d1b97b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filtering by parks helped reduce wrong crop guesses\n",
    "cm = confusion_matrix(y_true, y_local_parks_ndvi)\n",
    "disp = ConfusionMatrixDisplay(confusion_matrix=cm, display_labels=[\"Non-crop\", \"Crop\"])\n",
    "disp.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44fb12ad",
   "metadata": {},
   "source": [
    "# 9. Using different models on the Uganda validation set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "986c0170",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_dataset_metrics = {'Global': {'Rwanda': 0.6522, 'Kenya': 0.9328, 'Togo': 0.7614, 'Uganda': 0.4331}, 'Rwanda': {'Rwanda': 0.6647, 'Kenya': 0.8549, 'Togo': 0.7153, 'Uganda': 0.3438}, 'Kenya': {'Rwanda': 0.6187, 'Kenya': 0.9474, 'Togo': 0.2745, 'Uganda': 0.4359}, 'Togo': {'Rwanda': 0.4228, 'Kenya': 0.7443, 'Togo': 0.8129, 'Uganda': 0.3375}, 'Uganda': {'Rwanda': 0.6007, 'Kenya': 0.9211, 'Togo': 0.1778, 'Uganda': 0.4211}, 'Uganda_surrounding_5': {'Rwanda': 0.6689, 'Kenya': 0.9386, 'Togo': 0.1695, 'Uganda': 0.4387}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16dedf78",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(model_dataset_metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27568cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = sns.heatmap(df, annot=True)\n",
    "p.set_xlabel(\"Model\", fontsize = 20)\n",
    "p.set_ylabel(\"Validation Set\", fontsize = 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c913f7aa",
   "metadata": {},
   "source": [
    "# 10. What if make the Uganda dataset more like the others?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4f8405",
   "metadata": {},
   "outputs": [],
   "source": [
    "points[\"y_true\"].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97b7a569",
   "metadata": {},
   "outputs": [],
   "source": [
    "# How does the score change with the balance of the dataset\n",
    "amount = np.arange(0,398,30)\n",
    "xs = []\n",
    "ys = []\n",
    "for n in amount:\n",
    "    points_balanced = pd.concat([\n",
    "         points[points[\"y_true\"] == 1],\n",
    "        points[points[\"y_true\"] == 0].sample(n=n)\n",
    "    ])\n",
    "    # If we balance our dataset to be equal crop and non crop the f1-score rises substantially\n",
    "    y = f1_score(points_balanced[\"y_local\"].to_list(), points_balanced[\"y_true\"].to_list())\n",
    "    ys.append(y)\n",
    "    xs.append(60/(60+n))\n",
    "    \n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.plot(xs, ys)\n",
    "ax.set_title(\"F1-Score vs Crop Percentage\")\n",
    "ax.set_xlabel(\"Crop Percentage\")\n",
    "ax.set_ylabel(\"F1 Score\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0388741",
   "metadata": {},
   "source": [
    "The higher the crop percentage in the dataset the better we perform."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eec62517",
   "metadata": {},
   "source": [
    " # 11. Does Uganda just have a lot of clouds?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eb4484d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from s2cloudless import S2PixelCloudDetector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d05e2532",
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_detector = S2PixelCloudDetector(\n",
    "    threshold=0.4,\n",
    "    average_over=4,\n",
    "    dilation_size=2,\n",
    "    all_bands=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfb36077",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_band_values(inst):\n",
    "    tif = load_tif_from_inst(inst)\n",
    "    band_values = tif.sel(x=inst.instance_lon).sel(y=inst.instance_lat).values\n",
    "    return band_values\n",
    "\n",
    "def get_cloud_prob(band_values):\n",
    "    #tif = load_tif_from_inst(inst)\n",
    "    #band_values = tif.sel(x=inst.instance_lon).sel(y=inst.instance_lat).values\n",
    "    bands_for_model = np.expand_dims(band_values, axis=(1,2))\n",
    "    cloud_prob = cloud_detector.get_cloud_probability_maps(pixel_for_model).flatten()\n",
    "    return cloud_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "473cb0ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_band_values = [get_band_values(inst) for inst in tqdm(instances)]\n",
    "all_band_values_np = np.array(all_band_values)\n",
    "bands_for_model = np.expand_dims(all_band_values_np, axis=2)\n",
    "bands_for_model.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27444dde",
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_prob = cloud_detector.get_cloud_probability_maps(bands_for_model)\n",
    "cloud_prob_list = list(cloud_prob[:, :, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "196be6f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "points[\"months_with_clouds_05\"] = [(probabilities > 0.5).sum() for probabilities in cloud_prob_list]\n",
    "points[\"months_with_clouds_07\"] = [(probabilities > 0.7).sum() for probabilities in cloud_prob_list]\n",
    "points[\"months_with_clouds_09\"] = [(probabilities > 0.9).sum() for probabilities in cloud_prob_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e4eb232",
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_col = \"months_with_clouds_07\"\n",
    "correct = points[~points[\"error\"]][cloud_col].value_counts()\n",
    "wrong = points[points[\"error\"]][cloud_col].value_counts()\n",
    "\n",
    "correct_crop = points[~points[\"error\"] & points[\"y_true\"]][cloud_col].value_counts()\n",
    "wrong_crop = points[points[\"error\"] & points[\"y_true\"]][cloud_col].value_counts()\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(15, 5))\n",
    "ax.set_title(cloud_col)\n",
    "ax.bar(\n",
    "    x=list(correct.keys()), \n",
    "    height=list(correct), \n",
    "    width=0.3, \n",
    "    color=\"lightgreen\",\n",
    "    label=\"Correct\")\n",
    "\n",
    "ax.bar(\n",
    "    x=list(correct_crop.keys()), \n",
    "    height=list(correct_crop), \n",
    "    width=0.3, \n",
    "    color=\"forestgreen\",\n",
    "    label=\"Correct crop\")\n",
    "\n",
    "ax.bar(\n",
    "    x=[x+0.3 for x in list(wrong.keys())], \n",
    "    height=list(wrong), \n",
    "    width=0.3, \n",
    "    color=\"lightcoral\",\n",
    "    label=\"Wrong\")\n",
    "\n",
    "ax.bar(\n",
    "    x=[x+0.3 for x in list(wrong_crop.keys())], \n",
    "    height=list(wrong_crop), \n",
    "    width=0.3, \n",
    "    color=\"red\",\n",
    "    label=\"Wrong crop\")\n",
    "\n",
    "ax.legend()\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28aa6148",
   "metadata": {},
   "source": [
    "**Conclusion**: There does not seem to be a trend between cloud cover."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}