geospatial_machine_learning.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/mgamzec/GeospatialData-Python/blob/main/geospatial_machine_learning.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TV85F3PEufpF"
      },
      "source": [
        "# Geospatial machine learning\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. **Data** acquisition and preprocessing: we will get the data and preprocess it for machine learning.\n",
        "4. **Model fitting**: we will fit a model to the data and conduct hyperparameter search.\n",
        "5. **Model evaluation**: we will evaluate the model on the test data.\n",
        "6. **Inference**: we will predict the output for the test data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hjko0VXhvXEF"
      },
      "source": [
        "Before we start, let's install the necessary libraries:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p3LcpCLYvbjV"
      },
      "outputs": [],
      "source": [
        "%pip install rioxarray -q\n",
        "%pip install shap -q\n",
        "%pip install contextily -q"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oo4IQkXa2Bg0"
      },
      "source": [
        "Next, we need to download the necessary files to be used in this practical:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6WjhXrfx1_1-"
      },
      "outputs": [],
      "source": [
        "from pathlib import Path\n",
        "import requests\n",
        "import shutil\n",
        "import os\n",
        "\n",
        "# Set the URLs of the files\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",
        "\n",
        "def download_file(url, directory, file_name):\n",
        "\n",
        "    print(f\"Saving `{file_name}` to `{directory}/`\")\n",
        "\n",
        "    # Create directory if it doesn't exist\n",
        "    if not os.path.exists(directory):\n",
        "        os.makedirs(directory)\n",
        "\n",
        "    file_path = os.path.join(directory, file_name)\n",
        "\n",
        "    # Send a HTTP request to the URL of the file\n",
        "    response = requests.get(url, stream=True)\n",
        "\n",
        "    # Throw an error for bad status codes\n",
        "    response.raise_for_status()\n",
        "\n",
        "    # Write the file\n",
        "    with open(file_path, 'wb') as handle:\n",
        "        for block in response.iter_content(1024):\n",
        "            handle.write(block)\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(raster_url, data_dir, \"elevation.tiff\")\n",
        "download_file(crop_url, data_dir, \"df.feather\")\n",
        "download_file(bikes_url, data_dir, \"paris_bike_stations_mercator.gpkg\")\n",
        "download_file(paris_districts_url, data_dir, \"paris_districts_utm.geojson\")"
      ]
    },
    {
      "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 style=\"width:33%;\" 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": [
        "## Vector\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": "wcI2F3vGufpH"
      },
      "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": "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",
        "polygon_union"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4pRIn7kVufpI"
      },
      "outputs": [],
      "source": [
        "# Intersection of two polygons\n",
        "polygon_intersection = polygon1.intersection(polygon2)\n",
        "polygon_intersection"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xHr77oRxufpI"
      },
      "outputs": [],
      "source": [
        "# The difference\n",
        "polygon1.difference(polygon2)"
      ]
    },
    {
      "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": "QYbg03ZLze9T"
      },
      "source": [
        "### Exercises"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MtABnRFwze9T"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "We will start with exploring the bicycle station dataset (available as a GeoPackage file: `data/paris_bike_stations_mercator.gpkg`)\n",
        "    \n",
        "* Read the stations datasets into a GeoDataFrame called `stations`.\n",
        "* Check the type of the returned object\n",
        "* Check the first rows of the dataframes. What kind of geometries does this datasets contain?\n",
        "* How many features are there in the dataset?\n",
        "    \n",
        "<details><summary>Hints</summary>\n",
        "\n",
        "* Use `type(..)` to check any Python object type\n",
        "* The `geopandas.read_file()` function can read different geospatial file formats. You pass the file name as first argument.\n",
        "* Use the `.shape` attribute to get the number of features\n",
        "\n",
        "</details>\n",
        "    \n",
        "    \n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OwJLd72Xze9T"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "* Make a quick plot of the `stations` dataset.\n",
        "* Make the plot a bit larger by setting the figure size to (12, 6) (hint: the `plot` method accepts a `figsize` keyword).\n",
        "\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ey50qJOtze9T"
      },
      "source": [
        "A plot with points can be hard to interpret without any spatial context. Therefore, we will learn how to add a background map.\n",
        "\n",
        "We are going to make use of the [contextily](https://github.com/darribas/contextily) package. the `add_basemap()` function of this package makes it easy to add a background web map to our plot. We begin by plotting our data, then pass the matplotlib axes object to the `add_basemap()` function. `contextily` will then download the web tiles needed for the geographical extent of the plot.\n",
        "\n",
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "* Import `contextily`.\n",
        "* Re-do the figure of the previous exercise: make a plot of all the points in `stations`, but assign the result to an `ax` variable.\n",
        "* Set the marker size equal to 5 to reduce the size of the points (use the `markersize` keyword of the `plot()` method for this).\n",
        "* Use the `add_basemap()` function of `contextily` to add a background map: the first argument is the matplotlib axes object `ax`.\n",
        "\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KkhxCNKXze9T"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "* Make a histogram showing the distribution of the number of bike stands in the stations.\n",
        "\n",
        "<details>\n",
        "  <summary>Hints</summary>\n",
        "\n",
        "* Selecting a column can be done with the square brackets: `df['col_name']`\n",
        "* Single columns have a `hist()` method to plot a histogram of its values.\n",
        "    \n",
        "</details>\n",
        "    \n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jC7xhcZtze9T"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "Let's now visualize where the available bikes are actually stationed:\n",
        "    \n",
        "* Make a plot of the `stations` dataset (also with a (12, 6) figsize).\n",
        "* Use the `'available_bikes'` columns to determine the color of the points. For this, use the `column=` keyword.\n",
        "* Use the `legend=True` keyword to show a color bar.\n",
        "\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "feybB8-Vze9T"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "Next, we will explore the dataset on the administrative districts of Paris (available as a GeoJSON file: \"data/paris_districts_utm.geojson\")\n",
        "\n",
        "* Read the dataset into a GeoDataFrame called `districts`.\n",
        "* Check the first rows of the dataframe. What kind of geometries does this dataset contain?\n",
        "* How many features are there in the dataset? (hint: use the `.shape` attribute)\n",
        "* Make a quick plot of the `districts` dataset (set the figure size to (12, 6)).\n",
        "    \n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MxKJzRQVze9U"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "    \n",
        "What are the largest districts (biggest area)?\n",
        "\n",
        "* Calculate the area of each district.\n",
        "* Add this area as a new column to the `districts` dataframe.\n",
        "* Sort the dataframe by this area column for largest to smallest values (descending).\n",
        "\n",
        "<details><summary>Hints</summary>\n",
        "\n",
        "* Adding a column can be done by assigning values to a column using the same square brackets syntax: `df['new_col'] = values`\n",
        "* To sort the rows of a DataFrame, use the `sort_values()` method, specifying the colum to sort on with the `by='col_name'` keyword. Check the help of this method to see how to sort ascending or descending.\n",
        "\n",
        "</details>\n",
        "\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IroekjJpze9U"
      },
      "source": [
        "<div class=\"alert alert-success\">\n",
        "\n",
        "**EXERCISE**:\n",
        "\n",
        "* Add a column `'population_density'` representing the number of inhabitants per squared kilometer (Note: The area is given in squared meter, so you will need to multiply the result with `10**6`).\n",
        "* Plot the districts using the `'population_density'` to color the polygons. For this, use the `column=` keyword.\n",
        "* Use the `legend=True` keyword to show a color bar.\n",
        "\n",
        "</div>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LuQJNMJBufpJ"
      },
      "source": [
        "## Raster\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 the field by trying to solve a problem.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dvmQ-cJkufpJ"
      },
      "source": [
        "# Case study: *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": [
        "## Problem Scoping\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 style=\"width:50%;\" 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": "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": "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": "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": "code",
      "execution_count": null,
      "metadata": {
        "id": "mmRMY0PrufpR"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "\n",
        "np.any(df.isna())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LWEgNnHQufpR"
      },
      "source": [
        "Let's conduct the train/validation/test split:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "k_2IbyajufpR"
      },
      "outputs": [],
      "source": [
        "# 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.imgur.com/Brd7Jqe.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": "code",
      "execution_count": null,
      "metadata": {
        "id": "yQMpP8zAufpS"
      },
      "outputs": [],
      "source": [
        "from sklearn.metrics import log_loss\n",
        "\n",
        "# Calculate cross-entropy\n",
        "cross_entropy = log_loss(y_val, y_val_hat)\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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0bJl0QjRufpS"
      },
      "source": [
        "Let's implement the feature engineering function that would add additional vegetation indices of interest:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0Fc_Fbk_ufpS"
      },
      "outputs": [],
      "source": [
        "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": {
        "id": "TBXtxAafufpS"
      },
      "outputs": [],
      "source": [
        "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": {
        "id": "YSdcas3mufpS"
      },
      "outputs": [],
      "source": [
        "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": "code",
      "execution_count": null,
      "metadata": {
        "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",
        "\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",
        "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": [
        "---\n",
        "\n",
        "## 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": "b2qvh17OwOcY"
      },
      "source": [
        "---\n",
        "\n",
        "# Conclusion\n",
        "\n",
        "In this notebook, we explored the broad and exciting field of geospatial machine learning. We began with a brief introduction to geospatial data, both in terms of vector and raster information. We then defined our problem and proceeded to gather and preprocess the necessary data to address it.\n",
        "\n",
        "After preparing our data, we fitted a machine learning model and evaluated its performance against the 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.\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",
        "\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",
        "### Geospatial 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",
        "- [Gilles Q. HACHEME](https://www.linkedin.com/in/gilles-q-hacheme-a0956ab7/).\n",
        "- [Aisha Alaagib](https://sa.linkedin.com/in/aishaalaagib).\n",
        "- [Akram Zaytar](https://www.microsoft.com/en-us/research/people/akramzaytar/).\n",
        "- [Girmaw Abebe Tadesse](https://www.microsoft.com/en-us/research/people/gtadesse/)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dWZrAx7UufpT"
      },
      "source": [
        "---"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "collapsed_sections": [
        "QYbg03ZLze9T"
      ],
      "include_colab_link": true
    },
    "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.12"
    },
    "orig_nbformat": 4
  },
  "nbformat": 4,
  "nbformat_minor": 0
}