sentinel2_bandarithmetic.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "sentinel2_bandarithmetic.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "VhFf1ATGWx5-"
      },
      "source": [
        "# SENTINEL2_BandArithmetic -----------------------------------------------------\n",
        "# end to end implementation of the following indexes\n",
        "\n",
        "# BSI: The Bare Soil Index\n",
        "# FDI: The Floating Debris Index (see: Finding Plastic Patches in Coastal Waters using Optical Satellite Data, Lauren Biermann et al)\n",
        "# NBR: The Normalized Burn Ratio\n",
        "# NDVI: The Normalized Vegegation Index\n",
        "# NDWI: The Normalized Difference Water Index\n",
        "# TCI: True Color Image (this task directly from the tci image in newer sentinell2 payloads)\n",
        "\n",
        "# February - May 2021, RTS\n",
        "#-------------------------------------------------------------------------------\n",
        "# Notes:\n",
        "# Ensure the file sentinel2_helper.py is in the correct locatation (see below)\n",
        "# Assumes you already have downloaded ++ a single sentinel2 product ++ with 'Sentinel2_getdata'\n",
        "# and that this download resides in the 'sentinel' directory (see below)\n",
        "# Assumes you are using the same geojson file here as in 'Sentinel2_getdata'\n",
        "\n",
        "# Process:\n",
        "# First band operation: run all cells from top to bottom\n",
        "# Subsequent band operations: continue from 'operations for NDWI, NDVI, NBR, BSI, FDI, TCI'\n",
        "# (you can omit the previous steps as the requisite variables persist in the jupyter notebook session)\n",
        "# Last band operation: run the final cell to clean up intermediate files and \n",
        "# archive the sentinal2 download, if desired.\n",
        "#-------------------------------------------------------------------------------"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "V6XgMNN5WiW_"
      },
      "source": [
        "import sys\n",
        "#clear all variables (just in case)\n",
        "sys.modules[__name__].__dict__.clear()\n",
        "import os, sys, time, random, itertools, shutil\n",
        "\n",
        "#mount the google drive\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive')\n",
        "\n",
        "#maximize screen real estate\n",
        "from IPython.core.display import display, HTML\n",
        "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rFMU20yHXAb2"
      },
      "source": [
        "#all the required installations here\n",
        "!pip install rasterio --upgrade\n",
        "!pip install rioxarray --upgrade\n",
        "!pip install geopandas --upgrade\n",
        "!pip install earthpy --upgrade\n",
        "\n",
        "print('\\n\\nFINISHED installing rasterio rioxarray, geopandas, earthpy....')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZXKlEMbNRR4-"
      },
      "source": [
        "# variables specific to your CoLab Jupyter notebook setup ----------------------\n",
        "\n",
        "start = '/content/drive/MyDrive/your directory/'\n",
        "\n",
        "# update your path to include loctions of helper files, specific to your setup\n",
        "codepath = start + 'code/'\n",
        "sys.path.append(codepath)\n",
        "\n",
        "#set other directories (create them manually on your drive)\n",
        "datapath = start + 'data/'\n",
        "searchpath = datapath + 'sentinel/'\n",
        "archivepath = datapath + 'archive/'\n",
        "\n",
        "# place a valid geojson file in the .../data directory\n",
        "map = 'bali_kuta_beaches.geojson'\n",
        "geojsonpath = datapath"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0_5Cy6p5avJ3"
      },
      "source": [
        "# all imports here\n",
        "import rasterio\n",
        "import rioxarray\n",
        "from rasterio.mask import mask\n",
        "from rasterio.plot import show\n",
        "from skimage.io import imsave\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy\n",
        "import geopandas\n",
        "import earthpy.plot as ep\n",
        "import warnings\n",
        "warnings.filterwarnings('ignore')\n",
        "\n",
        "#  homebrew\n",
        "from sentinel2_helper import *"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RWaRn77iab9F"
      },
      "source": [
        "# <start here for the first pass or if you have remove itermediate files with the last module on this page>\n",
        "# these next steps assume you have already downloaded zipped satellite image package and that is resides in your searchpath\n",
        "unpack(searchpath)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "u95XduFUajDC"
      },
      "source": [
        "imageroot = get_imageroot(searchpath)\n",
        "print('Path to band data: ', imageroot)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TZYb3KLIbNY_"
      },
      "source": [
        "# <continue here if you already processed an operation and have not yet archived the downloaded sentinel package with the last module on this page>\n",
        "# prepare and perform operations for NDWI, NDVI, NBR, BSI, FDI, TCI\n",
        "\n",
        "band_operation = 'TCI'\n",
        "#-------------------------------------------------------------------------------\n",
        "prefix = band_operation.lower()\n",
        "img_type = prefix + '.tif'\n",
        "img_ccr_type = prefix + '_4326.tif'\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "if (band_operation == 'NDVI'):\n",
        "  band4=rasterio.open(imageroot + 'B04.jp2')      #red\n",
        "  band8=rasterio.open(imageroot + 'B08.jp2')      #nir\n",
        "  red = band4.read()\n",
        "  nir = band8.read()\n",
        "\n",
        "  refband = band4\n",
        "  result_img = ndvi(red, nir)\n",
        "  convert_jp2_tif(searchpath, img_type, refband, result_img)\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "elif (band_operation == 'NDWI'):\n",
        "  band3=rasterio.open(imageroot + 'B03.jp2')      #green\n",
        "  band8=rasterio.open(imageroot + 'B08.jp2')      #nir\n",
        "  green = band3.read()\n",
        "  nir = band8.read()\n",
        "\n",
        "  refband = band3\n",
        "  result_img  = ndwi(green, nir) \n",
        "  convert_jp2_tif(searchpath, img_type, refband, result_img)\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "elif (band_operation == 'NBR'):\n",
        "  band12=rasterio.open(imageroot + 'B12.jp2' )    #SWIR 20m\n",
        "  band8A=rasterio.open(imageroot + 'B8A.jp2')     #nir 20m\n",
        "  swir = band12.read()\n",
        "  nir = band8A.read()\n",
        "\n",
        "  refband = band8A\n",
        "  result_img  = nbr(nir, swir) \n",
        "  convert_jp2_tif(searchpath, img_type, refband, result_img)\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "elif (band_operation == 'BSI'):\n",
        "  b11 = 'B11.jp2'                                 #swir 1610nm, 20m\n",
        "  scale_factor = 2                                #upsample from 20m res to 10m\n",
        "\n",
        "  resampled_b11 = resample(searchpath, imageroot, b11, scale_factor)\n",
        "\n",
        "  band2=rasterio.open(imageroot + 'B02.jp2')      #blue, 10m\n",
        "  band8=rasterio.open(imageroot + 'B08.jp2')      #nir, 10m\n",
        "  band4=rasterio.open(imageroot + 'B04.jp2')      #red, 10m\n",
        "  sw = rasterio.open(searchpath + resampled_b11)\n",
        "\n",
        "  blue = band2.read()\n",
        "  nir = band8.read()\n",
        "  red = band4.read()\n",
        "  swir = sw.read()\n",
        "  \n",
        "  refband = band2\n",
        "  result_img  = bsi(red, blue, nir, swir) \n",
        "  convert_jp2_tif(searchpath, img_type, refband, result_img)\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "elif (band_operation == 'FDI'):\n",
        "  b11 = 'B11.jp2'                                 #swir 1610nm, 20m\n",
        "  b6 = 'B06.jp2'                                  #rededge2 740nm\n",
        "  scale_factor = 2                                #upsample from 20m res to 10m\n",
        "\n",
        "  resampled_b11 = resample(searchpath, imageroot, b11, scale_factor)\n",
        "  resampled_b6 = resample(searchpath, imageroot, b6, scale_factor)\n",
        "  band8 = rasterio.open(imageroot + 'B08.jp2')    #nir\n",
        "  nir = band8.read()\n",
        "  sw = rasterio.open(searchpath + resampled_b11)\n",
        "  swir = sw.read()\n",
        "  red = rasterio.open(searchpath + resampled_b6)\n",
        "  rededge = red.read()\n",
        "\n",
        "  lambda_nir = 832.8                              #nm\n",
        "  lambda_red = 664.6                              #nm\n",
        "  lambda_swir = 1613.7                            #nm\n",
        "\n",
        "  refband = band8\n",
        "  result_img = fdi(nir, rededge, swir, lambda_nir, lambda_red, lambda_swir)\n",
        "  convert_jp2_tif(searchpath, img_type, refband, result_img)\n",
        "\n",
        "#-------------------------------------------------------------------------------\n",
        "elif (band_operation == 'TCI'):\n",
        "  tci = get_tci(searchpath)\n",
        "  convert_tci_tif(tci, searchpath, img_type)\n",
        "#-------------------------------------------------------------------------------\n",
        "else:\n",
        "  print(' ... pick from one of the available operations: NDWI, NDVI, NBR, BSI, FDI, TCI')\n",
        "\n",
        "print(band_operation + ' tasks have completed... continue to next cells')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UKhRTGzQ4NO5"
      },
      "source": [
        "#convert coordinate system\n",
        "target_epsg = \"EPSG:4326\"\n",
        "convert_ccr(searchpath, img_type, target_epsg, img_ccr_type)\n",
        "print('CRS conversion complete')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HHgR8I3y5p-S"
      },
      "source": [
        "#check the geojson file (extent of the area we are interested in)\n",
        "gdf = geopandas.read_file(geojsonpath + map)\n",
        "gdf.plot();"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IABZXI9v558A"
      },
      "source": [
        "#apply mask, rename image, display and save\n",
        "current_4326_img = rasterio.open(searchpath + img_ccr_type) \n",
        "#apply the geojson file produced mask\n",
        "masked, mask_transform = mask(dataset=current_4326_img, shapes=gdf.geometry, crop=True)\n",
        "newname = create_name(prefix, imageroot, geojsonpath, map)\n",
        "\n",
        "#display and save\n",
        "if ('tci' in img_ccr_type):\n",
        "  show(masked, transform=mask_transform);\n",
        "  temp = rasterio.plot.reshape_as_image(masked)\n",
        "  imsave(datapath + newname +'.jpeg', temp)\n",
        "\n",
        "else:\n",
        "  colormap =   'gray'                                                               #'gray'    #'RdYlGn'  #'viridis' ... others\n",
        "  ax = ep.plot_bands(masked, cmap=colormap, cols=1, vmin=-1, vmax=1, cbar=False)    #cbar=False OR True\n",
        "  ax.figure.savefig(datapath + newname + '.jpeg')\n",
        "  ax.figure.show()\n",
        "\n",
        "print(newname + ' saved to searchpath')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qbLMYT24-rhq"
      },
      "source": [
        "#clean up\n",
        "REUSE = True\n",
        "ARCHIVE = True\n",
        "archive(searchpath, archivepath, REUSE, ARCHIVE)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}