practicals/geospatial_machine_learning.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TV85F3PEufpF"
      },
      "source": [
        "# **Let's map Africa!** Intro to Geospatial Machine Learning\n",
        "\n",
        "<img src=\"https://www.microsoft.com/en-us/research/uploads/prod/2022/01/Geospatial_Norway_header_01-2022_1920x720.jpg\" width=\"100%\" />\n",
        "\n",
        "<a href=\"https://colab.research.google.com/github/deep-learning-indaba/indaba-pracs-2023/blob/main/practicals/Indaba_2023_Prac_Template.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>\n",
        "\n",
        "© Deep Learning Indaba 2023. Apache License 2.0.\n",
        "\n",
        "**Authors:** [Akram Zaytar](https://www.linkedin.com/in/akramz/), [Gilles Q. Hacheme](https://www.linkedin.com/in/gilles-q-hacheme-a0956ab7/), [Aisha Alaagib](https://www.linkedin.com/in/aishaalaagib/), [Girmaw A. Tadesse](https://www.linkedin.com/in/girmaw-abebe-tadesse/).\n",
        "\n",
        "**Reviewers:**\n",
        "\n",
        "**Github link:**\n",
        "\n",
        "**Introduction:**\n",
        "\n",
        "In this notebook, we will introduce the field of geospatial machine learning by first going over the geospatial data primitives then solving a machine learning problem in an \"end-to-end\" fashion.\n",
        "\n",
        "We aim to cover the following:\n",
        "1. **Introduction to geospatial data**: vector and raster data primitives.\n",
        "2. **Problem framing**: introducing the problem that we are going to solve.\n",
        "3. **Approach 1**: Tabular Learning with `LightGBM`.\n",
        "    - **Data** acquisition and preprocessing: we will get the data and preprocess it for machine learning.\n",
        "    - **Model fitting**: we will fit a model to the data and conduct hyperparameter search.\n",
        "    - **Model evaluation**: we will evaluate the model on the test data.\n",
        "    - **Inference**: we will predict the output for the test data.\n",
        "4. **Approach 2**: Deep Learning with a `Sequence-to-One` model.\n",
        "\n",
        "**Topics:**\n",
        "\n",
        "Content: <font color='blue'>`Geospatial Data Analysis`</font>, <font color='blue'>`Computer Vision`</font>, <font color='blue'>`Tabular Data`</font>.\n",
        "\n",
        "Level: <font color='grey'>`Beginner`</font>, <font color='grey'>`Intermediate`</font>\n",
        "\n",
        "**Aims/Learning Objectives:**\n",
        "\n",
        "- Learn the basics of Geospatial: differentiating between the primary data types, such as vector and raster primitives. This knowledge will form the foundation for any geospatial analysis or modeling task.\n",
        "- Comprehensive Knowledge on Geospatial Machine Learning Workflow: learn how to frame a geospatial problem, acquire and preprocess relevant data, and fit a model.\n",
        "- Diversity in Modeling Approaches: By studying two distinct approaches - tabular learning with LightGBM and deep learning using a Sequence-to-One model you will appreciate the versatility of tools and techniques available in the geospatial machine learning domain, allowing to select the best approach for different types of problems.\n",
        "\n",
        "**Prerequisites:**\n",
        "\n",
        "Basic Machine Learning Concepts, Python Programming, Familiarity with Deep Learning, Hands-on Experience with Data Preprocessing.\n",
        "\n",
        "**Outline:**\n",
        "\n",
        ">[Unlocking Earth's Secrets: Intro to Geospatial Data & Machine Learning](#scrollTo=TV85F3PEufpF)\n",
        "\n",
        ">[Introduction to Geospatial Data](#scrollTo=5BBTx5ZwufpG)\n",
        "\n",
        ">>[Vector](#scrollTo=qub5T5y1ufpG)\n",
        "\n",
        ">>>[Basic Geometric Types](#scrollTo=xbrSesehufpG)\n",
        "\n",
        ">>>[Properties of Geometric Objects](#scrollTo=yE2Oz2FcufpI)\n",
        "\n",
        ">>[Raster](#scrollTo=LuQJNMJBufpJ)\n",
        "\n",
        ">[Crop Type Classification in Africa](#scrollTo=dvmQ-cJkufpJ)\n",
        "\n",
        ">>[Problem Scoping](#scrollTo=JYjJJkhLwOcT)\n",
        "\n",
        ">>>>[What type of machine learning problem is this?](#scrollTo=2qiNV7IhufpJ)\n",
        "\n",
        ">>>>[What inputs are we going to use?](#scrollTo=2zL1_uf5ufpJ)\n",
        "\n",
        ">>>>[What are we predicting?](#scrollTo=zE_EsHzzufpQ)\n",
        "\n",
        ">>>>[How will we validate the model?](#scrollTo=CGHmhNpnufpQ)\n",
        "\n",
        ">>>>[How will we measure performance?](#scrollTo=Lj_AF5bRufpQ)\n",
        "\n",
        ">>>>[Prediction](#scrollTo=Lj_AF5bRufpQ)\n",
        "\n",
        ">>>>[Target](#scrollTo=Lj_AF5bRufpQ)\n",
        "\n",
        ">>[Model 1: Tabular ML with LightGBM](#scrollTo=ejnyPSPMufpR)\n",
        "\n",
        ">>>[Data Preparation](#scrollTo=ejnyPSPMufpR)\n",
        "\n",
        ">>[Modeling](#scrollTo=KPILtkZvufpS)\n",
        "\n",
        ">>>[Frequency baseline](#scrollTo=1vob3uXJufpS)\n",
        "\n",
        ">>>[LightGBM](#scrollTo=tR3w5tP9ufpS)\n",
        "\n",
        ">>[Evaluation](#scrollTo=tiJVkLcVufpT)\n",
        "\n",
        ">>>[Investigating Class-Imbalances](#scrollTo=2XvrrCrEwOcW)\n",
        "\n",
        ">>>[XAI](#scrollTo=hKHoWqIkwOcX)\n",
        "\n",
        ">>[Inference](#scrollTo=4DJV_Hq0ufpT)\n",
        "\n",
        ">[Deep Learning](#scrollTo=A5w0WByCUpdb)\n",
        "\n",
        ">>[Data Preprocessing](#scrollTo=t60qAoCEbtm_)\n",
        "\n",
        ">>[Class Implementations](#scrollTo=xebQiqiVErI2)\n",
        "\n",
        ">>[Dataset](#scrollTo=28IXH4N4Updb)\n",
        "\n",
        ">>>[Key Methods](#scrollTo=28IXH4N4Updb)\n",
        "\n",
        ">>>[Why is this useful?](#scrollTo=28IXH4N4Updb)\n",
        "\n",
        ">>[Data Module](#scrollTo=Qs81pjALUpdb)\n",
        "\n",
        ">>>[Why is this helpful?](#scrollTo=Qs81pjALUpdb)\n",
        "\n",
        ">>[Model Architecture](#scrollTo=HXsq_Z_GUpdb)\n",
        "\n",
        ">>>[Main Components](#scrollTo=HXsq_Z_GUpdb)\n",
        "\n",
        ">>[Data Augmentation](#scrollTo=GuAZRsAgUpdb)\n",
        "\n",
        ">>>[Key Components](#scrollTo=GuAZRsAgUpdb)\n",
        "\n",
        ">>[Task](#scrollTo=oxsn6z15Amf7)\n",
        "\n",
        ">[Conclusion](#scrollTo=b2qvh17OwOcY)\n",
        "\n",
        ">[Resources](#scrollTo=Y1MbhybTwOcY)\n",
        "\n",
        ">>>[Tutorials](#scrollTo=Y1MbhybTwOcY)\n",
        "\n",
        ">>>[Geospatial Libraries](#scrollTo=Y1MbhybTwOcY)\n",
        "\n",
        ">>>[Credits](#scrollTo=Y1MbhybTwOcY)\n",
        "\n",
        "**Before you start:**\n",
        "\n",
        "For this practical, you will need to use a GPU to speed up training. To do this, go to the \"Runtime\" menu in Colab, select \"Change runtime type\" and then in the popup menu, choose \"GPU\" in the \"Hardware accelerator\" box."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XX1ufsqvfbNh"
      },
      "source": [
        "**Suggested experience level in this topic:**\n",
        "\n",
        "| Level         | Experience                            |\n",
        "| --- | --- |\n",
        "`Beginner`      | It is my first time being introduced to this work. |\n",
        "`Intermediate`  | I have done some basic courses/intros on this topic. |"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "lrXZXPJ9flmt"
      },
      "outputs": [],
      "source": [
        "# @title **Paths to follow:** What is your level of experience in the topics presented in this notebook? (Run Cell)\n",
        "experience = \"advanced\" #@param [\"beginner\", \"intermediate\", \"advanced\"]\n",
        "\n",
        "sections_to_follow=\"\"\n",
        "\n",
        "if experience == \"beginner\":\n",
        "  sections_to_follow=\"1. Introduction to Geospatial Data -> 2.1. Problem Scoping -> 2.2. Tabular ML with LightGBM -> Conclusion -> Feedback\"\n",
        "elif experience == \"intermediate\":\n",
        "  sections_to_follow=\"1. Introduction to Geospatial Data -> 2.1. Problem Scoping -> 2.2. Tabular ML with LightGBM -> 2.3. Deep Learning -> Conclusion -> Feedback\"\n",
        "else:\n",
        "  sections_to_follow=\"1. Introduction to Geospatial Data -> 2.1. Problem Scoping -> 2.2. Tabular ML with LightGBM -> 2.3. Deep Learning -> Conclusion -> Feedback\"\n",
        "\n",
        "print(f\"Based on your experience, it is advised you follow these -- {sections_to_follow} sections. Note this is just a guideline.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i977HwLZgEV9"
      },
      "source": [
        "## Installation and Imports"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "IKTm3dsvgMUK"
      },
      "outputs": [],
      "source": [
        "## Install and import anything required. Capture hides the output from the cell.\n",
        "# @title Install and import required packages. (Run Cell)\n",
        "\n",
        "%pip install rioxarray -q\n",
        "%pip install shap -q\n",
        "%pip install contextily -q\n",
        "%pip install torchgeo -q\n",
        "\n",
        "from IPython.core.debugger import set_trace\n",
        "from pathlib import Path\n",
        "import requests\n",
        "import zipfile\n",
        "import shutil\n",
        "import os\n",
        "import re\n",
        "\n",
        "import random\n",
        "import argparse\n",
        "from functools import partial\n",
        "from pathlib import Path\n",
        "from random import shuffle\n",
        "from typing import Callable, Optional, List, Dict, Any\n",
        "\n",
        "import joblib\n",
        "import kornia.augmentation as K\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import timm\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "from sklearn.preprocessing import label_binarize\n",
        "from lightning import LightningDataModule, LightningModule\n",
        "from lightning.pytorch import seed_everything, Trainer\n",
        "from lightning.pytorch.callbacks import EarlyStopping, ModelCheckpoint\n",
        "from lightning.pytorch.loggers import TensorBoardLogger\n",
        "from sklearn.model_selection import train_test_split\n",
        "from torch.optim.lr_scheduler import ReduceLROnPlateau\n",
        "from torch.utils.data import Dataset, DataLoader"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "6WjhXrfx1_1-"
      },
      "outputs": [],
      "source": [
        "#@hidden_cell# Set the URLs of the files\n",
        "# @title Next, we need to download the necessary files to be used in this practical:. (Run Cell)\n",
        "raster_url = \"https://drive.google.com/uc?export=download&id=1CrizA11Ri3jBtlMLu-58MDkM_aELt9RQ\"\n",
        "crop_url = \"https://drive.google.com/uc?export=download&id=1w1pvR0ESImXhgoCdy3QO6dV03vXbmvHH\"\n",
        "bikes_url = \"https://drive.google.com/uc?export=download&id=161vAJpEnau9pXuJEi0omDmNu38cudJq1\"\n",
        "paris_districts_url = \"https://drive.google.com/uc?export=download&id=1XyM6U-rO963zRDmtHAUx918De3qmMqzz\"\n",
        "dl_arrays_url = \"https://drive.google.com/uc?export=download&id=1vHws-qgnzA9JOO5CCPFZadKNnsckJ2Yd\"\n",
        "\n",
        "def download_file_from_google_drive(url, directory, file_name):\n",
        "    # Extract file_id from the URL\n",
        "    file_id = re.search('id=([a-zA-Z0-9_-]+)', url).group(1)\n",
        "\n",
        "    # Google Drive URL\n",
        "    base_url = \"https://drive.google.com/uc?export=download\"\n",
        "\n",
        "    session = requests.Session()\n",
        "    response = session.get(base_url, params={'id': file_id}, stream=True)\n",
        "\n",
        "    # Check the content type\n",
        "    content_type = response.headers.get('Content-Type')\n",
        "\n",
        "    # If it's an HTML page, it might be a confirmation page\n",
        "    if 'text/html' in content_type:\n",
        "        confirm_token_match = re.search('confirm=([0-9A-Za-z_]+)', response.text)\n",
        "        if confirm_token_match:\n",
        "            confirm_token = confirm_token_match.group(1)\n",
        "            response = session.get(base_url, params={'id': file_id, 'confirm': confirm_token}, stream=True)\n",
        "\n",
        "    # Save the content to the destination\n",
        "    save_response_content(response, os.path.join(directory, file_name))\n",
        "\n",
        "def save_response_content(response, destination):\n",
        "    CHUNK_SIZE = 32768\n",
        "\n",
        "    with open(destination, \"wb\") as f:\n",
        "        for chunk in response.iter_content(CHUNK_SIZE):\n",
        "            if chunk:\n",
        "                f.write(chunk)\n",
        "\n",
        "# Set the path of the directory\n",
        "data_dir = Path(\"./files\")\n",
        "if data_dir.exists(): shutil.rmtree(data_dir, ignore_errors=True)\n",
        "data_dir.mkdir(exist_ok=True)\n",
        "\n",
        "# Download\n",
        "download_file_from_google_drive(raster_url, data_dir, \"elevation.tiff\")\n",
        "download_file_from_google_drive(crop_url, data_dir, \"df.feather\")\n",
        "download_file_from_google_drive(bikes_url, data_dir, \"paris_bike_stations_mercator.gpkg\")\n",
        "download_file_from_google_drive(paris_districts_url, data_dir, \"paris_districts_utm.geojson\")\n",
        "download_file_from_google_drive(dl_arrays_url, data_dir, \"deep_learning_arrays.zip\")\n",
        "with zipfile.ZipFile(os.path.join(data_dir, \"deep_learning_arrays.zip\"), \"r\") as zip_ref:\n",
        "    zip_ref.extractall(data_dir)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5BBTx5ZwufpG"
      },
      "source": [
        "---\n",
        "\n",
        "## Introduction to Geospatial Data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3_t7B_edufpG"
      },
      "source": [
        "Geospatial data refers to **information that can be associated with geographic locations on Earth**. It comes with spatial attributes such as location and geometry shape.\n",
        "\n",
        "Examples of geospatial data:\n",
        "- Population density.\n",
        "- Weather and climate information.\n",
        "- Transportation networks.\n",
        "\n",
        "Typically, geospatial data is represented in two ways:\n",
        "\n",
        "<div style=\"text-align:center;\"> <figure> <img width=\"500px\" src=\"https://i0.wp.com/pangeography.com/wp-content/uploads/2022/05/Raster_vector_tikz.png\" /> <figcaption style=\"font-size:small;\">Image credit: <a href=\"https://pangeography.com/geographic-data-structure-vector-data-and-raster-data/\">Pan Geography</a></figcaption> </figure> </div>\n",
        "\n",
        "- **Vector**: points, lines, polygons, etc. Each vector object is a geometry that can have multiple attributes. Such data is typically saved as a vector file (`Shapefile` (.shp) and `GeoJSON`, among many).\n",
        "- **Raster**: similar to images, it is represented as a grid of cells or pixels, each cell holds a value representing a value or measurement. Raster data is typically stored in formats such as `GeoTIFF` or `NetCDF`.\n",
        "\n",
        "Next, let's get an overview of **vector** and **raster** data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qub5T5y1ufpG"
      },
      "source": [
        "### 1.1 Vector <font color='blue'>`Beginner`</font>\n",
        "\n",
        "- You'd want to use `shapely` to create and manipulate geometry objects in Python.\n",
        "- `shapely` provides support for gegraphic information systems operations, such as **spatial relationships**, **geometric operations**, and **coordinate transformations**.\n",
        "\n",
        "... but in most cases, using the higher-level library **geopandas** would be enough."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xbrSesehufpG"
      },
      "source": [
        "#### Basic Geometric Types\n",
        "\n",
        "`shapely` allows us to create basic geometric objects that are commonly used in geospatial analysis. Examples:\n",
        "- `Point`: Represents a single point in 2-3 D space.\n",
        "- `LineString`: Represents a sequence of connected points forming a line.\n",
        "- `Polygon`: Represents a filled area defined by a sequence of points that form a closed ring."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fsaws6UQufpH"
      },
      "source": [
        "Let's start by creating a `point`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3MHUhaPUXmMi"
      },
      "outputs": [],
      "source": [
        "from shapely.geometry import Point\n",
        "\n",
        "# Create a point\n",
        "point = Point(1,2)\n",
        "point"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LDLadwElufpH"
      },
      "source": [
        "... how about a polygon (enclosed list of coordinates)?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "H3oQAUL4ufpH"
      },
      "outputs": [],
      "source": [
        "from shapely.geometry import Polygon\n",
        "\n",
        "points = [(0,0), (0,2), (2,2), (2,0)]\n",
        "polygon = Polygon(points)\n",
        "polygon"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yLsZT5l7ufpH"
      },
      "source": [
        "We can also create multiple geometries:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "phjfSAzeufpH"
      },
      "outputs": [],
      "source": [
        "from shapely.geometry import MultiPoint, MultiLineString, MultiPolygon, LineString\n",
        "\n",
        "points = [Point(0, 0), Point(1, 1), Point(2, 2)]\n",
        "multipoint = MultiPoint(points)\n",
        "multipoint"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GbUtAFJeufpI"
      },
      "outputs": [],
      "source": [
        "linestrings = [LineString([(0, 0), (1, 1)]), LineString([(1, 1), (0, 2)])]\n",
        "multiline = MultiLineString(linestrings)\n",
        "multiline"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qPOr7I3_ufpI"
      },
      "outputs": [],
      "source": [
        "polygons = [Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)]), Polygon([(1, 1), (1, 2), (2, 2), (2, 1), (1, 1)])]\n",
        "multipolygon = MultiPolygon(polygons)\n",
        "multipolygon"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4QLI7xl5PUER"
      },
      "source": [
        "<font color=\"red\"> **Exercise 1:**\n",
        "Create linestring and polygon using 3 points?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PUqks--zQEuu"
      },
      "outputs": [],
      "source": [
        "#@title Exercise 1: Create linestring and polygon using 3 points?\n",
        "point1 = #update me\n",
        "point2 = #update me\n",
        "point3 = #update me\n",
        "\n",
        "# create line\n",
        "line = #update me\n",
        "line_from_tuples = #update me\n",
        "#create plygon\n",
        "poly = #update me\n",
        "poly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "lSPLEhtLPTBY"
      },
      "outputs": [],
      "source": [
        "# @title Answer to Exercise 1 (Try not to peek until you've given it a good try!')\n",
        "point1 = Point(2.2, 4.2)\n",
        "point2 = Point(7.2, -25.1)\n",
        "point3 = Point(9.26, -2.456)\n",
        "\n",
        "line = LineString([point1, point2, point3])\n",
        "line_from_tuples = LineString([(2.2, 4.2), (7.2, -25.1), (9.26, -2.456)])\n",
        "line"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yE2Oz2FcufpI"
      },
      "source": [
        "#### Properties of Geometric Objects\n",
        "\n",
        "After we create geometric objects, we can access attributes that can be useful in pratice. We highlight a few important ones:\n",
        "\n",
        "- `area`: Returns the area of a `Polygon` or `MultiPolygon` object.\n",
        "- `bounds`: Returns the bounding box of a geometric object as a tuple `(min_x, min_y, max_x, max_y)`.\n",
        "- `centroid`: Returns the geometric centroid of a geometric object as a `Point`."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XYfBiVjZufpI"
      },
      "outputs": [],
      "source": [
        "# Area of a Polygon\n",
        "polygon_area = polygon.area\n",
        "polygon_area"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "WyMPPIQiufpI"
      },
      "outputs": [],
      "source": [
        "# Bounds of a Point\n",
        "point_bounds = point.bounds\n",
        "print(f\"Point bounds: {point_bounds}\")\n",
        "print(f\"Point coordinates: {point.x, point.y}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wwHU80ysufpI"
      },
      "outputs": [],
      "source": [
        "# Centroid of a Polygon\n",
        "polygon_centroid = polygon.centroid\n",
        "polygon_centroid"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4TWhJNH2ufpI"
      },
      "source": [
        "On top of attributes, we have spatial operations and relationships that can take multiple geometries and produce new ones. We highlight the following:\n",
        "- `union`: Computes the geometric union of two objects.\n",
        "- `intersection`: Computes the geometric intersection of two objects.\n",
        "- `difference`: Computes the geometric difference of two objects.\n",
        "- `contains`: whether one geometric object contains another.\n",
        "- `intersects`: whether two geometric objects intersect."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eppoS0UNufpI"
      },
      "outputs": [],
      "source": [
        "# Create two polygon objects\n",
        "polygon1 = Polygon([(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])\n",
        "polygon2 = Polygon([(1, 1), (1, 3), (3, 3), (3, 1), (1, 1)])\n",
        "\n",
        "# Union of two polygons\n",
        "polygon_union = polygon1.union(polygon2)\n",
        "print(f\"polygon1 area: {polygon1.area}, polygon2 area: {polygon2.area}, polygon_union area: {polygon_union.area}\")\n",
        "polygon_union"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AVz0Y1x1UWHT"
      },
      "source": [
        "<font color=\"red\">**Exercise 2:**\n",
        "Computes the geometric intersection of two objects.</font> ( Hint: check ploygon_name.intersection / polygon_name.different functions)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ULVbcWrDT82z"
      },
      "outputs": [],
      "source": [
        "#@title Exercise 2: Computes the geometric intersection of two objects.\n",
        "# IMplement the Intersection of two polygons\n",
        "polygon_intersection = # update me\n",
        "print(f\"polygon1 area: {polygon1.area}, polygon2 area: {polygon2.area}, polygon_intersection area: {# update me}\")\n",
        "\n",
        "polygon_intersection"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title Answer to Exercise 2 (Try not to run until you've given it a good try!')\n",
        "\n",
        "# Intersection of two polygons\n",
        "polygon_intersection = polygon1.intersection(polygon2)\n",
        "print(f\"polygon1 area: {polygon1.area}, polygon2 area: {polygon2.area}, polygon_intersection area: {polygon_intersection.area}\")\n",
        "polygon_intersection"
      ],
      "metadata": {
        "cellView": "form",
        "id": "5XLmw6vF3Wzd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OLGmDTN2V7Q4"
      },
      "outputs": [],
      "source": [
        "#@title Implement The difference\n",
        "# update me"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "xHr77oRxufpI"
      },
      "outputs": [],
      "source": [
        "# @title Answer to Exercise 2 \"the difference\"\n",
        "polygon1.difference(polygon2)"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can also use predicates to get useful information about geospatial object relations:"
      ],
      "metadata": {
        "id": "VgxFj-gD3f_B"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9tGxdifRufpI"
      },
      "outputs": [],
      "source": [
        "# A bunch of points\n",
        "point1 = Point(0, 0)\n",
        "point2 = Point(0, 0)\n",
        "point3 = Point(1, 1)\n",
        "\n",
        "# One linestring (connected points)\n",
        "linestring = LineString([(0, 0), (1, 1), (2, 2)])\n",
        "\n",
        "# Two polygons\n",
        "polygon1 = Polygon([(0, 0), (0, 2), (2, 2), (2, 0), (0, 0)])\n",
        "polygon2 = Polygon([(1, 1), (1, 3), (3, 3), (3, 1), (1, 1)])\n",
        "\n",
        "# Spatial relationship predicates\n",
        "print(\"point1 equals point2:\", point1.equals(point2))\n",
        "print(\"point1 within polygon1:\", point1.within(polygon1))\n",
        "print(\"polygon1 intersects polygon2:\", polygon1.intersects(polygon2))\n",
        "print(\"polygon1 touches polygon2:\", polygon1.touches(polygon2))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z_C_282MufpI"
      },
      "source": [
        "In most cases, however, we'll be using `geopandas` to read vector files and analyze the data:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2e-Lv1mtufpI"
      },
      "outputs": [],
      "source": [
        "import geopandas as gpd\n",
        "\n",
        "# Load the countries dataframe using geopandas\n",
        "countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FIixrpY0ufpJ"
      },
      "source": [
        "We can use methods like `head()`, `info()`, and `describe()` to inspect and explore GeoDataFrames, similar to how you would use them with Pandas DataFrames:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "APSo1SW4ufpJ"
      },
      "outputs": [],
      "source": [
        "countries.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "C3D8_oiTufpJ"
      },
      "outputs": [],
      "source": [
        "countries.info()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eTQahTRLufpJ"
      },
      "source": [
        "Cooordinate reference systems (CRS) can take you from the geometric coordinates (numbers) to the earth's surface. `GeoPandas` allows us to inspect the CRS and reproject it if necessary:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4qjF-6FVufpJ"
      },
      "outputs": [],
      "source": [
        "# Check the CRS\n",
        "countries.crs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ff6txH2ufpJ"
      },
      "source": [
        "We can use `shapely` attributes and operations to get geometries of interest. Example:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VxB0cdG_ufpJ"
      },
      "outputs": [],
      "source": [
        "# Plot the union of all african countries\n",
        "countries[countries[\"continent\"] == \"Africa\"].unary_union"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LuQJNMJBufpJ"
      },
      "source": [
        "### 1.2 Raster <font color='blue'>`Beginner`</font>\n",
        "\n",
        "To read raster data, we need a library capable of reading geo file formats (e.g. GeoTIFF, NetCDF, etc.) and preserving their metadata.\n",
        "\n",
        "A library that can do that and is integrated with the Python ecosystem (NumPy, Geopandas, Rasterio, etc.) is `rioxarray`. `rioxarray` supports many critical tasks that we might want to do such as:\n",
        "- Reading/writing raster files.\n",
        "- Visualizing the raster data.\n",
        "- Processing the data with HPC capabilities (`NumPy` + `Dask`).\n",
        "- Geospatial operations such as...\n",
        "    - Resampling.\n",
        "    - Clipping.\n",
        "    - Reprojecting."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TmOxISJOufpJ"
      },
      "outputs": [],
      "source": [
        "import rioxarray as rxr\n",
        "\n",
        "# Read the raster Tiff file\n",
        "ds = rxr.open_rasterio(data_dir / \"elevation.tiff\")\n",
        "ds"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "flhGXEQfufpJ"
      },
      "source": [
        "We can see that the data is essentially a multi-dimensional array of values. That has the following components:\n",
        "- **Dimensions**: `band` (only elevation), `y` (latitude), `x` (longitude).\n",
        "- **Coordinates**: specify the dimension tick values on the multi-dimensional array.\n",
        "- **Attributes**: that come with the data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9Xg4yefEufpJ"
      },
      "outputs": [],
      "source": [
        "# Visualize the array\n",
        "_ = ds.plot(robust=True, cmap=\"terrain\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NrA-ASykufpJ"
      },
      "source": [
        "After this quick introduction to geospatial data, let's move to learning more about Geospatial ML by trying to solve a problem.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dvmQ-cJkufpJ"
      },
      "source": [
        "## *Crop Type Classification in Africa*\n",
        "\n",
        "In this project, we will use Geospatial machine leraning to classify farm-level crop types in Kenya using Sentinel-2 satellite imagery.\n",
        "\n",
        "<div style=\"text-align:center;\">\n",
        "    <figure>\n",
        "        <img style=\"width:50%;\" src=\"https://zindi-public-release.s3.eu-west-2.amazonaws.com/uploads/competition/image/42/header_e7a684a3-b7c0-4f53-81ca-7406f148fc5e.png\" />\n",
        "        <figcaption style=\"font-size:small;\">Question: <b>What is the farmed crop for each field?</b> (Image credit: <a href=\"https://zindi.africa/\">Zindi</a>)</figcaption>\n",
        "    </figure>\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JYjJJkhLwOcT"
      },
      "source": [
        "### 2.1 Problem Scoping <font color='blue'>`Beginner`</font>\n",
        "\n",
        "Let's answer a few fundamental questions about the problem before we begin:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2qiNV7IhufpJ"
      },
      "source": [
        "#### **What type of machine learning problem is this?**\n",
        "\n",
        "This is a supervised multiclass classification problem."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2zL1_uf5ufpJ"
      },
      "source": [
        "#### **What inputs are we going to use?**\n",
        "\n",
        "We will use pixel-level Sentinel-2 satellite imagery as input to our model:\n",
        "\n",
        "<div style=\"text-align:center;\">\n",
        "    <figure>\n",
        "        <img width=\"750px\" src=\"https://images.ctfassets.net/qfhr9fiom9gi/l1YijPOaC0jOT4pIs4FFq/c308480fc0a7ecef82c6724a4113ed7c/pasted_image_0.png\" />\n",
        "        <figcaption style=\"font-size:small;\">Sentinel-2 Bands (reference <a href=\"https://www.mdpi.com/2072-4292/8/3/166\">paper</a>)</figcaption>\n",
        "    </figure>\n",
        "</div>\n",
        "\n",
        "- The input includes **12 bands of observations from Sentinel-2 L2A**: observations in the ultra-blue, blue, green, red; visible and near-infrared (VNIR); and short wave infrared (SWIR) spectra, as well as a cloud probability layer.\n",
        "- Each pixel has measurements for **13 dates** that cover the whole farming season.\n",
        "\n",
        "Details about the bands:\n",
        "- The twelve bands are `[B01, B02, B03, B04, B05, B06, B07, B08, B8A, B09, B11, B12]`.\n",
        "    - B01 (Coastal aerosol)\n",
        "    - B02 (Blue)\n",
        "    - B03 (Green)\n",
        "    - B04 (Red)\n",
        "    - B05 (Red Edge 1)\n",
        "    - B06 (Red Edge 2)\n",
        "    - B07 (Red Edge 3)\n",
        "    - B08 (NIR - Near Infrared)\n",
        "    - B8A (Red Edge 4)\n",
        "    - B09 (Water vapor)\n",
        "    - B11 (SWIR - Shortwave Infrared 1)\n",
        "    - B12 (SWIR - Shortwave Infrared 2)\n",
        "- The cloud probability layer is a product of the Sentinel-2 atmospheric correction algorithm (Sen2Cor) and provides an estimated cloud probability (0-100%) per pixel."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SEDMlfO1UqPi"
      },
      "source": [
        "#### **Can you Frame your Approach?**\n",
        "\n",
        "<div style=\"text-align:center;\">\n",
        "    <figure>\n",
        "        <img style=\"750px;\" src=\"https://i.postimg.cc/DZLB2ytm/Indaba-bd.png\" />\n",
        "        <figcaption style=\"font-size:small;\">Approach Overview.</figcaption>\n",
        "    </figure>\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zE_EsHzzufpQ"
      },
      "source": [
        "#### **What are we predicting?**\n",
        "\n",
        "We need to **classify each farm** into one of the following categories:\n",
        "\n",
        "```\n",
        "Crop ID   Crop Type\n",
        "   1      Maize\n",
        "   2      Cassava\n",
        "   3      Common Bean\n",
        "   4      Maize & Common Bean (intercropping)\n",
        "   5      Maize & Cassava (intercropping)\n",
        "   6      Maize & Soybean (intercropping)\n",
        "   7      Cassava & Common Bean (intercropping)\n",
        "```"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CGHmhNpnufpQ"
      },
      "source": [
        "#### **How will we validate the model?**\n",
        "\n",
        "We will conduct a random train-validation split by farm IDs."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lj_AF5bRufpQ"
      },
      "source": [
        "#### **How will we measure performance?**\n",
        "\n",
        "The evaluation metric is **cross-entropy**. For each farm `field ID`, we are expected to predict the probability that the farm has a crop of that type. Example:\n",
        "\n",
        "#### Prediction\n",
        "```\n",
        "FieldID     CropId_1  CropId_2  CropId_3  CropId_4  CropId_5  CropId_6  CropId_7\n",
        "<integer>   <float>   <float>   <float>   <float>   <float>   <float>   <float>\n",
        "  1184       0.14       0.14      0.14      0.14      0.14      0.14      0.16\n",
        "```\n",
        "\n",
        "#### Target\n",
        "```\n",
        "FieldID     CropId_1  CropId_2  CropId_3  CropId_4  CropId_5  CropId_6  CropId_7\n",
        "<integer>   <float>   <float>   <float>   <float>   <float>   <float>   <float>\n",
        "  1184         0         0         1         0         0         0         0\n",
        "```\n",
        "\n",
        "---\n",
        "\n",
        "Next, we want to prepare the dataset for machine learning."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1jdvHQQ5bTDU"
      },
      "source": [
        "### 2.2 Model 1: Tabular ML with `LightGBM` <font color='orange'>`Intermediate`</font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ejnyPSPMufpR"
      },
      "source": [
        "#### Data Preparation\n",
        "\n",
        "In this section, we want to do the following:\n",
        "\n",
        "1. Remove the pixels where the cloud probability value is greater than `50%`\n",
        "2. Split the data into train/validation/test.\n",
        "3. Verify that no data leakage is present in the train/validation/test data.\n",
        "4. Check the distribution of each channel or band.\n",
        "5. Plot the farms by their labels in a map.\n",
        "6. Visualize a single farm's NDVI as it changes through time (13 dates)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ROgs27VNufpR"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "\n",
        "df = pd.read_feather(data_dir / \"df.feather\")\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#@title plot the data in the context of Kenya\n",
        "import warnings\n",
        "import geopandas as gpd\n",
        "from shapely.geometry import Point\n",
        "\n",
        "# Sample one pixel per field to simplify visualization\n",
        "report = df.copy()\n",
        "report = report.sample(frac=1)\n",
        "report = report[[\"field\", \"lat\", \"lon\"]].drop_duplicates()\n",
        "report = gpd.GeoDataFrame(report, geometry=[Point(xy) for xy in zip(report['lon'], report['lat'])])\n",
        "report = report[[\"field\", \"geometry\"]].drop_duplicates(subset=\"field\")\n",
        "\n",
        "# Get Kenya\n",
        "with warnings.catch_warnings():\n",
        "  warnings.filterwarnings(\"ignore\", category=FutureWarning)\n",
        "  world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
        "kenya = world[world['name'] == 'Kenya']\n",
        "\n",
        "# Plot Kenya and our fields\n",
        "fig, ax = plt.subplots()\n",
        "_ = kenya.boundary.plot(ax=ax, color=\"blue\")\n",
        "_ = report.plot(ax=ax, color=\"red\")\n",
        "ax.set_title(\"Kenya (blue). Fields (red)\")\n",
        "plt.show()"
      ],
      "metadata": {
        "cellView": "form",
        "id": "wlQweiC3NXMt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#@title plot the data globally\n",
        "\n",
        "# Plot country boundaries and our fields\n",
        "fig, ax = plt.subplots(figsize=(10, 7))\n",
        "_ = world.boundary.plot(ax=ax, color=\"blue\")\n",
        "_ = report.plot(ax=ax, color=\"red\")\n",
        "ax.set_title(\"Fields (Red)\")\n",
        "plt.show()"
      ],
      "metadata": {
        "cellView": "form",
        "id": "Q7Vt_rbLM0Fj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "37YL6qx0ufpR"
      },
      "source": [
        "Each `(pixel, time)` is a row. Let's start by removing the pixels that are cloudy:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "VHQsqbjrufpR"
      },
      "outputs": [],
      "source": [
        "# Drop pixels that have a cloud cover greater than 50\n",
        "df = df[df[\"CLD\"] < 50]\n",
        "\n",
        "# No need to keep the `CLD` column anymore\n",
        "df = df.drop(columns=[\"CLD\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FZYQqJoQufpR"
      },
      "source": [
        "Let's quickly check if we have missing values:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b2hB-iYYZh5F"
      },
      "source": [
        "<font color=\"red\">Exercise 3 update the code to check the missing values? </font> (*Hint: make sure you convert your df to numpy array to check if you have any missing values*)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KZQecee0Y1p3"
      },
      "outputs": [],
      "source": [
        "#@title Exercise 3\n",
        "import numpy as np\n",
        "\n",
        "answer = # update me\n",
        "\n",
        "assert answer == False\n",
        "answer"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "mmRMY0PrufpR"
      },
      "outputs": [],
      "source": [
        "# @title Answer to Exercise 3 (Try not to peek until you've given it a good try!')\n",
        "answer = np.any(df.isna())\n",
        "assert answer == False\n",
        "answer"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "k_2IbyajufpR"
      },
      "outputs": [],
      "source": [
        "#@title Split the data to train, val and test using the field column\n",
        "# Set the seed for reproducibility\n",
        "import numpy as np\n",
        "np.random.seed(42)\n",
        "\n",
        "# NOTE: the `deploy` labels are hidden and marked with `crop == 0`\n",
        "# We use the label `crop == 0` to get the `deploy` frame and ignore it\n",
        "deploy_crop_id = 0\n",
        "deploy = df[df[\"crop\"] == deploy_crop_id]\n",
        "\n",
        "# Train/Validation/Test are the remaining rows\n",
        "train_val_test = df[~df[\"field\"].isin(deploy[\"field\"])]\n",
        "\n",
        "# Get the unique field IDs from the train/validation/Test rows\n",
        "train_val_test_field_ids = train_val_test[\"field\"].sample(frac=1).unique()\n",
        "\n",
        "# Randomly select 80/10/10 split for train/val/test\n",
        "val_field_ids = np.random.choice(train_val_test_field_ids, size=int(len(train_val_test_field_ids) * 0.1), replace=False)\n",
        "test_field_ids = np.random.choice(list(set(train_val_test_field_ids) - set(val_field_ids)), size=int(len(train_val_test_field_ids) * 0.1), replace=False)\n",
        "train_field_ids = list(set(train_val_test_field_ids) - set(val_field_ids) - set(test_field_ids))\n",
        "\n",
        "# Create `train`, `val`, and `test` sets based on the validation field IDs\n",
        "train = train_val_test[train_val_test[\"field\"].isin(train_field_ids)]\n",
        "val = train_val_test[train_val_test[\"field\"].isin(val_field_ids)]\n",
        "test = train_val_test[train_val_test[\"field\"].isin(test_field_ids)]\n",
        "\n",
        "# print the shapes of the train/val/test sets\n",
        "train.shape, val.shape, test.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZP1aqqDoufpR"
      },
      "source": [
        "Let's verify that no data leakage is happening. We define leakage as follows:\n",
        "\n",
        "> A validation or test farm pixels in the training dataframe (or the reverse)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ViklB6HOufpR"
      },
      "outputs": [],
      "source": [
        "# Verify that the sets of field IDs from `train`, `val`, and `test` are mutually exclusive\n",
        "assert len(set(train[\"field\"].unique()).intersection(set(val[\"field\"].unique()))) == 0\n",
        "assert len(set(train[\"field\"].unique()).intersection(set(test[\"field\"].unique()))) == 0\n",
        "assert len(set(val[\"field\"].unique()).intersection(set(test[\"field\"].unique()))) == 0"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TI4sK_ssufpR"
      },
      "source": [
        "Next, let's check the distribution of the band values we have:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FPSA7Wk-ufpR"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "\n",
        "g = sns.displot(data=train.drop(columns=[\"time\", \"lat\", \"lon\", \"field\", \"crop\"]).melt().sample(100_000), x='value', hue=\"band\", multiple='stack')\n",
        "plt.title('Distribution of Input for each Band')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "e4AKUnrkufpR"
      },
      "source": [
        "Let's visualize the spatial distribution of the farms by their crop IDs:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JTj4rmNzufpR"
      },
      "outputs": [],
      "source": [
        "import contextily as ctx\n",
        "from shapely.geometry import Polygon\n",
        "\n",
        "# Create a new GeoDataFrame\n",
        "d = train[[\"field\", \"lon\", \"lat\", \"crop\"]].copy()\n",
        "\n",
        "# Map crop IDs to names\n",
        "id_to_name = {\n",
        "    1: 'Maize',\n",
        "    2: 'Cassava',\n",
        "    3: 'Common Bean',\n",
        "    4: 'Maize & Common Bean (intercropping)',\n",
        "    5: 'Maize & Cassava (intercropping)',\n",
        "    6: 'Maize & Soybean (intercropping)',\n",
        "    7: 'Cassava & Common Bean (intercropping)',\n",
        "}\n",
        "\n",
        "# Replace the 'crop' column with mapped names\n",
        "d['crop'] = d['crop'].map(id_to_name)\n",
        "\n",
        "# Group by field and crop, and create polygons from point coordinates\n",
        "polygons = d.groupby(['field', 'crop']).apply(lambda df: Polygon(zip(df.lon, df.lat))).reset_index()\n",
        "polygons.columns = ['field', 'crop', 'geometry']\n",
        "\n",
        "# Create a GeoDataFrame\n",
        "gdf = gpd.GeoDataFrame(polygons, geometry='geometry', crs=\"EPSG:4326\")\n",
        "\n",
        "# Create a figure and an axes object\n",
        "fig, ax = plt.subplots(figsize=(10, 10))\n",
        "\n",
        "# Plot the GeoDataFrame using the 'crop' column to color the polygons\n",
        "gdf.plot(column=\"crop\", legend=True, ax=ax, cmap=\"Accent\")\n",
        "\n",
        "# TODO: comment to check the spatial coverage of all training labels\n",
        "ax.set_xlim([34.2, 34.3])\n",
        "ax.set_ylim([.5, .6])\n",
        "\n",
        "# Add a basemap\n",
        "ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.Stamen.Terrain)\n",
        "\n",
        "# Show the plot\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wZkQtB2WufpR"
      },
      "source": [
        "Every field has a history of 13 dates across the growing farming season. We one field's NDVI evolution.\n",
        "\n",
        "To emphasize vegetation, a common technique in remote sensing is to use the Normalized Difference Vegetation Index (NDVI). NDVI is a measure of the amount and condition of green vegetation present. The NDVI is calculated from the visible and near-infrared light reflected by vegetation. The formula for NDVI is:\n",
        "\n",
        "$$NDVI = \\frac{NIR-Red}{NIR+Red}$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "J7ImEds1ufpR"
      },
      "outputs": [],
      "source": [
        "from random import choice\n",
        "\n",
        "d = train.copy()\n",
        "field_id = choice(d[\"field\"].unique().tolist())\n",
        "d = d.loc[d[\"field\"] == field_id]\n",
        "d = d[[\"time\", \"lat\", \"lon\", \"B08\", \"B04\"]]\n",
        "d = d.melt(id_vars=[\"time\", \"lat\", \"lon\"], value_vars=[\"B08\", \"B04\"], var_name=\"band\", value_name=\"value\")\n",
        "d = d.set_index([\"time\", \"lat\", \"lon\", \"band\"])\n",
        "d = d.to_xarray()\n",
        "\n",
        "# Calculate NDVI and assign to 'value'\n",
        "d = (d.sel(band='B08') - d.sel(band='B04')) / (d.sel(band='B08') + d.sel(band='B04'))\n",
        "\n",
        "# Normalize NDVI to the range 0-255\n",
        "d = (d.apply(lambda x: (x - x.min()) / (x.max() - x.min())) * 255).astype(\"uint8\")\n",
        "\n",
        "# Plot the field (each column should represent a time)\n",
        "_ = d[\"value\"].plot.imshow(col=\"time\", x=\"lon\", y=\"lat\", col_wrap=5, figsize=(13, 7))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UA2oDNHtufpR"
      },
      "source": [
        "There are many more things that we can explore with data. For now, let's skip ahead to the **modeling** section.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KPILtkZvufpS"
      },
      "source": [
        "#### Modeling\n",
        "\n",
        "In this section, we aim to train a `LightGBM` model to predict each farm's crop type by summarizing the historical band information. We will go over the following:\n",
        "\n",
        "<div style=\"text-align:center;\">\n",
        "    <figure>\n",
        "        <img style=\"width:66%;\" src=\"https://i.postimg.cc/8z3qzdJ0/Screenshot-2023-08-12-at-20-15-41.png\" />\n",
        "        <figcaption style=\"font-size:small;\">Processing steps</figcaption>\n",
        "    </figure>\n",
        "</div>\n",
        "\n",
        "- Establishing the validation metric of a **frequency based model** that always predicts crop type frequencies derived from `y_train`.\n",
        "- Feature engineering: we will calculate the following S2-based indidces:\n",
        "$$\n",
        "\\begin{align*}\n",
        "\\text{NDVI} & = \\frac{{B08 - B04}}{{B08 + B04}} \\\\\n",
        "\\text{RDNDVI1} & = \\frac{{B08 - B05}}{{B08 + B05}} \\\\\n",
        "\\text{RDNDVI2} & = \\frac{{B08 - B06}}{{B08 + B06}} \\\\\n",
        "\\text{GCVI} & = \\frac{{B08}}{{B03}} - 1 \\\\\n",
        "\\text{RDGCVI1} & = \\frac{{B08}}{{B05}} - 1 \\\\\n",
        "\\text{RDGCVI2} & = \\frac{{B08}}{{B06}} - 1 \\\\\n",
        "\\text{MTCI} & = \\frac{{B08 - B05}}{{B05 - B04}} \\\\\n",
        "\\text{MTCI2} & = \\frac{{B06 - B05}}{{B05 - B04}} \\\\\n",
        "\\text{REIP} & = 700 + 40 \\left( \\frac{{(B04 + B07)/2 - B05}}{{B07 - B05}} \\right) \\\\\n",
        "\\text{NBR1} & = \\frac{{B08 - B11}}{{B08 + B11}} \\\\\n",
        "\\text{NBR2} & = \\frac{{B08 - B12}}{{B08 + B12}} \\\\\n",
        "\\text{NDTI} & = \\frac{{B11 - B12}}{{B11 + B12}} \\\\\n",
        "\\text{CRC} & = \\frac{{B11 - B03}}{{B11 + B03}} \\\\\n",
        "\\text{STI} & = \\frac{{B11}}{{B12}}\n",
        "\\end{align*}\n",
        "$$\n",
        "- **Spatial median-aggregation** by field `ID` and `time`.\n",
        "- Conduct **period-based temporal aggregation** and for each band and index, create period-based columns using the following temporal groups:\n",
        "    - `period 1`\n",
        "        - *2019-06-06*\n",
        "    - `period 2`\n",
        "        - *2019-07-01*\n",
        "        - *2019-07-06*\n",
        "        - *2019-07-11*\n",
        "        - *2019-07-21*\n",
        "    - `period 3`\n",
        "        - *2019-08-05*\n",
        "        - *2019-08-15*\n",
        "        - *2019-08-25*\n",
        "    - `period 4`\n",
        "        - *2019-09-09*\n",
        "        - *2019-09-19*\n",
        "        - *2019-09-24*\n",
        "        - *2019-10-04*\n",
        "    - `period 5`\n",
        "        - *2019-11-03*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1vob3uXJufpS"
      },
      "source": [
        "#### Frequency baseline\n",
        "\n",
        "A simple algorithm can:\n",
        "\n",
        "1. Calculate the frequency of each crop type from the training data.\n",
        "2. Predicts the same frequencies for each validation field.\n",
        "\n",
        "The resulting metric serves to filter out any subsequent models."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KKJiRoqQufpS"
      },
      "outputs": [],
      "source": [
        "def prepare_Xy(df):\n",
        "    d = df.copy()\n",
        "    d = d.groupby([\"field\", \"time\"], as_index=False).mean()\n",
        "    d = d.drop(\"time\", axis=1).groupby(\"field\", as_index=False).mean()\n",
        "    return d.drop(['field', 'crop'], axis=1), d['crop']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Zq3TjEFXufpS"
      },
      "outputs": [],
      "source": [
        "X_train, y_train = prepare_Xy(train)\n",
        "X_val, y_val = prepare_Xy(val)\n",
        "X_train.shape, y_train.shape, X_val.shape, y_val.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3GZkD3ZBufpS"
      },
      "outputs": [],
      "source": [
        "# Calculate the class frequencies from `y_train` in order to generate the baseline predictions\n",
        "y_val_hat = np.repeat(y_train.value_counts(normalize=True).sort_index().values[None,...], y_val.shape[0], axis=0)\n",
        "y_val_hat.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eWcYrPG9pa-c"
      },
      "source": [
        "<font color=\"red\">Exercise 4 Calculate the cross entropy loss? </font> (*Hint: import the cross entropy loss from sklearn*)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "9rueQOjUpLig"
      },
      "outputs": [],
      "source": [
        "#@title Exercise 4: Calculate the cross entrop loss\n",
        "\n",
        "# Calculate cross-entropy\n",
        "cross_entropy = ...\n",
        "\n",
        "assert cross_entropy == 1.4435641708151055\n",
        "print(f'Cross-entropy is {cross_entropy}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "yQMpP8zAufpS"
      },
      "outputs": [],
      "source": [
        "#@title Answer to Exercise 4 (Try not to peek until you've given it a good try!')\n",
        "from sklearn.metrics import log_loss\n",
        "\n",
        "# Calculate cross-entropy\n",
        "cross_entropy = log_loss(y_val, y_val_hat)\n",
        "\n",
        "assert cross_entropy == 1.4435641708151055\n",
        "\n",
        "print(f'Cross-entropy is {cross_entropy}')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9hLqBAitufpS"
      },
      "outputs": [],
      "source": [
        "# .. to be used for comparison\n",
        "baseline_ce = cross_entropy\n",
        "baseline_ce"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TKCsZnxkufpS"
      },
      "source": [
        "Any model that we construct should have a validation cross-entropy less than the baseline cross-entropy."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tR3w5tP9ufpS"
      },
      "source": [
        "#### `LightGBM`\n",
        "\n",
        "We will create functions that cover the data preparation steps in the original section description.\n",
        "\n",
        "Let's implement the feature engineering function that would add additional vegetation indices of interest:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "0Fc_Fbk_ufpS"
      },
      "outputs": [],
      "source": [
        "# @title Define feature engineering function `calculate_indices(df)` (Run Cell)\n",
        "def calculate_indices(df):\n",
        "    \"\"\"\n",
        "    Compute various spectral indices commonly used in remote sensing for vegetation monitoring.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    df : pandas.DataFrame\n",
        "        Input DataFrame containing columns for the different band values.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    pandas.DataFrame\n",
        "        The DataFrame with added columns for the calculated indices.\n",
        "    \"\"\"\n",
        "    # Make a copy of the dataframe to avoid SettingWithCopyWarning\n",
        "    df = df.copy()\n",
        "\n",
        "    # Normalized Difference Vegetation Index (NDVI)\n",
        "    df['NDVI'] = (df['B08'] - df['B04']) / (df['B08'] + df['B04'])\n",
        "\n",
        "    # Red-edge Normalized Difference Vegetation Index (RDNDVI)\n",
        "    df['RDNDVI1'] = (df['B08'] - df['B05']) / (df['B08'] + df['B05'])\n",
        "    df['RDNDVI2'] = (df['B08'] - df['B06']) / (df['B08'] + df['B06'])\n",
        "\n",
        "    # Green Chlorophyll Vegetation Index (GCVI)\n",
        "    df['GCVI'] = df['B08'] / df['B03'] - 1\n",
        "\n",
        "    # Red-edge GCVI\n",
        "    df['RDGCVI1'] = df['B08'] / df['B05'] - 1\n",
        "    df['RDGCVI2'] = df['B08'] / df['B06'] - 1\n",
        "\n",
        "    # Meris Terrestrial Chlorophyll Index (MTCI)\n",
        "    df['MTCI'] = (df['B08'] - df['B05']) / (df['B05'] - df['B04'])\n",
        "    df['MTCI2'] = (df['B06'] - df['B05']) / (df['B05'] - df['B04'])\n",
        "\n",
        "    # Red-edge Inflection Point (REIP)\n",
        "    df['REIP'] = 700 + 40 * (((df['B04'] + df['B07']) / 2) - df['B05']) / (df['B07'] - df['B05'])\n",
        "\n",
        "    # Normalized Burn Ratio (NBR)\n",
        "    df['NBR1'] = (df['B08'] - df['B11']) / (df['B08'] + df['B11'])\n",
        "    df['NBR2'] = (df['B08'] - df['B12']) / (df['B08'] + df['B12'])\n",
        "\n",
        "    # Normalized Difference Tillage Index (NDTI)\n",
        "    df['NDTI'] = (df['B11'] - df['B12']) / (df['B11'] + df['B12'])\n",
        "\n",
        "    # Canopy Chlorophyll Content Index (CRC)\n",
        "    df['CRC'] = (df['B11'] - df['B03']) / (df['B11'] + df['B03'])\n",
        "\n",
        "    # Soil Tillage Index (STI)\n",
        "    df['STI'] = df['B11'] / df['B12']\n",
        "\n",
        "    return df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MNLOv16tufpS"
      },
      "source": [
        "We also need function for spatial and temporal aggregation to reduce the dimensionality of the dataset:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "TBXtxAafufpS"
      },
      "outputs": [],
      "source": [
        "#@title Define spatial `spatial_median_aggregation(df, bands)` and temporal `period_based_aggregation(df, bands)` aggregation function (Run cell)\n",
        "def spatial_median_aggregation(df, bands):\n",
        "    \"\"\"\n",
        "    Aggregate data by field and time, using the median of band values.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    df : pandas.DataFrame\n",
        "        Input DataFrame with 'field', 'time', and band columns.\n",
        "    bands : list\n",
        "        List of band columns to be aggregated.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    pandas.DataFrame\n",
        "        Aggregated DataFrame with median band values.\n",
        "    \"\"\"\n",
        "    # Calculate median of band values for each unique 'field' and 'time'\n",
        "    agg_df = df.groupby(['field', 'time'])[bands].median().reset_index()\n",
        "\n",
        "    # Drop duplicate entries for each unique 'field' and 'time', and remove band columns\n",
        "    unique_df = df.drop_duplicates(['field', 'time']).drop(bands, axis=1)\n",
        "\n",
        "    # Merge aggregated DataFrame with unique DataFrame\n",
        "    return pd.merge(agg_df, unique_df, on=['field', 'time'])\n",
        "\n",
        "\n",
        "def period_based_aggregation(df, bands):\n",
        "    \"\"\"\n",
        "    Aggregate data by field and defined time periods, using the mean of band values.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    df : pandas.DataFrame\n",
        "        Input DataFrame with 'field', 'time', and band columns.\n",
        "    bands : list\n",
        "        List of band columns to be aggregated.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    pandas.DataFrame\n",
        "        Aggregated DataFrame with mean band values for each time period.\n",
        "    \"\"\"\n",
        "    # Define time periods\n",
        "    periods = {\n",
        "        'p1': pd.to_datetime(['2019-06-06']),\n",
        "        'p2': pd.to_datetime(['2019-07-01', '2019-07-06', '2019-07-11', '2019-07-21']),\n",
        "        'p3': pd.to_datetime(['2019-08-05', '2019-08-15', '2019-08-25']),\n",
        "        'p4': pd.to_datetime(['2019-09-09', '2019-09-19', '2019-09-24', '2019-10-04']),\n",
        "        'p5': pd.to_datetime(['2019-11-03'])\n",
        "    }\n",
        "\n",
        "    # Assign period labels based on 'time'\n",
        "    for period, dates in periods.items():\n",
        "        df.loc[df['time'].isin(dates), 'period'] = period\n",
        "\n",
        "    # Calculate mean of band values for each unique 'field' and 'period'\n",
        "    period_agg_df = df.groupby(['field', 'period'])[bands].mean().reset_index()\n",
        "\n",
        "    # Drop duplicate entries for each unique 'field' and 'period', and remove band columns\n",
        "    unique_df = df.drop_duplicates(['field', 'period']).drop(bands, axis=1)\n",
        "\n",
        "    # Merge aggregated DataFrame with unique DataFrame\n",
        "    return pd.merge(period_agg_df, unique_df, on=['field', 'period'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hR78FnclufpS"
      },
      "source": [
        "Finally, we create functions to pivot the table (making periods into columns) and another function that runs the steps and splits the dataframe into `X` and `y`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "YSdcas3mufpS"
      },
      "outputs": [],
      "source": [
        "#@title Define helper functions\n",
        "def pivot_dataframe(df):\n",
        "    \"\"\"\n",
        "    Pivot the DataFrame so that each time period becomes a separate column.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    df : pandas.DataFrame\n",
        "        Input DataFrame with 'field', 'period', and other columns.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    pandas.DataFrame\n",
        "        Pivoted DataFrame with each 'period' as a separate column.\n",
        "    \"\"\"\n",
        "    return df.pivot(index=['field', 'crop', 'lat', 'lon'], columns='period').fillna(-1).reset_index()\n",
        "\n",
        "\n",
        "def process_dataframe(df, bands):\n",
        "    \"\"\"\n",
        "    Process the DataFrame by calculating indices, aggregating data, and pivoting.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    df : pandas.DataFrame\n",
        "        Input DataFrame with 'field', 'time', and band columns.\n",
        "    bands : list\n",
        "        List of band columns to be processed.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    X : pandas.DataFrame\n",
        "        Processed DataFrame with features for machine learning.\n",
        "    y : pandas.Series\n",
        "        Target labels for machine learning.\n",
        "    \"\"\"\n",
        "    # Calculate spectral indices\n",
        "    df = calculate_indices(df)\n",
        "\n",
        "    # Aggregate data by field and time using spatial median\n",
        "    df = spatial_median_aggregation(df, bands)\n",
        "\n",
        "    # Aggregate data by field and time period using mean\n",
        "    df = period_based_aggregation(df, bands)\n",
        "\n",
        "    # Calculate average latitude and longitude for each field\n",
        "    lat_lon_agg = df.groupby('field')[['lat', 'lon']].mean().reset_index()\n",
        "\n",
        "    # Merge aggregated DataFrame with latitude and longitude DataFrame\n",
        "    df = pd.merge(df.drop(columns=['lat', 'lon']), lat_lon_agg, on='field', how='left')\n",
        "\n",
        "    # Pivot DataFrame to have each period as a separate column\n",
        "    df = pivot_dataframe(df)\n",
        "\n",
        "    # Flatten multi-level column names\n",
        "    df.columns = [''.join(col).strip() if isinstance(col, tuple) else col for col in df.columns.values]\n",
        "\n",
        "    # Select columns to keep\n",
        "    columns_to_keep = ['field', 'lat', 'lon', 'crop'] + [col for col in df.columns if col.endswith(('p1', 'p2', 'p3', 'p4', 'p5')) and not col.startswith(('time', 'lat', 'lon', 'crop'))]\n",
        "    df = df[columns_to_keep]\n",
        "\n",
        "    # Split DataFrame into features (X) and target labels (y)\n",
        "    X, y = df.drop([\"crop\"], axis=1), df[\"crop\"]\n",
        "\n",
        "    return X, y"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hkQ_RIVqufpT"
      },
      "source": [
        "Let's prepare the training, validation, and test arrays:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xpGiMfIKufpT"
      },
      "outputs": [],
      "source": [
        "# Set the band columns'\n",
        "bands = ['B04', 'B05', 'B06', 'B07', 'B08', 'B11', 'B12', 'NDVI', 'RDNDVI1', 'RDNDVI2', 'GCVI', 'RDGCVI1', 'RDGCVI2', 'MTCI', 'MTCI2', 'REIP', 'NBR1', 'NBR2', 'NDTI', 'CRC', 'STI']\n",
        "\n",
        "# Prepare the dataset\n",
        "print(\"Processing `train` ...\")\n",
        "X_train, y_train = process_dataframe(train, bands)\n",
        "\n",
        "print(\"Processing `val` ...\")\n",
        "X_val, y_val = process_dataframe(val, bands)\n",
        "\n",
        "print(\"Processing `test` ...\")\n",
        "X_test, y_test = process_dataframe(test, bands)\n",
        "\n",
        "# Print the shapes\n",
        "X_train.shape, y_train.shape, X_val.shape, y_val.shape, X_test.shape, y_test.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NeoIeHOrufpT"
      },
      "source": [
        "Now, let's conduct random hyperparameter search with cross-validation using the `LightGBM` estimator:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CvrS6SMTufpT"
      },
      "outputs": [],
      "source": [
        "from sklearn.model_selection import RandomizedSearchCV\n",
        "from sklearn.metrics import make_scorer\n",
        "import lightgbm as lgb\n",
        "\n",
        "# Define the LightGBM model\n",
        "model = lgb.LGBMClassifier(objective=\"multiclass\", verbose=-1, num_class=7)\n",
        "banned_cols = [\"field\"]\n",
        "\n",
        "# Define the hyperparameters space\n",
        "param_dist = {\n",
        "    'num_leaves': [31, 127, 200, 300],\n",
        "    'reg_alpha': [0.1, 0.5],\n",
        "    'min_data_in_leaf': [30, 50, 100, 300, 400],\n",
        "    'lambda_l1': [0, 1, 1.5],\n",
        "    'lambda_l2': [0, 1]\n",
        "}\n",
        "\n",
        "# Define the scorer\n",
        "scorer = make_scorer(log_loss, greater_is_better=False, needs_proba=True)\n",
        "\n",
        "# Randomized Search for hyperparameter tuning\n",
        "random_search = RandomizedSearchCV(model, param_distributions=param_dist, n_iter=15, scoring=scorer, cv=3, verbose=1, n_jobs=-1)\n",
        "random_search.fit(X_train.drop(banned_cols, axis=1), y_train)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tiJVkLcVufpT"
      },
      "source": [
        "#### Evaluation\n",
        "\n",
        "We re-train the best estimator on the training data and get the validation cross-entropy:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "R5Jc_nprufpT"
      },
      "outputs": [],
      "source": [
        "# Create the LightGBM model instance with the best hyperparameters\n",
        "model = lgb.LGBMClassifier(objective=\"multiclass\", num_class=7, verbose=-1, **random_search.best_params_)\n",
        "\n",
        "# Fit the model to the training set\n",
        "model.fit(X_train.drop(banned_cols, axis=1), y_train)\n",
        "\n",
        "# Predict the validation set results\n",
        "y_val_hat = model.predict_proba(X_val.drop(banned_cols, axis=1))\n",
        "\n",
        "# Report cross-entropy\n",
        "print(f\"Cross-entropy with best hyperparameters is {log_loss(y_val, y_val_hat):.5f}\")\n",
        "print(f\"It is {100*(log_loss(y_val, y_val_hat) - baseline_ce)/baseline_ce:.2f}% better than the baseline\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2XvrrCrEwOcW"
      },
      "source": [
        "#### Investigating Class-Imbalances\n",
        "\n",
        "Let's report the following metrics on the combination of validation + test points:\n",
        "- `Precision`\n",
        "- `Recall`\n",
        "- `F1`\n",
        "- `Confusion matrix`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8F01HPUWwOcW"
      },
      "outputs": [],
      "source": [
        "# Predict the validation set results\n",
        "y_test_hat = model.predict(pd.concat([X_val, X_test]).drop(banned_cols, axis=1))\n",
        "y_test_arr = pd.concat([y_val, y_test]).values\n",
        "y_test_hat.shape, y_test_arr.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FI8SB5x_d4ST"
      },
      "source": [
        "<font color=\"red\">Exercise 5 Calculate precision, recall, and F1 score?</font> (*Hint: for average variable use \"average = 'weighted*, for more information you can check example: precision_score?)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CAMwMqZwe0Lv"
      },
      "outputs": [],
      "source": [
        "#@title Exercise 5 Calculate precision, recall, and F1 score?\n",
        "from sklearn.metrics import precision_score, recall_score, f1_score, confusion_matrix\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Calculate precision, recall, and F1 score\n",
        "precision = #update me\n",
        "recall = #update me\n",
        "f1 = #update me\n",
        "\n",
        "assert precision == 0.45818711816733754\n",
        "assert recall == 0.5503048780487805\n",
        "assert f1 ==  0.48797569337092017\n",
        "\n",
        "print(f\"Precision: {precision}\")\n",
        "print(f\"Recall: {recall}\")\n",
        "print(f\"F1 Score: {f1}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "zV7XpGLSwOcX"
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import precision_score, recall_score, f1_score, confusion_matrix\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "# @title Answer to Exercise 5 (Try not to run until you've given it a good try!')\n",
        "# Calculate precision, recall, and F1 score\n",
        "precision = precision_score(y_test_arr, y_test_hat, average='weighted')\n",
        "recall = recall_score(y_test_arr, y_test_hat, average='weighted')\n",
        "f1 = f1_score(y_test_arr, y_test_hat, average='weighted')\n",
        "\n",
        "assert precision == 0.45818711816733754\n",
        "assert recall == 0.5503048780487805\n",
        "assert f1 ==  0.48797569337092017\n",
        "\n",
        "print(f\"Precision: {precision}\")\n",
        "print(f\"Recall: {recall}\")\n",
        "print(f\"F1 Score: {f1}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TBcD5ZFXwOcX"
      },
      "outputs": [],
      "source": [
        "# Calculate confusion matrix\n",
        "cm = confusion_matrix(y_test_arr, y_test_hat, normalize=\"true\")\n",
        "\n",
        "# Plot confusion matrix\n",
        "plt.figure(figsize=(10, 7))\n",
        "sns.heatmap(cm, annot=True, xticklabels=id_to_name.values(), yticklabels=id_to_name.values())\n",
        "plt.xlabel('Predicted')\n",
        "plt.ylabel('True')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T2BSN1fuwOcX"
      },
      "source": [
        "Except for `maize` (which is majority class), we are not doing well classifiying the other minority crop classes. Since cross-entropy does not mitigate against class imbalance, we still get a good score.\n",
        "\n",
        "*Hint: try changing `LGBMClassifier`'s `class_weight` attribute to `balanced`.*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hKHoWqIkwOcX"
      },
      "source": [
        "#### XAI\n",
        "\n",
        "What are the most important periods and indices?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KJQC7gl6wOcX"
      },
      "outputs": [],
      "source": [
        "import shap\n",
        "\n",
        "# Prepare the validation + test data for the model\n",
        "X_vt = pd.concat([X_val, X_test])\n",
        "\n",
        "# explain the model's predictions using SHAP\n",
        "explainer = shap.TreeExplainer(model)\n",
        "shap_values = explainer.shap_values(X_vt.drop(banned_cols, axis=1))\n",
        "\n",
        "shap.summary_plot(shap_values, X_vt.drop(banned_cols, axis=1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jX_Fx_jjwOcX"
      },
      "source": [
        "Let's figure out which periods are the most important:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_Kq9mVzvwOcX"
      },
      "outputs": [],
      "source": [
        "# Compute the absolute SHAP values for each feature\n",
        "abs_shap_values = np.sum(np.abs(shap_values), axis=(0, 1))\n",
        "\n",
        "# Get the feature names\n",
        "feature_names = X_vt.drop(banned_cols, axis=1).columns\n",
        "\n",
        "# Create a DataFrame linking feature names to their importance\n",
        "feature_importances = pd.DataFrame({\n",
        "    'feature': feature_names,\n",
        "    'importance': abs_shap_values\n",
        "})\n",
        "\n",
        "# Sort the DataFrame by importance in descending order\n",
        "feature_importances = feature_importances.sort_values('importance', ascending=False)\n",
        "\n",
        "# Drop `lat` and `lon` from the dataset\n",
        "feature_importances = feature_importances[feature_importances[\"feature\"].isin([\"lat\", \"lon\"]) == False]\n",
        "\n",
        "# Normalize the feature importances to sum to one\n",
        "feature_importances['importance'] = feature_importances['importance'] / feature_importances['importance'].sum()\n",
        "\n",
        "# Split the feature name into `index` and `period` (period is the last two characters)\n",
        "feature_importances['period'] = feature_importances['feature'].str[-2:]\n",
        "feature_importances['index'] = feature_importances['feature'].str[:-2]\n",
        "feature_importances = feature_importances.drop(\"feature\", axis=1)\n",
        "\n",
        "# Get the most important periods separately by aggregating their importance\n",
        "periods = feature_importances.drop(\"index\", axis=1).groupby('period').sum().sort_values('importance', ascending=False)\n",
        "\n",
        "# Get the most important bands separately by aggregating their importance\n",
        "bands = feature_importances.drop(\"period\", axis=1).groupby('index').sum().sort_values('importance', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X53rh077wOcX"
      },
      "outputs": [],
      "source": [
        "periods"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "39eZAKSbwOcX"
      },
      "outputs": [],
      "source": [
        "bands.iloc[:10]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EWuTq2lUwOcX"
      },
      "source": [
        "We highlight the following top indices:\n",
        "\n",
        "- `NBR2`: Normalized Burn Ratio 2. It is used in remote sensing to identify burned areas. It uses the Near Infrared (NIR) and Shortwave Infrared (SWIR) portions of the electromagnetic spectrum.\n",
        "- `CRC`: Canopy Chlorophyll Content Index. It is used to estimate the chlorophyll content in plant canopies.\n",
        "- `B01`: Band 1 of Sentinel-2. This is a coastal aerosol band and captures light in the blue portion of the electromagnetic spectrum.\n",
        "- `B08`: Band 8 of Sentinel-2. This is a NIR (Near Infrared) band. It is often used in vegetation analysis as it reflects maximum light in healthy vegetation.\n",
        "- `GCVI`: Green Chlorophyll Vegetation Index. It is used to measure the chlorophyll content of vegetation by using the Green and Near-Infrared bands.\n",
        "- `REIP`: Red Edge Inflection Point. It is used in vegetation studies to indicate the boundary between the red and NIR region of the spectrum where vegetation has a strong reflection.\n",
        "- `B8A`: Band 8a of Sentinel-2. This is a Narrow NIR (Near Infrared) band. It is used to study water bodies and vegetation.\n",
        "- `B06`: Band 6 of Sentinel-2. It is a red-edge band, which is useful for vegetation studies.\n",
        "- `RDNDVI2`: Ratio Divergence Normalized Difference Vegetation Index 2. This is presumably a custom vegetation index that uses a ratio and divergence calculation to normalize the vegetation index.\n",
        "- `B09`: Band 9 of Sentinel-2. It is a water vapor band, useful for atmospheric studies."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4DJV_Hq0ufpT"
      },
      "source": [
        "#### Inference\n",
        "\n",
        "In this section, we will report the final metrics on the validation set and visualize the farms with their crop types:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KXh97bVmufpT"
      },
      "outputs": [],
      "source": [
        "# Predict on the test set\n",
        "y_test_pred = model.predict_proba(X_test.drop(banned_cols, axis=1))\n",
        "y_test_pred.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CI1h-n6hufpT"
      },
      "outputs": [],
      "source": [
        "# Export the results\n",
        "report = X_test[[\"field\"]].copy()\n",
        "\n",
        "# Create the Crop_ID_1,Crop_ID_2,Crop_ID_3,Crop_ID_4,Crop_ID_5,Crop_ID_6,Crop_ID_7 columns and assign the predictions\n",
        "cols = ['Crop_ID_1','Crop_ID_2','Crop_ID_3','Crop_ID_4','Crop_ID_5','Crop_ID_6','Crop_ID_7']\n",
        "report[cols] = y_test_pred\n",
        "report"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uZuh89XEufpT"
      },
      "source": [
        "Let's visualize the predicted test farms:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "t83GTLUBufpT"
      },
      "outputs": [],
      "source": [
        "from shapely.geometry import Point, LineString, Polygon\n",
        "\n",
        "def create_geometry(df):\n",
        "    coords = list(zip(df.lon, df.lat))\n",
        "    if len(coords) == 1: return Point(coords[0])\n",
        "    elif len(coords) == 2: return LineString(coords)\n",
        "    else: return Polygon(coords)\n",
        "\n",
        "# Create the polygons from the test set\n",
        "d = test.copy()\n",
        "cols = [\"field\", \"lat\", \"lon\"]\n",
        "d = d[cols].drop_duplicates()\n",
        "d = d.groupby('field').apply(create_geometry).reset_index().rename(columns={0: \"geometry\"})\n",
        "\n",
        "# Create the dataframe to hold the pixel locations and the predicted crop types\n",
        "report = X_test.copy()\n",
        "report = report[[\"field\"]]\n",
        "report[\"crop\"] = y_test_pred.argmax(axis=1) + 1\n",
        "\n",
        "# Merge the two dataframes\n",
        "report = report.merge(d, on=\"field\", how=\"left\").rename(columns={0: \"geometry\"})\n",
        "report = gpd.GeoDataFrame(report, geometry=\"geometry\")\n",
        "\n",
        "# Replace the 'crop' column with mapped names\n",
        "report['crop'] = report['crop'].map(id_to_name)\n",
        "\n",
        "# Plot\n",
        "fig, ax = plt.subplots(figsize=(10, 10))\n",
        "\n",
        "# Plot the GeoDataFrame using the 'crop' column to color the polygons\n",
        "report.plot(column=\"crop\", legend=True, ax=ax, cmap=\"Accent\")\n",
        "\n",
        "# Add a basemap\n",
        "ctx.add_basemap(ax, crs=\"EPSG:4326\", source=ctx.providers.Stamen.Terrain)\n",
        "\n",
        "# Show the plot\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A5w0WByCUpdb"
      },
      "source": [
        "---\n",
        "\n",
        "### 2.3 Model 2: Deep Learning <font color='orange'>`Intermediate`</font>\n",
        "\n",
        "<center><img src=\"https://drive.google.com/u/0/uc?id=1pd_-Azeunh7O7SoxXWhrZ5_f9O4wq1tU&export=download\" width=\"750px;\" /></center>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t60qAoCEbtm_"
      },
      "source": [
        "#### Data Preprocessing\n",
        "\n",
        "We have conducted the following steps to produce the `X.npy` and `y.npy` files:\n",
        "1. **Data normalization**: squared-root all bands in the Sentinel-2 imagery.\n",
        "2. **Farm masking**: Added an extra binary band that shows where the farm pixels are.\n",
        "3. **Farm Extraction**: For each field polygon, we get its center location and crop a 32x32 window from the original image to create the patches.\n",
        "4. **Patch Standardization**: Each patch was standardized using the temporal band-wise mean and standard deviation.\n",
        "5. **Creating a single Cube**: stacked all patches and their labels into 2 unified cubes and saved them as `X.npy` and `y.npy`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xebQiqiVErI2"
      },
      "source": [
        "#### Class Implementations\n",
        "\n",
        "As we now have the NumPy arrays for training, we want to implement a sequence-to-one classification model.\n",
        "\n",
        "Here are the choices that we are going to implement:\n",
        "- Encode each image using a pre-trained encoder (ResNet18).\n",
        "- Pass the sequence of encodings to a 3-layer Bi-directional GRU.\n",
        "- Take the final concatenated output representation from the GRU and pass it through a fully-connected layer to predict the final class probabilities (7 classes).\n",
        "- Use cross entropy as the loss function.\n",
        "- Conduct data augmentation to regularize the model.\n",
        "- Export the validation results"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "28IXH4N4Updb"
      },
      "source": [
        "#### Dataset\n",
        "\n",
        "Let's start with the `FieldSequenceDataset` class.\n",
        "\n",
        "We have a series of images taken of a field over time, and each of these sequences of images corresponds to a `label` (e.g., the type of crop in the field). The `FieldSequenceDataset` class is designed to manage and provide easy access to these sequences of images and their corresponding labels.\n",
        "\n",
        "##### **Key Methods**\n",
        "- `__len__`: This returns the number of items in the dataset.\n",
        "- `__getitem__`: Given an index, it provides the images sequence and its label.\n",
        "- `plot`: lets you visualize a sequence of field images.\n",
        "\n",
        "##### **Why is this useful?**\n",
        "By structuring the data in this way, it becomes much easier to:\n",
        "- Feed data into machine learning models.\n",
        "- Apply consistent modifications to sequences.\n",
        "- Visualize and understand the data you're working with."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5ra4rNyhUpdb"
      },
      "outputs": [],
      "source": [
        "#@title Implement the `__len__` and `__getitem__` methods\n",
        "class FieldSequenceDataset(Dataset):\n",
        "    \"\"\"\n",
        "    A dataset class for sequences of field images.\n",
        "\n",
        "    Attributes:\n",
        "    - X: Numpy array containing image sequences.\n",
        "    - y: Labels associated with each image sequence.\n",
        "    - classes: List of class names/labels.\n",
        "    - transforms: Optional data augmentation operations.\n",
        "\n",
        "    Methods:\n",
        "    - __len__ : Returns the length of the dataset.\n",
        "    - __getitem__ : Fetches a data sample for a given index.\n",
        "    - plot: Plots an image sequence from a given sample.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        X,\n",
        "        y,\n",
        "        field_ids: List[int],\n",
        "        transforms: Optional[Callable] = None\n",
        "    ) -> None:\n",
        "        \"\"\"\n",
        "        Initializes the dataset object.\n",
        "\n",
        "        Parameters:\n",
        "        - X: Numpy array containing image sequences of shape (num_samples, num_images, height, width, bands).\n",
        "        - y: Numpy array containing labels for each sequence.\n",
        "        - field_ids: List of indices to subset the dataset. Defaults to None (use all data).\n",
        "        - transforms: Optional data augmentation operations.\n",
        "        \"\"\"\n",
        "\n",
        "        # Define class labels\n",
        "        self.classes = [str(i) for i in range(1, 8)]\n",
        "\n",
        "        # Instead of slicing the data, store the indices\n",
        "        self.field_ids = field_ids\n",
        "        self.X = X\n",
        "        self.y = y\n",
        "\n",
        "        # Set the data augmentation transforms\n",
        "        self.transforms = transforms\n",
        "\n",
        "    def __len__(self) -> int:\n",
        "        \"\"\"Returns the number of samples in the dataset.\"\"\"\n",
        "        pass\n",
        "\n",
        "    def __getitem__(self, index: int) -> dict[str, torch.Tensor]:\n",
        "        \"\"\"\n",
        "        Returns a data sample given an index.\n",
        "\n",
        "        Parameters:\n",
        "        - index: Index of the sample to fetch.\n",
        "\n",
        "        Returns:\n",
        "        Dictionary containing the image sequence and its associated label.\n",
        "        \"\"\"\n",
        "        pass\n",
        "\n",
        "    def plot(\n",
        "        self,\n",
        "        sample: Dict[str, Any],\n",
        "        show_titles: bool = True,\n",
        "        suptitle: Optional[str] = None,\n",
        "    ) -> plt.Figure:\n",
        "        \"\"\"\n",
        "        Plots an image sequence from a sample.\n",
        "\n",
        "        Parameters:\n",
        "        - sample: Dictionary containing an image sequence and its label.\n",
        "        - show_titles: Whether to display titles on the plots.\n",
        "        - suptitle: Optional overarching title for the entire plot.\n",
        "\n",
        "        Returns:\n",
        "        Matplotlib figure object.\n",
        "        \"\"\"\n",
        "\n",
        "        # Extract and normalize image sequence\n",
        "        sequence = sample['image'].numpy()[:, [3, 2, 1], :, :]\n",
        "        label = sample['label'].item()\n",
        "        min_vals = sequence.min(axis=(0, 2, 3), keepdims=True)\n",
        "        max_vals = sequence.max(axis=(0, 2, 3), keepdims=True)\n",
        "        sequence = (sequence - min_vals) / (max_vals - min_vals)\n",
        "\n",
        "        # Calculate layout for plotting multiple images\n",
        "        num_images = sequence.shape[0]\n",
        "        num_rows = int(np.ceil(num_images / 4.0))\n",
        "\n",
        "        # Create a figure and plot each image in the sequence\n",
        "        fig, axarr = plt.subplots(num_rows, 4, figsize=(15, 4 * num_rows))\n",
        "        if num_rows == 1:\n",
        "            axarr = np.expand_dims(axarr, axis=0)\n",
        "        for i in range(num_rows):\n",
        "            for j in range(4):\n",
        "                idx = i * 4 + j\n",
        "                if idx < num_images:\n",
        "                    ax = axarr[i, j]\n",
        "                    ax.imshow(sequence[idx].transpose(1, 2, 0))\n",
        "                    ax.axis('off')\n",
        "                    if show_titles and idx == num_images - 1:\n",
        "                        ax.set_title(f'Label: {self.classes[label]}')\n",
        "                else:\n",
        "                    axarr[i, j].axis('off')\n",
        "\n",
        "        # Set the optional overarching title\n",
        "        if suptitle:\n",
        "            fig.suptitle(suptitle, fontsize=16)\n",
        "\n",
        "        return fig"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "88vcKR1n6sqO",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title Answer to code task (Try not to peek until you've given it a good try!')\n",
        "class FieldSequenceDataset(Dataset):\n",
        "    \"\"\"\n",
        "    A dataset class for sequences of field images.\n",
        "\n",
        "    Attributes:\n",
        "    - X: Numpy array containing image sequences.\n",
        "    - y: Labels associated with each image sequence.\n",
        "    - classes: List of class names/labels.\n",
        "    - transforms: Optional data augmentation operations.\n",
        "\n",
        "    Methods:\n",
        "    - __len__ : Returns the length of the dataset.\n",
        "    - __getitem__ : Fetches a data sample for a given index.\n",
        "    - plot: Plots an image sequence from a given sample.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        X,\n",
        "        y,\n",
        "        field_ids: List[int],\n",
        "        transforms: Optional[Callable] = None\n",
        "    ) -> None:\n",
        "        \"\"\"\n",
        "        Initializes the dataset object.\n",
        "\n",
        "        Parameters:\n",
        "        - X: Numpy array containing image sequences of shape (num_samples, num_images, height, width, bands).\n",
        "        - y: Numpy array containing labels for each sequence.\n",
        "        - field_ids: List of indices to subset the dataset. Defaults to None (use all data).\n",
        "        - transforms: Optional data augmentation operations.\n",
        "        \"\"\"\n",
        "\n",
        "        # Define class labels\n",
        "        self.classes = [str(i) for i in range(1, 8)]\n",
        "\n",
        "        # Instead of slicing the data, store the indices\n",
        "        self.field_ids = field_ids\n",
        "        self.X = X\n",
        "        self.y = y\n",
        "\n",
        "        # Set the data augmentation transforms\n",
        "        self.transforms = transforms\n",
        "\n",
        "    def __len__(self) -> int:\n",
        "        \"\"\"Returns the number of samples in the dataset.\"\"\"\n",
        "        return len(self.field_ids)\n",
        "\n",
        "    def __getitem__(self, index: int) -> dict[str, torch.Tensor]:\n",
        "        \"\"\"\n",
        "        Returns a data sample given an index.\n",
        "\n",
        "        Parameters:\n",
        "        - index: Index of the sample to fetch.\n",
        "\n",
        "        Returns:\n",
        "        Dictionary containing the image sequence and its associated label.\n",
        "        \"\"\"\n",
        "\n",
        "        # Use the field_ids to fetch the relevant data\n",
        "        sequence = self.X[self.field_ids[index]]\n",
        "        label = self.y[self.field_ids[index]]\n",
        "\n",
        "        # Convert them to PyTorch tensors\n",
        "        sample = {'image': torch.tensor(sequence, dtype=torch.float32), 'label': torch.tensor(label, dtype=torch.long)}\n",
        "\n",
        "        return sample\n",
        "\n",
        "    def plot(\n",
        "        self,\n",
        "        sample: Dict[str, Any],\n",
        "        show_titles: bool = True,\n",
        "        suptitle: Optional[str] = None,\n",
        "    ) -> plt.Figure:\n",
        "        \"\"\"\n",
        "        Plots an image sequence from a sample.\n",
        "\n",
        "        Parameters:\n",
        "        - sample: Dictionary containing an image sequence and its label.\n",
        "        - show_titles: Whether to display titles on the plots.\n",
        "        - suptitle: Optional overarching title for the entire plot.\n",
        "\n",
        "        Returns:\n",
        "        Matplotlib figure object.\n",
        "        \"\"\"\n",
        "\n",
        "        # Extract and normalize image sequence\n",
        "        sequence = sample['image'].numpy()[:, [3, 2, 1], :, :]\n",
        "        label = sample['label'].item()\n",
        "        min_vals = sequence.min(axis=(0, 2, 3), keepdims=True)\n",
        "        max_vals = sequence.max(axis=(0, 2, 3), keepdims=True)\n",
        "        sequence = (sequence - min_vals) / (max_vals - min_vals)\n",
        "\n",
        "        # Calculate layout for plotting multiple images\n",
        "        num_images = sequence.shape[0]\n",
        "        num_rows = int(np.ceil(num_images / 4.0))\n",
        "\n",
        "        # Create a figure and plot each image in the sequence\n",
        "        fig, axarr = plt.subplots(num_rows, 4, figsize=(15, 4 * num_rows))\n",
        "        if num_rows == 1:\n",
        "            axarr = np.expand_dims(axarr, axis=0)\n",
        "        for i in range(num_rows):\n",
        "            for j in range(4):\n",
        "                idx = i * 4 + j\n",
        "                if idx < num_images:\n",
        "                    ax = axarr[i, j]\n",
        "                    ax.imshow(sequence[idx].transpose(1, 2, 0))\n",
        "                    ax.axis('off')\n",
        "                    if show_titles and idx == num_images - 1:\n",
        "                        ax.set_title(f'Label: {self.classes[label]}')\n",
        "                else:\n",
        "                    axarr[i, j].axis('off')\n",
        "\n",
        "        # Set the optional overarching title\n",
        "        if suptitle:\n",
        "            fig.suptitle(suptitle, fontsize=16)\n",
        "\n",
        "        return fig"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qs81pjALUpdb"
      },
      "source": [
        "#### Data Module\n",
        "\n",
        "In machine learning, it's common to split our data into three parts:\n",
        "- **Training set**: to train our model.\n",
        "- **Validation set**: to tune and optimize our model's parameters.\n",
        "- **Test set**: to check how well our model will perform in real-world scenarios.\n",
        "\n",
        "The next class we are going to implement (`FieldDataModule`) is a utility to handle and organize this process, especially for our field image sequences.\n",
        "\n",
        "##### **Why is this helpful?**\n",
        "- Organization: Keeps data management clean and structured.\n",
        "- Flexibility: Easily change the size of each dataset or the batch size.\n",
        "- Efficiency: Facilitates parallel data loading, so you spend less time waiting and more time training."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hC0_IDJKUpdb"
      },
      "outputs": [],
      "source": [
        "#@title Implement the `setup()`, `train_dataloader()`, `val_dataloader()`, and `test_dataloader()` methods\n",
        "class FieldDataModule(LightningDataModule):\n",
        "    \"\"\"\n",
        "    PyTorch Lightning data module for handling field sequence data.\n",
        "\n",
        "    This class helps in loading and splitting the dataset into train, validation, and test sets.\n",
        "\n",
        "    Attributes:\n",
        "    - root: The path to the root directory containing the data.\n",
        "    - batch_size: Size of the batches during training.\n",
        "    - workers: Number of workers for data loading.\n",
        "    - X: Numpy array containing image sequences.\n",
        "    - y: Numpy array containing labels for each sequence.\n",
        "    - train_ids, val_ids, test_ids: Lists containing indices for the train, validation, and test splits respectively.\n",
        "    - train_ds, val_ds, test_ds: Dataset objects for the train, validation, and test sets.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        root: str,\n",
        "        train_size: float = 0.8,\n",
        "        val_size: float = 0.1,\n",
        "        test_size: float = 0.1,\n",
        "        batch_size: int = 8,\n",
        "        workers: int = 4,\n",
        "    ):\n",
        "        super().__init__()\n",
        "\n",
        "        # Define directory path and loading configurations\n",
        "        self.root = Path(root)\n",
        "        self.batch_size = batch_size\n",
        "        self.workers = workers\n",
        "\n",
        "        # Load the dataset into memory\n",
        "        self.X = np.load(self.root / \"X.npy\")\n",
        "        self.y = np.load(self.root / \"y.npy\")\n",
        "\n",
        "        # Randomly shuffle field IDs for dataset split\n",
        "        all_field_ids = list(range(3280))\n",
        "        shuffle(all_field_ids)\n",
        "\n",
        "        # Split the dataset into train, validation, and test sets based on provided ratios\n",
        "        self.train_ids, temp_ids = train_test_split(all_field_ids, test_size=1 - train_size, random_state=42)\n",
        "        self.val_ids, self.test_ids = train_test_split(temp_ids, test_size=test_size / (test_size + val_size), random_state=42)\n",
        "\n",
        "        # Setup datasets\n",
        "        self.setup()\n",
        "\n",
        "    def setup(self, stage=None):\n",
        "        \"\"\"\n",
        "        Prepare datasets for training, validation, and testing.\n",
        "\n",
        "        Uses the field IDs generated during initialization to subset the full dataset.\n",
        "        \"\"\"\n",
        "        pass\n",
        "\n",
        "    def train_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the training dataset.\"\"\"\n",
        "        pass\n",
        "\n",
        "    def val_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the validation dataset.\"\"\"\n",
        "        pass\n",
        "\n",
        "    def test_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the test dataset.\"\"\"\n",
        "        pass"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "T0xfWX6_6sqO",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title Answer to code task (Try not to peek until you've given it a good try!')\n",
        "class FieldDataModule(LightningDataModule):\n",
        "    \"\"\"\n",
        "    PyTorch Lightning data module for handling field sequence data.\n",
        "\n",
        "    This class helps in loading and splitting the dataset into train, validation, and test sets.\n",
        "\n",
        "    Attributes:\n",
        "    - root: The path to the root directory containing the data.\n",
        "    - batch_size: Size of the batches during training.\n",
        "    - workers: Number of workers for data loading.\n",
        "    - X: Numpy array containing image sequences.\n",
        "    - y: Numpy array containing labels for each sequence.\n",
        "    - train_ids, val_ids, test_ids: Lists containing indices for the train, validation, and test splits respectively.\n",
        "    - train_ds, val_ds, test_ds: Dataset objects for the train, validation, and test sets.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(\n",
        "        self,\n",
        "        root: str,\n",
        "        train_size: float = 0.8,\n",
        "        val_size: float = 0.1,\n",
        "        test_size: float = 0.1,\n",
        "        batch_size: int = 8,\n",
        "        workers: int = 4,\n",
        "    ):\n",
        "        super().__init__()\n",
        "\n",
        "        # Define directory path and loading configurations\n",
        "        self.root = Path(root)\n",
        "        self.batch_size = batch_size\n",
        "        self.workers = workers\n",
        "\n",
        "        # Load the dataset into memory\n",
        "        self.X = np.load(self.root / \"X.npy\")\n",
        "        self.y = np.load(self.root / \"y.npy\")\n",
        "\n",
        "        # Randomly shuffle field IDs for dataset split\n",
        "        all_field_ids = list(range(3280))\n",
        "        shuffle(all_field_ids)\n",
        "\n",
        "        # Split the dataset into train, validation, and test sets based on provided ratios\n",
        "        self.train_ids, temp_ids = train_test_split(all_field_ids, test_size=1 - train_size, random_state=42)\n",
        "        self.val_ids, self.test_ids = train_test_split(temp_ids, test_size=test_size / (test_size + val_size), random_state=42)\n",
        "\n",
        "        # Setup datasets\n",
        "        self.setup()\n",
        "\n",
        "    def setup(self, stage=None):\n",
        "        \"\"\"\n",
        "        Prepare datasets for training, validation, and testing.\n",
        "\n",
        "        Uses the field IDs generated during initialization to subset the full dataset.\n",
        "        \"\"\"\n",
        "        self.train_ds = FieldSequenceDataset(self.X, self.y, self.train_ids)\n",
        "        self.val_ds = FieldSequenceDataset(self.X, self.y, self.val_ids)\n",
        "        self.test_ds = FieldSequenceDataset(self.X, self.y, self.test_ids)\n",
        "\n",
        "    def train_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the training dataset.\"\"\"\n",
        "        return DataLoader(self.train_ds, batch_size=self.batch_size, num_workers=self.workers, shuffle=True)\n",
        "\n",
        "    def val_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the validation dataset.\"\"\"\n",
        "        return DataLoader(self.val_ds, batch_size=self.batch_size, num_workers=self.workers, shuffle=False)\n",
        "\n",
        "    def test_dataloader(self):\n",
        "        \"\"\"Returns a DataLoader object for the test dataset.\"\"\"\n",
        "        return DataLoader(self.test_ds, batch_size=self.batch_size, num_workers=self.workers, shuffle=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HXsq_Z_GUpdb"
      },
      "source": [
        "#### Model Architecture\n",
        "\n",
        "Now, to the model architecture. Imagine you have sequences of images captured over time, and you want to classify each sequence into a category. In our case, the sequence are images of a field, and we want to identify the crop grown in that field based on the sequence.\n",
        "\n",
        "##### **Main Components**\n",
        "- `ResNet18 Encoder`: Every image in the sequence is first processed by this encoder, which extracts important features from each image. Think of this as converting a detailed image into a summarized version that retains all the essential information.\n",
        "- `Bidirectional GRU`: Once we have the features for each image in the sequence, this component helps the model understand the order and relationship between these images. It looks at the sequence forwards and backwards, ensuring it captures patterns that emerge over time.\n",
        "- `Fully Connected Layer`: After understanding the sequence, this part of the model makes the final decision. It takes the output of the GRU and classifies the entire sequence into one of the categories."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eINHRkwMUpdb"
      },
      "outputs": [],
      "source": [
        "#@title Implement the `forward()` method\n",
        "class SequenceClassificationModel(nn.Module):\n",
        "    \"\"\"\n",
        "    Neural network model for sequence classification tasks.\n",
        "\n",
        "    This model consists of a ResNet18 encoder, a bidirectional GRU, and a fully connected classifier.\n",
        "    Given an input sequence of images, it outputs class probabilities for each sequence.\n",
        "\n",
        "    Attributes:\n",
        "    - encoder: ResNet18 encoder for feature extraction from each image in the sequence.\n",
        "    - gru: Bidirectional GRU to model temporal dependencies in the sequence of features.\n",
        "    - fc: Fully connected layer to produce class probabilities.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self, in_channels: int, input_size: int, hidden_size: int, num_layers: int, num_classes: int):\n",
        "        \"\"\"\n",
        "        Initialize the model.\n",
        "\n",
        "        Parameters:\n",
        "        - in_channels (int): Number of channels in the input images.\n",
        "        - input_size (int): Size of the encoded image features, which serves as the input size for the GRU.\n",
        "        - hidden_size (int): Number of units in the hidden layer of the GRU.\n",
        "        - num_layers (int): Number of recurrent layers in the GRU.\n",
        "        - num_classes (int): Number of output classes for classification.\n",
        "        \"\"\"\n",
        "        super(SequenceClassificationModel, self).__init__()\n",
        "\n",
        "        # Create ResNet18 encoder using timm, configured for the specified number of input channels.\n",
        "        self.encoder = timm.create_model(\n",
        "            \"resnet18\",\n",
        "            num_classes=0,  # Setting to 0 removes the classification head.\n",
        "            in_chans=in_channels,\n",
        "            pretrained=True\n",
        "        )\n",
        "\n",
        "        # Bidirectional GRU for modeling sequences.\n",
        "        self.gru = nn.GRU(input_size, hidden_size, num_layers, bidirectional=True, batch_first=True)\n",
        "\n",
        "        # Fully connected layer for outputting class probabilities.\n",
        "        self.fc = nn.Linear(2 * hidden_size, num_classes)  # 2 * hidden_size because the GRU is bidirectional.\n",
        "\n",
        "    def forward(self, x: torch.Tensor) -> torch.Tensor:\n",
        "        \"\"\"\n",
        "        Forward pass of the model.\n",
        "\n",
        "        Parameters:\n",
        "        - x (torch.Tensor): Input tensor of shape (batch_size, sequence_length, channels, height, width).\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Output tensor of class probabilities with shape (batch_size, num_classes).\n",
        "        \"\"\"\n",
        "        pass"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0VhOgvtR6sqO",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title Answer to code task (Try not to peek until you've given it a good try!')\n",
        "class SequenceClassificationModel(nn.Module):\n",
        "    \"\"\"\n",
        "    Neural network model for sequence classification tasks.\n",
        "\n",
        "    This model consists of a ResNet18 encoder, a bidirectional GRU, and a fully connected classifier.\n",
        "    Given an input sequence of images, it outputs class probabilities for each sequence.\n",
        "\n",
        "    Attributes:\n",
        "    - encoder: ResNet18 encoder for feature extraction from each image in the sequence.\n",
        "    - gru: Bidirectional GRU to model temporal dependencies in the sequence of features.\n",
        "    - fc: Fully connected layer to produce class probabilities.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self, in_channels: int, input_size: int, hidden_size: int, num_layers: int, num_classes: int):\n",
        "        \"\"\"\n",
        "        Initialize the model.\n",
        "\n",
        "        Parameters:\n",
        "        - in_channels (int): Number of channels in the input images.\n",
        "        - input_size (int): Size of the encoded image features, which serves as the input size for the GRU.\n",
        "        - hidden_size (int): Number of units in the hidden layer of the GRU.\n",
        "        - num_layers (int): Number of recurrent layers in the GRU.\n",
        "        - num_classes (int): Number of output classes for classification.\n",
        "        \"\"\"\n",
        "        super(SequenceClassificationModel, self).__init__()\n",
        "\n",
        "        # Create ResNet18 encoder using timm, configured for the specified number of input channels.\n",
        "        self.encoder = timm.create_model(\n",
        "            \"resnet18\",\n",
        "            num_classes=0,  # Setting to 0 removes the classification head.\n",
        "            in_chans=in_channels,\n",
        "            pretrained=True\n",
        "        )\n",
        "\n",
        "        # Bidirectional GRU for modeling sequences.\n",
        "        self.gru = nn.GRU(input_size, hidden_size, num_layers, bidirectional=True, batch_first=True)\n",
        "\n",
        "        # Fully connected layer for outputting class probabilities.\n",
        "        self.fc = nn.Linear(2 * hidden_size, num_classes)  # 2 * hidden_size because the GRU is bidirectional.\n",
        "\n",
        "    def forward(self, x: torch.Tensor) -> torch.Tensor:\n",
        "        \"\"\"\n",
        "        Forward pass of the model.\n",
        "\n",
        "        Parameters:\n",
        "        - x (torch.Tensor): Input tensor of shape (batch_size, sequence_length, channels, height, width).\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Output tensor of class probabilities with shape (batch_size, num_classes).\n",
        "        \"\"\"\n",
        "        batch_size, seq_len, C, H, W = x.shape\n",
        "\n",
        "        # Reshape to pass individual images through the encoder.\n",
        "        x = x.view(batch_size * seq_len, C, H, W)\n",
        "\n",
        "        # Pass through encoder to obtain encoded features.\n",
        "        x = self.encoder(x)\n",
        "\n",
        "        # Reshape back to sequence format for GRU.\n",
        "        x = x.view(batch_size, seq_len, -1)\n",
        "\n",
        "        # Pass through GRU. We only need the hidden state of the last layer.\n",
        "        _, h_n = self.gru(x)\n",
        "\n",
        "        # Reshape hidden states and concatenate the hidden states of the last layer's forward and backward passes.\n",
        "        h_n = h_n.view(self.gru.num_layers, 2, batch_size, self.gru.hidden_size)\n",
        "        final_hidden = torch.cat((h_n[-1, 0], h_n[-1, 1]), dim=-1)\n",
        "\n",
        "        # Pass through the classifier to get class probabilities.\n",
        "        output = self.fc(final_hidden)\n",
        "        return output"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GuAZRsAgUpdb"
      },
      "source": [
        "#### Data Augmentation\n",
        "\n",
        "We design the `SequenceAugmentationPipeline` class to address the challenge of consistent augmentation across sequences. We also want to introduce variability (augmentation) in our data for better model training.\n",
        "\n",
        "##### **Key Components**\n",
        "- Random Horizontal Flip: This might mirror your image as if you're looking at its reflection in a pond.\n",
        "- Random Vertical Flip: Imagine flipping your photo upside-down.\n",
        "- Random Rotation: This slightly rotates your image, just as if you'd tilted your camera a bit."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FBBlhXwCUpdc"
      },
      "outputs": [],
      "source": [
        "#@title Implement the `forward()` method\n",
        "class SequenceAugmentationPipeline(nn.Module):\n",
        "    \"\"\"\n",
        "    A data augmentation pipeline for sequences of images.\n",
        "\n",
        "    This module defines a set of transformations that are applied consistently across\n",
        "    all images in a sequence. This ensures that the spatial relationship between\n",
        "    images in a sequence remains consistent after augmentation.\n",
        "\n",
        "    Attributes:\n",
        "    - hflip: Random horizontal flip transformation.\n",
        "    - vflip: Random vertical flip transformation.\n",
        "    - rotate: Random rotation transformation.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self) -> None:\n",
        "        \"\"\"\n",
        "        Initialize the augmentation pipeline with desired transformations.\n",
        "        \"\"\"\n",
        "        super(SequenceAugmentationPipeline, self).__init__()\n",
        "\n",
        "        self.hflip = K.RandomHorizontalFlip()\n",
        "        self.vflip = K.RandomVerticalFlip()\n",
        "        self.rotate = K.RandomRotation(degrees=30)\n",
        "\n",
        "    def forward(self, input: torch.Tensor) -> torch.Tensor:\n",
        "        \"\"\"\n",
        "        Apply the transformations consistently across each image in the sequence.\n",
        "\n",
        "        Parameters:\n",
        "        - input (torch.Tensor): Input tensor of shape (batch_size, sequence_length, bands, height, width).\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Augmented tensor with the same shape as input.\n",
        "        \"\"\"\n",
        "        pass"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iOSfrfF76sqP",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title Answer to code task (Try not to peek until you've given it a good try!')\n",
        "class SequenceAugmentationPipeline(nn.Module):\n",
        "    \"\"\"\n",
        "    A data augmentation pipeline for sequences of images.\n",
        "\n",
        "    This module defines a set of transformations that are applied consistently across\n",
        "    all images in a sequence. This ensures that the spatial relationship between\n",
        "    images in a sequence remains consistent after augmentation.\n",
        "\n",
        "    Attributes:\n",
        "    - hflip: Random horizontal flip transformation.\n",
        "    - vflip: Random vertical flip transformation.\n",
        "    - rotate: Random rotation transformation.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self) -> None:\n",
        "        \"\"\"\n",
        "        Initialize the augmentation pipeline with desired transformations.\n",
        "        \"\"\"\n",
        "        super(SequenceAugmentationPipeline, self).__init__()\n",
        "\n",
        "        self.hflip = K.RandomHorizontalFlip()\n",
        "        self.vflip = K.RandomVerticalFlip()\n",
        "        self.rotate = K.RandomRotation(degrees=30)\n",
        "\n",
        "    def forward(self, input: torch.Tensor) -> torch.Tensor:\n",
        "        \"\"\"\n",
        "        Apply the transformations consistently across each image in the sequence.\n",
        "\n",
        "        Parameters:\n",
        "        - input (torch.Tensor): Input tensor of shape (batch_size, sequence_length, bands, height, width).\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Augmented tensor with the same shape as input.\n",
        "        \"\"\"\n",
        "        # Extract the shape parameters for the transformations from the first image\n",
        "        # in the sequence. This ensures consistent augmentation across all images.\n",
        "        hflip_params = self.hflip.forward_parameters(input[:, 0, ...].shape)\n",
        "        vflip_params = self.vflip.forward_parameters(input[:, 0, ...].shape)\n",
        "        rotate_params = self.rotate.forward_parameters(input[:, 0, ...].shape)\n",
        "\n",
        "        # Apply the transformations to each image in the sequence.\n",
        "        transformed_seq = []\n",
        "        for image in input.unbind(dim=1):\n",
        "            image = self.hflip(image, hflip_params)\n",
        "            image = self.vflip(image, vflip_params)\n",
        "            image = self.rotate(image, rotate_params)\n",
        "            transformed_seq.append(image)\n",
        "\n",
        "        # Combine the transformed images back into the sequence format.\n",
        "        output = torch.stack(transformed_seq, dim=1)\n",
        "\n",
        "        return output"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oxsn6z15Amf7"
      },
      "source": [
        "#### Task\n",
        "\n",
        "The following class helps us manage the training process using Pytorch-Lightning:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "btgeI9viUpdc"
      },
      "outputs": [],
      "source": [
        "#@title Create the task class\n",
        "class SequenceClassificationTask(LightningModule):\n",
        "    \"\"\"\n",
        "    Lightning module for the sequence classification task.\n",
        "\n",
        "    This module wraps the SequenceClassificationModel for training, validation, and testing.\n",
        "    It also handles data augmentation using the SequenceAugmentationPipeline.\n",
        "\n",
        "    Attributes:\n",
        "    - model: The sequence classification model.\n",
        "    - loss_fn: Loss function for classification.\n",
        "    - learning_rate: Learning rate for the optimizer.\n",
        "    - aug: Data augmentation pipeline for training sequences.\n",
        "    \"\"\"\n",
        "\n",
        "    def __init__(self, input_size, hidden_size, in_channels=14, num_layers=3, num_classes=7, learning_rate=0.001):\n",
        "        \"\"\"\n",
        "        Initialize the lightning module.\n",
        "\n",
        "        Parameters:\n",
        "        - input_size (int): Size of the input to the GRU.\n",
        "        - hidden_size (int): Size of the GRU hidden state.\n",
        "        - in_channels (int, optional): Number of input channels to the model. Defaults to 14.\n",
        "        - num_layers (int, optional): Number of GRU layers. Defaults to 3.\n",
        "        - num_classes (int, optional): Number of classification classes. Defaults to 7.\n",
        "        - learning_rate (float, optional): Learning rate for the optimizer. Defaults to 0.001.\n",
        "        \"\"\"\n",
        "        super(SequenceClassificationTask, self).__init__()\n",
        "\n",
        "        self.model = SequenceClassificationModel(in_channels, input_size, hidden_size, num_layers, num_classes)\n",
        "        self.loss_fn = nn.CrossEntropyLoss()\n",
        "        self.learning_rate = learning_rate\n",
        "\n",
        "        # Define the data augmentation pipeline for training.\n",
        "        self.aug = SequenceAugmentationPipeline()\n",
        "\n",
        "    def forward(self, x):\n",
        "        \"\"\"\n",
        "        Forward pass through the model.\n",
        "\n",
        "        Parameters:\n",
        "        - x (torch.Tensor): Input tensor.\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Model predictions.\n",
        "        \"\"\"\n",
        "        return self.model(x)\n",
        "\n",
        "    def training_step(self, batch, batch_idx):\n",
        "        \"\"\"\n",
        "        Defines a single step during training.\n",
        "\n",
        "        Parameters:\n",
        "        - batch (dict): Batch of data.\n",
        "        - batch_idx (int): Index of the batch.\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Training loss.\n",
        "        \"\"\"\n",
        "        x, y = batch[\"image\"], batch[\"label\"]\n",
        "\n",
        "        # Apply data augmentation to the training data.\n",
        "        x = self.aug(x)\n",
        "\n",
        "        logits = self(x)\n",
        "        loss = self.loss_fn(logits, y)\n",
        "\n",
        "        # Log training loss to TensorBoard.\n",
        "        self.log(\"train_loss\", loss)\n",
        "\n",
        "        return loss\n",
        "\n",
        "    def validation_step(self, batch, batch_idx):\n",
        "        \"\"\"\n",
        "        Defines a single step during validation.\n",
        "\n",
        "        Parameters:\n",
        "        - batch (dict): Batch of data.\n",
        "        - batch_idx (int): Index of the batch.\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Validation loss.\n",
        "        \"\"\"\n",
        "        x, y = batch[\"image\"], batch[\"label\"]\n",
        "        logits = self(x)\n",
        "        loss = self.loss_fn(logits, y)\n",
        "\n",
        "        # Log validation loss to TensorBoard.\n",
        "        self.log(\"val_loss\", loss)\n",
        "\n",
        "        return loss\n",
        "\n",
        "    def test_step(self, batch, batch_idx):\n",
        "        \"\"\"\n",
        "        Defines a single step during testing.\n",
        "\n",
        "        Parameters:\n",
        "        - batch (dict): Batch of data.\n",
        "        - batch_idx (int): Index of the batch.\n",
        "\n",
        "        Returns:\n",
        "        - torch.Tensor: Testing loss.\n",
        "        \"\"\"\n",
        "        x, y = batch[\"image\"], batch[\"label\"]\n",
        "        logits = self(x)\n",
        "        loss = self.loss_fn(logits, y)\n",
        "\n",
        "        # Log testing loss to TensorBoard.\n",
        "        self.log(\"test_loss\", loss)\n",
        "\n",
        "        return loss\n",
        "\n",
        "    def configure_optimizers(self):\n",
        "        \"\"\"\n",
        "        Configures the optimizer(s) and learning rate scheduler(s).\n",
        "\n",
        "        Returns:\n",
        "        - Dict: Contains optimizer and learning rate scheduler information.\n",
        "        \"\"\"\n",
        "        optimizer = torch.optim.AdamW(self.model.parameters(), lr=self.learning_rate)\n",
        "\n",
        "        # Define a learning rate scheduler that reduces the learning rate when the validation loss plateaus.\n",
        "        scheduler = ReduceLROnPlateau(optimizer, patience=5)\n",
        "\n",
        "        return {\n",
        "            \"optimizer\": optimizer,\n",
        "            \"lr_scheduler\": {\n",
        "                \"scheduler\": scheduler,\n",
        "                \"monitor\": \"val_loss\",\n",
        "            },\n",
        "        }"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1NxkGKgTUpdc"
      },
      "source": [
        "#### **Training**\n",
        "\n",
        "We train the model:"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Set training config\n",
        "# root = \"./files\"\n",
        "# experiment_name = \"seq2one-poc\"\n",
        "# gpu = 0\n",
        "# min_epochs, max_epochs = 3, 30\n",
        "\n",
        "# # Set the hyperparameters\n",
        "# batch_size = 16\n",
        "# hidden_size = 256\n",
        "# num_layers = 3\n",
        "# learning_rate = 0.00021291878584343366\n",
        "# early_stopping_patience = 10\n",
        "\n",
        "# Set training config\n",
        "root = \"./files\"\n",
        "experiment_name = \"seq2one-poc\"\n",
        "gpu = 0\n",
        "min_epochs, max_epochs = 3, 30\n",
        "\n",
        "# Set the hyperparameters\n",
        "batch_size = 64\n",
        "learning_rate = 0.001\n",
        "hidden_size = 128\n",
        "num_layers = 3\n",
        "early_stopping_patience = 15"
      ],
      "metadata": {
        "id": "is_2cE1-YdM0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Create the data module\n",
        "dm = FieldDataModule(root=root, batch_size=batch_size, workers=2)"
      ],
      "metadata": {
        "id": "mRg6EeSBYeRJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "e86jAftmyGIZ"
      },
      "outputs": [],
      "source": [
        "# Create the task with the sampled hyperparameters\n",
        "task = SequenceClassificationTask(input_size=512,\n",
        "                                  hidden_size=hidden_size,\n",
        "                                  num_layers=num_layers,\n",
        "                                  learning_rate=learning_rate)\n",
        "\n",
        "# Create a dedicated models' directory for saving the trial's best models\n",
        "models_path = Path(f\"./models/{experiment_name}/\")\n",
        "models_path.mkdir(parents=True, exist_ok=True)\n",
        "\n",
        "# Set the callbacks\n",
        "checkpoint_callback = ModelCheckpoint(\n",
        "    monitor=\"val_loss\",\n",
        "    dirpath=models_path,\n",
        "    filename=f\"model-{{epoch:02d}}-{{val_loss:.2f}}\",\n",
        "    save_top_k=1,\n",
        "    mode=\"min\",\n",
        ")\n",
        "early_stopping_callback = EarlyStopping(monitor=\"val_loss\",\n",
        "                                        mode=\"min\",\n",
        "                                        patience=early_stopping_patience)\n",
        "\n",
        "# Create a TensorBoard logger\n",
        "logger = TensorBoardLogger(\"./tb_logs\", name=experiment_name)\n",
        "\n",
        "# Trainer definition\n",
        "trainer = Trainer(\n",
        "    logger=logger,\n",
        "    accelerator='gpu',\n",
        "    devices=[gpu],\n",
        "    max_epochs=max_epochs,\n",
        "    min_epochs=min_epochs,\n",
        "    callbacks=[checkpoint_callback, early_stopping_callback],\n",
        "    precision=16\n",
        ")\n",
        "\n",
        "trainer.fit(model=task, datamodule=dm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5WJXRk-u0the"
      },
      "source": [
        "Let's check the best validation score:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mQDYjQMz0sv5"
      },
      "outputs": [],
      "source": [
        "checkpoint_callback.best_model_score.item()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K2IsZPbBUpdc"
      },
      "source": [
        "Let's get the prediction of the deep learning model on the test set:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "h25bbCpPUpdc"
      },
      "outputs": [],
      "source": [
        "# Load your model\n",
        "model = SequenceClassificationTask.load_from_checkpoint(trainer.checkpoint_callback.best_model_path,\n",
        "                                                        input_size=512,\n",
        "                                                        hidden_size=hidden_size)\n",
        "model.eval()\n",
        "model.freeze()\n",
        "\n",
        "# Get the validation data loader\n",
        "test_dl = dm.test_dataloader()\n",
        "\n",
        "# Predict\n",
        "all_logits = []\n",
        "y_tests = []\n",
        "\n",
        "with torch.no_grad():\n",
        "    for batch in test_dl:\n",
        "        inputs = batch['image']\n",
        "        y_test = batch['label']\n",
        "        logits = model(inputs)\n",
        "        all_logits.append(logits)\n",
        "        y_tests.append(y_test)\n",
        "\n",
        "# Concatenate all the results\n",
        "all_logits = torch.cat(all_logits, dim=0)\n",
        "y_test = torch.cat(y_tests, dim=0)\n",
        "\n",
        "# Get the probabilities\n",
        "y_test_hat = torch.nn.functional.softmax(all_logits, dim=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DUVyMobbUpdc"
      },
      "source": [
        "Let's calculate the final metrics:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vLVShFa6Updc"
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import log_loss\n",
        "\n",
        "# Get the arrays\n",
        "y_test_np = y_test.cpu().numpy()\n",
        "y_test_hat_np = y_test_hat.cpu().numpy()\n",
        "\n",
        "# Convert y_val to a binary label indicator format\n",
        "y_test_bin = label_binarize(y_test_np, classes=[0, 1, 2, 3, 4, 5, 6])\n",
        "\n",
        "cross_entropy = log_loss(y_test_bin, y_test_hat_np)\n",
        "print(\"Cross Entropy:\", cross_entropy)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "b2qvh17OwOcY"
      },
      "source": [
        "---\n",
        "\n",
        "## Conclusion\n",
        "\n",
        "In this notebook, we explored the exciting field of geospatial machine learning. We began with a brief introduction to geospatial data, both in terms of vector and raster data. We then defined our problem and acquired the dataset.\n",
        "\n",
        "After preparing our data, we trained a first machine learning model (`LightGBM`) and evaluated its performance against a human baseline. We then used our model to make predictions, allowing us to see the power and potential of applying machine learning techniques to geospatial data. We repeated the same steps using a Seq2One neural network.\n",
        "\n",
        "Finally, we expanded our horizons by discussing the use of more advanced deep learning models for solving similar problems. This exploration underlines the fact that the techniques and approaches we've covered here are just the beginning. There's a wealth of more complex and powerful tools available in the geospatial machine learning space, and we encourage you to continue exploring!\n",
        "\n",
        "Remember, the journey of learning never ends, and each step brings us closer to making meaningful contributions to the field. Happy learning!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y1MbhybTwOcY"
      },
      "source": [
        "\n",
        "## Resources\n",
        "\n",
        "### Tutorials\n",
        "\n",
        "- [Geospatial Primer](https://github.com/Akramz/geospatial-primer).\n",
        "- [Introduction to Geospatial Data](https://colab.research.google.com/drive/1-85h5tEB0AJYT8xQ5H1wtSnCafXuLTHo#scrollTo=JDT5jUmCiTH-).\n",
        "- [Geospatial Data Analysis](https://colab.research.google.com/drive/1Yfkm63OV3eCtR3IVB-4owi2DJgj2Wd84).\n",
        "- [Geospatial Deep learning: Getting started with TorchGeo](https://pytorch.org/blog/geospatial-deep-learning-with-torchgeo/).\n",
        "- [Automating GIS-processes Course]((https://autogis-site.readthedocs.io/en/latest/))\n",
        "- [Geospatial Data with Python: Shapely and Fiona](https://macwright.com/2012/10/31/gis-with-python-shapely-fiona.html)\n",
        "- [Introduction to Raster Data Processing in Open Source Python](https://www.earthdatascience.org/courses/use-data-open-source-python/intro-raster-data-python/raster-data-processing/).\n",
        "- [XArray fundamental](https://rabernat.github.io/research_computing_2018/xarray.html).\n",
        "- [XArray tutorials](https://github.com/xarray-contrib/xarray-tutorial).\n",
        "- [Visualization: contextily tutorial](https://geopandas.org/en/stable/gallery/plotting_basemap_background.html).\n",
        "\n",
        "\n",
        "### Libraries\n",
        "\n",
        "- [Shapely](https://github.com/shapely/shapely).\n",
        "- [GeoPandas](https://github.com/geopandas/geopandas).\n",
        "- [Contextily](https://github.com/geopandas/contextily).\n",
        "- [Rasterio](https://github.com/rasterio/rasterio).\n",
        "- [Xarray](https://github.com/pydata/xarray).\n",
        "- [RioXarray](https://github.com/corteva/rioxarray).\n",
        "- [TorchGeo](https://github.com/microsoft/torchgeo).\n",
        "\n",
        "### Credits\n",
        "\n",
        "- [CV4A Kenya Crop Type Competition (source dataset)](https://mlhub.earth/data/ref_african_crops_kenya_02).\n",
        "- [Related Zindi Competition](https://zindi.africa/competitions/iclr-workshop-challenge-2-radiant-earth-computer-vision-for-crop-recognition/).\n",
        "- Paper: Jin, Zhenong, et al. \"[Smallholder maize area and yield mapping at national scales with Google Earth Engine.](https://web.stanford.edu/~mburke/papers/jin%20et%20al%202019.pdf)\" Remote Sensing of Environment 228 (2019): 115-128.\n",
        "\n",
        "If you want to keep up with Geospatial machine learning research in Africa, join our [GeoAI Africa](https://www.linkedin.com/company/99894375/admin/feed/posts/) community!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dHFlWWJAh1Us"
      },
      "source": [
        "## Feedback\n",
        "\n",
        "Please provide feedback that we can use to improve our practicals in the future."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "FHraUWVbh3yF"
      },
      "outputs": [],
      "source": [
        "# @title Generate Feedback Form. (Run Cell)\n",
        "from IPython.display import HTML\n",
        "\n",
        "HTML(\n",
        "    \"\"\"\n",
        "<iframe\n",
        "\tsrc=\"https://forms.gle/Cg9aoa7czoZCYqxF7\",\n",
        "  width=\"80%\"\n",
        "\theight=\"1200px\" >\n",
        "\tLoading...\n",
        "</iframe>\n",
        "\"\"\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mI4Q-35Oh6NP"
      },
      "source": [
        "<img src=\"https://baobab.deeplearningindaba.com/static/media/indaba-logo-dark.d5a6196d.png\" width=\"50%\" />"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dWZrAx7UufpT"
      },
      "source": [
        "---"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}