{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "### Image Processing with Python\n",
        "\n",
        "Now that you have downloaded some Planetscope imagery, what can you do with it? This workshop includes an interactive overview of basic image processing with rasters and vectors in Python with spatial libraries like rasterio and geopandas.\n",
        "\n",
        "\n",
        "#### Background\n",
        "Lake Lagunita on Stanford Campus is not a permanantly filled lake unless high rainfall fills it by diversion from San Francisquito Creek. Let's use Python libraries to inspect this area and the timing of a flooding event in the winter of 2022/2023 ([link to Explorer](https://www.planet.com/explorer/?s=VQxuExgIQwOCE3LaspaYgQ)) with Planetscope imagery."
      ],
      "metadata": {
        "id": "bDmU8jay2QEo"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ZDW1q_010z-c",
        "outputId": "541e85da-f64b-4059-a1d3-9a0bf4f80a1a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Import spatial and plotting libraries\n",
        "\n",
        "Below we will import a handful of raster and vector processing libraries.  Notable processing libraries include **rasterio** and **geopandas**.\n",
        "\n",
        "**Rasterio** is a Python library for working with geospatial raster data. Rasterio makes it easy to open, read, and write raster datasets in various formats, such as GeoTIFF and JPEG, and perform common geospatial operations.\n",
        "\n",
        "**GeoPandas** is an open-source Python library that extends the capabilities of the popular Pandas library to handle geospatial data. Many geopandas datasets will take the form of a dataframe with a \"geometry\" column that allows for spatial operations like plotting, joining, or clipping and aggregating raster data. It simplifies the manipulation and analysis of geographic datasets, allowing users to work with both vector and raster data within a unified framework."
      ],
      "metadata": {
        "id": "V2pxTpi92U7S"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Install Libraries"
      ],
      "metadata": {
        "id": "qR0d9bcp2j2q"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install rasterio"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QrOp6qu92g9m",
        "outputId": "796487e9-d5d2-4183-eb7a-04c78a8a40f6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting rasterio\n",
            "  Downloading rasterio-1.3.9-cp310-cp310-manylinux2014_x86_64.whl (20.6 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m20.6/20.6 MB\u001b[0m \u001b[31m46.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting affine (from rasterio)\n",
            "  Downloading affine-2.4.0-py3-none-any.whl (15 kB)\n",
            "Requirement already satisfied: attrs in /usr/local/lib/python3.10/dist-packages (from rasterio) (23.1.0)\n",
            "Requirement already satisfied: certifi in /usr/local/lib/python3.10/dist-packages (from rasterio) (2023.7.22)\n",
            "Requirement already satisfied: click>=4.0 in /usr/local/lib/python3.10/dist-packages (from rasterio) (8.1.7)\n",
            "Requirement already satisfied: cligj>=0.5 in /usr/local/lib/python3.10/dist-packages (from rasterio) (0.7.2)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.23.5)\n",
            "Collecting snuggs>=1.4.1 (from rasterio)\n",
            "  Downloading snuggs-1.4.7-py3-none-any.whl (5.4 kB)\n",
            "Requirement already satisfied: click-plugins in /usr/local/lib/python3.10/dist-packages (from rasterio) (1.1.1)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.10/dist-packages (from rasterio) (67.7.2)\n",
            "Requirement already satisfied: pyparsing>=2.1.6 in /usr/local/lib/python3.10/dist-packages (from snuggs>=1.4.1->rasterio) (3.1.1)\n",
            "Installing collected packages: snuggs, affine, rasterio\n",
            "Successfully installed affine-2.4.0 rasterio-1.3.9 snuggs-1.4.7\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Import Libraries"
      ],
      "metadata": {
        "id": "ITE_I5Z92pub"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pfvt_UmA1CIB"
      },
      "outputs": [],
      "source": [
        "# Import libraries\n",
        "import os\n",
        "import glob\n",
        "\n",
        "\n",
        "import geopandas as gpd\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import rasterio\n",
        "from rasterio import features\n",
        "from skimage import exposure\n",
        "\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\")\n",
        "\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "data_dir = \"/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/\"\n",
        "glob.glob(os.path.join(data_dir, \"*\"))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "MskQ5XfN21Vb",
        "outputId": "6306a94e-187b-4832-84d2-90fb9524ade6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/geoms',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data']"
            ]
          },
          "metadata": {},
          "execution_count": 9
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "glob.glob(os.path.join(data_dir, \"data\", \"*composite.tif\"))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QL2w-JCR2_Pn",
        "outputId": "ad7ae5de-a8cc-4752-f2ed-314ab79678ee"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-02-15_strip_6290924_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-30_strip_6682704_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-18_strip_6153179_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-30_strip_6743722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6458341_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-16_strip_6717370_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-25_strip_6384901_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080437_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-26_strip_6104067_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-30_strip_6615371_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-14_strip_6647556_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-30_strip_6250675_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-02_strip_6551855_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325614_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6457722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437404_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-17_strip_6219068_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080597_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-27_strip_6800033_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325499_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437085_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-13_strip_6770796_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-05-17_strip_6515202_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-02_strip_6117199_composite.tif']"
            ]
          },
          "metadata": {},
          "execution_count": 11
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "files = glob.glob(os.path.join(data_dir, \"data\", \"*composite.tif\"))"
      ],
      "metadata": {
        "id": "LcVvTSxR2Bqt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's select a single file to explore its contents. We'll revisit the rest later."
      ],
      "metadata": {
        "id": "2elZ2aBc3H8d"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "test_file = os.path.join(data_dir, \"data\", \"2023-04-15_strip_6437085_composite.tif\")\n",
        "os.path.exists(test_file)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DsGbaKdV3KvW",
        "outputId": "0e7eeb78-5c17-4f39-a695-992b609fca71"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "True"
            ]
          },
          "metadata": {},
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Read the file with Rasterio\n",
        "with rasterio.open(test_file) as src:\n",
        "    print(src.meta)\n",
        "    # From the open src image, let's save the array as a variable.\n",
        "    arr = src.read()/10000"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iQfviNOK3Kxy",
        "outputId": "63288d07-fb7f-4dd6-f4ed-9a45fe63c418"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': 0.0, 'width': 438, 'height': 343, 'count': 8, 'crs': CRS.from_epsg(32610), 'transform': Affine(3.0, 0.0, 572298.0,\n",
            "       0.0, -3.0, 4142520.0)}\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "arr.shape"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Ja7yWeeq3K0C",
        "outputId": "36a3beef-17f6-4e3d-a4f1-29ed2b427df5"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(8, 343, 438)"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "As we read the array, we must divide reflectance by 10,000 because Planet imagery stores reflectance as a value scaled by 10,000 to adhere to the required data type (uint16) while maintaining decimal precision (described in Atmospheric Corrections on page 18 of the [Product Spec](https://assets.planet.com/docs/Planet_PSScene_Imagery_Product_Spec_letter_screen.pdf))."
      ],
      "metadata": {
        "id": "JL911m8J3Woq"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Plotting individual bands in greyscale\n",
        "\n",
        "Next, we will plot individual bands in greyscale to check absorption features of the land cover types we do know."
      ],
      "metadata": {
        "id": "trzgx5I23dyC"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Order of the bands for 8-band planetscope imagery is from lowest to highest wavelength\n",
        "band_order = [\"coastal_blue\", \"blue\", \"green_I\", \"green\", \"yellow\", \"red\", \"red_edge\", \"nir\"]"
      ],
      "metadata": {
        "id": "6OYJ01Wz3K2i"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# With the list above, we can identify the index of a color by name and use it to access that layer in the array.\n",
        "band_order.index(\"red\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ZFKjhQV93XyJ",
        "outputId": "88981374-f933-4855-abe0-192993dd2729"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "5"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "arr[5].shape"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VDde38MN3Xz8",
        "outputId": "8bbdd764-372c-4ea1-9923-d8aba1d55bf0"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(343, 438)"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We'll use the ordered list of bands to loop over and plot each band in greyscale, and title the plot according to the color name."
      ],
      "metadata": {
        "id": "KaWWEg6e3k4K"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Define the layout of our figure\n",
        "nrow = 2\n",
        "ncol = 4\n",
        "f, axes = plt.subplots(nrow, ncol, figsize=(5*ncol, 3.5*nrow))\n",
        "\n",
        "# Loop over every band name and axis available in the plot together\n",
        "for band_name, ax in zip(band_order, axes.flatten()):\n",
        "    # Extract the band from our array\n",
        "    band_index = band_order.index(band_name)\n",
        "    band_arr = arr[band_index]\n",
        "\n",
        "    # Plot in greyscale to the subplot\n",
        "    im = ax.imshow(band_arr, cmap=\"Greys_r\")\n",
        "    f.colorbar(im, ax=ax)\n",
        "\n",
        "    # Set the title for the subplot\n",
        "    ax.set_title(band_name)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 404
        },
        "id": "yeu5rkzZ3X2M",
        "outputId": "54b6f3e5-04b3-4051-ca8e-2fcc525e23e3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 2000x700 with 16 Axes>"
            ],
            "image/png": "[encoded data, 618252 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We've plotted all of the bands!\n",
        "\n",
        "Note reflectance and absorption of different wavelengths by features we know:\n",
        "- water absorbing most wavelengths more than other features\n",
        "- buildings reflecting most wavelengths, so they're bright in most subplots\n",
        "- grassy field absorb blues and reds, and in increasing intensity reflect yellow, greens, red edge, and NIR"
      ],
      "metadata": {
        "id": "oQTrUR6L3oNj"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Creating a function to plot an RGB\n",
        "\n",
        "Let's plot an RGB image next!\n",
        "\n",
        "First, we can test how best to stretch and plot the image before writing a function."
      ],
      "metadata": {
        "id": "TWzMPJ4X3rIp"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Save R, G, and B as variables\n",
        "red = arr[band_order.index(\"red\")]\n",
        "green = arr[band_order.index(\"green\")]\n",
        "blue = arr[band_order.index(\"blue\")]\n",
        "\n",
        "# Stack them into a single array\n",
        "img = np.stack([red, green, blue], axis=2)"
      ],
      "metadata": {
        "id": "yrTsLkIr3qUU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "When plotting an RGB image you are plotting 3 bands together and decisions must be made about how best to balance the color contrast. This is unlike individual bands plotted above where matplotlib can dynamically stretch to fit the range of individual bands.\n",
        "\n",
        "One common approach to color balancing an RGB is rescaling all bands to a range defined by their 2nd and 98th percentiles."
      ],
      "metadata": {
        "id": "CW8hy3Z83udC"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Rescale surface reflectance values for plotting with color balance\n",
        "p2, p98 = np.percentile(img, (2, 98))\n",
        "img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))"
      ],
      "metadata": {
        "id": "WAT8TNMw3X4U"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "f, ax = plt.subplots()\n",
        "ax.hist(red[red!=0], label=\"red band before rescaling\", histtype=\"step\", bins = 20);\n",
        "ax.hist(img_rescale[:, :, 0][red!=0], label=\"red band after rescaling\", histtype=\"step\")\n",
        "ax.legend();"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 430
        },
        "id": "bqoHAi9N3X6c",
        "outputId": "cde62611-04c2-4644-d113-95124ef76d27"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "[encoded data, 22540 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "f, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
        "axes[0].imshow(img)\n",
        "axes[0].set_title(\"RGB before rescaling\")\n",
        "axes[1].imshow(img_rescale)\n",
        "axes[1].set_title(\"RGB after rescaling\");"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 303
        },
        "id": "Yx44NnTw3X8u",
        "outputId": "c82e3a8c-8188-465a-9829-c6e8dc9c12e3"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x300 with 2 Axes>"
            ],
            "image/png": "[encoded data, 311004 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can define a function, and even play with different band orders:"
      ],
      "metadata": {
        "id": "LFGyailG31qL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def arr_to_rgb(arr, use_bands=[\"red\", \"green\", \"blue\"]):\n",
        "    band_order = [\n",
        "        \"coastal_blue\", \"blue\", \"green_I\", \"green\",\n",
        "        \"yellow\", \"red\", \"red_edge\", \"nir\"\n",
        "    ]\n",
        "    b1 = arr[band_order.index(use_bands[0])]\n",
        "    b2 = arr[band_order.index(use_bands[1])]\n",
        "    b3 = arr[band_order.index(use_bands[2])]\n",
        "\n",
        "    img = np.stack([b1, b2, b3], axis=2)\n",
        "\n",
        "    p2, p98 = np.percentile(img, (2, 98))\n",
        "    img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))\n",
        "    return img_rescale"
      ],
      "metadata": {
        "id": "i9IVEuFP3FiA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "f, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
        "axes[0].imshow(arr_to_rgb(arr))\n",
        "axes[0].set_title(\"Red-Green-Blue (RGB)\\nvisual inspection\")\n",
        "axes[1].imshow(arr_to_rgb(arr,  use_bands=[\"nir\", \"red\", \"green\"]))\n",
        "axes[1].set_title(\"False Color Composite (NIR-R-G)\\nidentifying vegetation\");"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 324
        },
        "id": "rJ2yOPkr33Lh",
        "outputId": "3e075c83-4b42-4ead-9a24-8a8f9cebc18e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x300 with 2 Axes>"
            ],
            "image/png": "[encoded data, 375720 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's turn these steps into a function, and loop over every image to get a better view of when water began to flood Lake Lagunita."
      ],
      "metadata": {
        "id": "texxy65r35xU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Plotting several RGB images\n",
        "\n",
        "First, we'll sort our files in alphabetical order.  This will allow us to iterate over the list in an order that is chronological."
      ],
      "metadata": {
        "id": "YGuUKEmA37vV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# These are out of order\n",
        "files"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8FVR7rpe34SP",
        "outputId": "fa2ecbd7-d0c1-4076-ac5b-fd4a51601bb2"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-02-15_strip_6290924_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-30_strip_6682704_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-18_strip_6153179_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-30_strip_6743722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6458341_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-16_strip_6717370_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-25_strip_6384901_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080437_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-26_strip_6104067_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-30_strip_6615371_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-14_strip_6647556_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-30_strip_6250675_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-02_strip_6551855_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325614_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6457722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437404_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-17_strip_6219068_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080597_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-27_strip_6800033_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325499_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437085_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-13_strip_6770796_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-05-17_strip_6515202_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-02_strip_6117199_composite.tif']"
            ]
          },
          "metadata": {},
          "execution_count": 28
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# These are now in alphabetical order\n",
        "files.sort()\n",
        "files"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "S7MkDMwY39Op",
        "outputId": "2b095dac-b439-4e42-c8eb-2731a3ba8640"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "['/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080437_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080597_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-26_strip_6104067_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-02_strip_6117199_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-18_strip_6153179_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-17_strip_6219068_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-30_strip_6250675_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-02-15_strip_6290924_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325499_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325614_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-25_strip_6384901_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437085_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437404_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6457722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6458341_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-05-17_strip_6515202_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-02_strip_6551855_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-30_strip_6615371_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-14_strip_6647556_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-30_strip_6682704_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-16_strip_6717370_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-30_strip_6743722_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-13_strip_6770796_composite.tif',\n",
              " '/content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-27_strip_6800033_composite.tif']"
            ]
          },
          "metadata": {},
          "execution_count": 29
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# We can decide on how many rows and columns would display a total of 24 images\n",
        "len(files)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "o6YoNAdx4HA4",
        "outputId": "a9b89677-d296-4209-ca06-09b063ffb3f1"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "24"
            ]
          },
          "metadata": {},
          "execution_count": 31
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "nrow = 4\n",
        "ncol = 6\n",
        "\n",
        "f, axes = plt.subplots(nrow, ncol, figsize=(2*ncol, 2*nrow))\n",
        "for file, ax in zip(files, axes.flatten()):\n",
        "    with rasterio.open(file) as src:\n",
        "        arr = src.read()/10000\n",
        "    rgb_im = arr_to_rgb(arr)\n",
        "    #rgb_im = arr_to_rgb(arr,  use_bands=[\"nir\", \"red\", \"green\"])\n",
        "    ax.imshow(rgb_im)\n",
        "\n",
        "    date = os.path.basename(file).split(\"_\")[0]\n",
        "    ax.set_title(date)\n",
        "\n",
        "for ax in axes.flatten():\n",
        "    ax.axis(\"off\")\n",
        "plt.tight_layout()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 715
        },
        "id": "gmC5sUHE4JyD",
        "outputId": "ca3aa20c-a3cd-4733-ac16-0ce2b036b625"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1200x800 with 24 Axes>"
            ],
            "image/png": "[encoded data, 2107824 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can visually assess these images and guess that flooding began at the end of December."
      ],
      "metadata": {
        "id": "9Z3ziaNy4VDe"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Raster Analysis with Band Math: Calculating an Index\n",
        "\n",
        "Next, we'll calculate a spectral index.  Spectral indices are mathematical algorithms used in remote sensing to highlight specific features or characteristics of the Earth's surface.\n",
        "\n",
        "The Normalized Difference Water Index (NDWI) is derived from the Near-Infrared (NIR) and Green (G) channels. This formula is used to detect water bodies:\n",
        "\n",
        "NDWI = (G-NIR)/(G+NIR)\n",
        "\n",
        "Higher NDWI values typically indicate the presence of water, while lower values correspond to land surfaces.  The NDWI values correspond to the following ranges:\n",
        "\n",
        "| Range      | Description                  |\n",
        "|------------|------------------------------|\n",
        "| 0.2 – 1   | Water surface      |\n",
        "| 0.0 – 0.2   | Flooding, humidity      |\n",
        "| -0.3 – 0.0    | Moderate drought, non-aqueous surfaces      |\n",
        "| -1 – -0.3    | Drought, non-aqueous surfaces      |"
      ],
      "metadata": {
        "id": "E_0nvf1e4Vju"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Read the test file in with rasterio, and pull Green and NIR\n",
        "with rasterio.open(test_file) as src:\n",
        "    arr = src.read()/10000\n",
        "green = arr[band_order.index(\"green\")]\n",
        "nir = arr[band_order.index(\"nir\")]"
      ],
      "metadata": {
        "id": "ouwWTBUF4Sw-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "ndwi = (green - nir)/(green + nir)"
      ],
      "metadata": {
        "id": "Zuxw_c1l4aEH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "f, axes = plt.subplots(1, 2, figsize=(8, 3))\n",
        "\n",
        "# plot the RGB in the first subplot\n",
        "rgb = arr_to_rgb(arr)\n",
        "axes[0].imshow(rgb)\n",
        "\n",
        "# Plot NDWI in the second subplot, high values will indicate the presence of water.\n",
        "im = axes[1].imshow(ndwi)\n",
        "f.colorbar(im, ax=axes[1])\n",
        "\n",
        "plt.tight_layout();"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 317
        },
        "id": "5uovXY_t4bqS",
        "outputId": "dacfb18a-6435-4887-aafe-dd1ebdea8b83"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 800x300 with 3 Axes>"
            ],
            "image/png": "[encoded data, 381980 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Spectral indices like NDWI help us to objectively declare an area we're looking at is likely \"water\""
      ],
      "metadata": {
        "id": "K5ZhIJr84ebN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "f, ax = plt.subplots()\n",
        "im = plt.imshow(ndwi > 0)\n",
        "f.colorbar(im, ax=ax)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 438
        },
        "id": "ue2sQOB04cqA",
        "outputId": "22f6647a-2f78-4101-9609-2ecffb64a963"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.colorbar.Colorbar at 0x78c3b43a70a0>"
            ]
          },
          "metadata": {},
          "execution_count": 36
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 2 Axes>"
            ],
            "image/png": "[encoded data, 32504 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Vector Data: How does NDWI look inside of our Lake Lagunita Boundary?\n",
        "\n",
        "Next, we'll read a vector dataset containing boundaries around Lake Lagunita.\n",
        "\n",
        "In this section, we will load a geojson into a geodataframe."
      ],
      "metadata": {
        "id": "AhzwV2Ez4k48"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# boundaries which surround the lake\n",
        "geojson_filename = os.path.join(data_dir, 'geoms', 'lake-lagunita.geojson')\n",
        "gdf = gpd.read_file(geojson_filename)\n",
        "gdf"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 81
        },
        "id": "UjmAbS9M4h5j",
        "outputId": "eec95002-a43b-4d64-83f9-68d76327f504"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "                                            geometry\n",
              "0  POLYGON ((-122.17712 37.42450, -122.17800 37.4..."
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-c8430e9d-b9d8-48b0-bd98-a0e9bcfedb06\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>geometry</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>POLYGON ((-122.17712 37.42450, -122.17800 37.4...</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-c8430e9d-b9d8-48b0-bd98-a0e9bcfedb06')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-c8430e9d-b9d8-48b0-bd98-a0e9bcfedb06 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-c8430e9d-b9d8-48b0-bd98-a0e9bcfedb06');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ]
          },
          "metadata": {},
          "execution_count": 38
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "A geodataframe is like a pandas dataframe, but contains a standardized spatial attribute called \"geometry\" which can be used to map a vector dataset spatially or to extract data from spatial geotiffs."
      ],
      "metadata": {
        "id": "lSzxAGQx4r2X"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "gdf.plot()\n",
        "plt.title(\"Mapping Lake Lagunita using the GeoDataFrame's 'geometry' attribute\\n\");"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 506
        },
        "id": "Vyeypsel4mPX",
        "outputId": "f32f1f24-ba27-4cfa-def9-2c2c09bdc54a"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "[encoded data, 43164 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Masking raster data with vector data\n",
        "To join these two, we can create an array out of our vector dataset that is True where Lake Lagunita is present.\n",
        "\n",
        "Converting the vector dataset to a raster requires using the CRS and Tranform of the source raster."
      ],
      "metadata": {
        "id": "7kRJECFC4vCz"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "with rasterio.open(test_file) as src:\n",
        "    # reading metadata\n",
        "    metadata = src.meta\n",
        "\n",
        "# Cast to the same CRS as the rasterio image\n",
        "gdf_crs = gdf.to_crs(metadata[\"crs\"])\n",
        "\n",
        "# Create an array with a geometry mask\n",
        "aoi_mask = features.geometry_mask(gdf_crs.geometry, ndwi.shape, metadata[\"transform\"])\n",
        "# Flipping the mask so that True values are over the field\n",
        "lake_mask = ~aoi_mask"
      ],
      "metadata": {
        "id": "m1lD0ksj4yJ3"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "f, axes = plt.subplots(1, 2, figsize=(10, 5))\n",
        "axes[0].imshow(ndwi)\n",
        "axes[0].set_title(\"UDM2: Clear Data\")\n",
        "axes[1].imshow(lake_mask)\n",
        "axes[1].set_title(\"Lake Lag Location\");"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 358
        },
        "id": "j_6xorCN48aW",
        "outputId": "04fe6936-ee0b-4914-faae-fcc7ea815eb4"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 1000x500 with 2 Axes>"
            ],
            "image/png": "[encoded data, 225028 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Yellow are Trues in the right image.\n",
        "\n",
        "Below, we can use the lake_mask to mask values in the NDWI array"
      ],
      "metadata": {
        "id": "vp8C9Jm_4-R-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# all NDWI values in the lake\n",
        "ndwi[lake_mask]"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "H0Go2NYb4tHP",
        "outputId": "e352a056-e2f1-4440-e90f-76226a023262"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([-0.59990242, -0.60466849, -0.58866684, ..., -0.49433249,\n",
              "       -0.52315914, -0.52229498])"
            ]
          },
          "metadata": {},
          "execution_count": 43
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "The median of the masked values will be the median for the lake in the NDWI array."
      ],
      "metadata": {
        "id": "7kRJHhKd5C2G"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "median_ndwi_in_lake = ndwi[lake_mask].mean()\n",
        "print(f\"The median NDWI in measured in the lake is {median_ndwi_in_lake}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YMr8eTVI5CBt",
        "outputId": "a05ac6a8-f92b-4ea8-b54e-8d9aa4d85ec6"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The median NDWI in measured in the lake is -0.043109648565362306\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "median_ndwi_in_lake = ndwi[~np.isnan(ndwi) * ~lake_mask].mean()\n",
        "print(f\"The median NDWI in measured outside of the lake is {median_ndwi_in_lake}\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iL1XVLlf4_wb",
        "outputId": "d2b7ba25-2985-4716-daf7-fa9c846c288b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The median NDWI in measured outside of the lake is -0.5508857888236413\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Calculating a timeseries with an area of interest\n",
        "\n",
        "Next, we'll scale up this workflow to run over every image and plot a timeseries of NDWI in Lake Lagunita."
      ],
      "metadata": {
        "id": "eymnB97Z5G5L"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def calc_ndwi(image_filepath):\n",
        "    \"\"\"\n",
        "    Wrapper to calculate NDWI for an image\n",
        "    \"\"\"\n",
        "    with rasterio.open(image_filepath) as src:\n",
        "        arr = src.read()/10000\n",
        "    green = arr[band_order.index(\"green\")]\n",
        "    nir = arr[band_order.index(\"nir\")]\n",
        "    ndwi = (green - nir)/(green + nir)\n",
        "\n",
        "    return ndwi\n",
        "\n",
        "def median_ndwi_in_aoi(image_filepath, geojson_filepath):\n",
        "    \"\"\"\n",
        "    Wrapper to report median NDWI within a geojson geometry for an input image.\n",
        "    \"\"\"\n",
        "\n",
        "    ndwi = calc_ndwi(image_filepath)\n",
        "    with rasterio.open(image_filepath) as src:\n",
        "        metadata = src.meta\n",
        "\n",
        "    # Create an array with True values over the field location\n",
        "    gdf = gpd.read_file(geojson_filepath)\n",
        "    aoi_mask = features.geometry_mask(\n",
        "            gdf.to_crs(metadata[\"crs\"]).geometry, ndwi.shape, metadata[\"transform\"])\n",
        "\n",
        "    median_ndwi_in_lake = ndwi[lake_mask].mean()\n",
        "\n",
        "    return median_ndwi_in_lake"
      ],
      "metadata": {
        "id": "5TUDBAb_5Soe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Below, we'll create empty lists for `dates` and `medians` which we'll fill with values as we loop through each file."
      ],
      "metadata": {
        "id": "ZN7Aoz2p6tgV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "dates = []\n",
        "medians = []\n",
        "for im_filepath in files:\n",
        "    # From the filename, we'll pull dates.\n",
        "    filename = os.path.basename(im_filepath)\n",
        "    date = filename.split(\"_\")[0]\n",
        "    print(f\"Original filename: {im_filepath}\")\n",
        "    print(f\"Date: {date}\")\n",
        "    dates.append(date)\n",
        "\n",
        "    # Calculate medians\n",
        "    median = median_ndwi_in_aoi(im_filepath, geojson_filename)\n",
        "    medians.append(median)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nhurQ2-i6rge",
        "outputId": "acb0a2cf-d828-4197-d4e2-9f31ed3a86cc"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080437_composite.tif\n",
            "Date: 2022-11-15\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-15_strip_6080597_composite.tif\n",
            "Date: 2022-11-15\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-11-26_strip_6104067_composite.tif\n",
            "Date: 2022-11-26\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-02_strip_6117199_composite.tif\n",
            "Date: 2022-12-02\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2022-12-18_strip_6153179_composite.tif\n",
            "Date: 2022-12-18\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-17_strip_6219068_composite.tif\n",
            "Date: 2023-01-17\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-01-30_strip_6250675_composite.tif\n",
            "Date: 2023-01-30\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-02-15_strip_6290924_composite.tif\n",
            "Date: 2023-02-15\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325499_composite.tif\n",
            "Date: 2023-03-02\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-02_strip_6325614_composite.tif\n",
            "Date: 2023-03-02\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-03-25_strip_6384901_composite.tif\n",
            "Date: 2023-03-25\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437085_composite.tif\n",
            "Date: 2023-04-15\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-15_strip_6437404_composite.tif\n",
            "Date: 2023-04-15\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6457722_composite.tif\n",
            "Date: 2023-04-24\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-04-24_strip_6458341_composite.tif\n",
            "Date: 2023-04-24\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-05-17_strip_6515202_composite.tif\n",
            "Date: 2023-05-17\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-02_strip_6551855_composite.tif\n",
            "Date: 2023-06-02\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-06-30_strip_6615371_composite.tif\n",
            "Date: 2023-06-30\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-14_strip_6647556_composite.tif\n",
            "Date: 2023-07-14\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-07-30_strip_6682704_composite.tif\n",
            "Date: 2023-07-30\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-16_strip_6717370_composite.tif\n",
            "Date: 2023-08-16\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-08-30_strip_6743722_composite.tif\n",
            "Date: 2023-08-30\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-13_strip_6770796_composite.tif\n",
            "Date: 2023-09-13\n",
            "Original filename: /content/drive/MyDrive/Colab Data/GIS Day 2023/Workshop 2/data/2023-09-27_strip_6800033_composite.tif\n",
            "Date: 2023-09-27\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Next, we can plot all dates and values as a scatterplot."
      ],
      "metadata": {
        "id": "Y2UmOMUJ6zXA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "f, ax = plt.subplots()\n",
        "\n",
        "ax.scatter(pd.to_datetime(dates), medians)\n",
        "for label in ax.get_xticklabels():\n",
        "    label.set_rotation(45)\n",
        "ax.set_ylabel(\"NDWI\");"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 469
        },
        "id": "VBNVN6916u_Q",
        "outputId": "514a1ccf-3f69-4c1c-c702-116ae3c67a41"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "[encoded data, 25956 characters]\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "From this plot, we can guess that Lake Lagunita began flooding in mid-December.\n",
        "\n",
        "Congratulations! You've finished the Processing with Python workshop."
      ],
      "metadata": {
        "id": "re2YCRnc7J82"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "84Q0sDG960kY"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}