examples/Advanced_Band_Math_Tutorial_Using_SWIR_Imagery.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial covers using SWIR imagery for creating spectral indices. SWIR stands for **S**hort **W**avelength **I**nfra**R**ed and represents an additional swath of spectrum that can be used for analyzing satellite imagery. The use of the SWIR band is not necessarily more complicated than other bands when creating spectral indices, but accessing SWIR imagery can be.  Instead of using the imagery tab to the left of the notebook, you'll need to directly query the GBDX API to find and order SWIR imagery.  \n",
    "\n",
    "Specifically we'll need to look for Worldview-3 (WV03_SWIR) images since it is the DigitalGlobe satellite which has a SWIR sensor to generate the imagery we need for our analysis.  In addition to identifying the appopriate sensor and satellite we'll also need to find an AOI.  To create an AOI we can use the imagery tab to generate a bounding box that can be cut and pasted for our code below. In addition to our SWIR search we will also look for VNIR (**V**ery **N**ear **I**nfra**R**ed) since several spectral indices require bands from both image products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import json\n",
    "\n",
    "try:\n",
    "    from shapely.geometry import box\n",
    "except:\n",
    "    print(\"Installing fiona...\")\n",
    "    !pip install shapely -qq\n",
    "    from shapely.geometry import box\n",
    "\n",
    "from gbdxtools import Interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gbdx = Interface()\n",
    "\n",
    "bbox = map(float, '-73.9651, 40.7897, -73.9458, 40.8025'.split(\",\"))\n",
    "\n",
    "aoi = box(*bbox).wkt\n",
    "query = \"(item_type:WV03_VNIR OR item_type:WV03_SWIR)\"\n",
    "query += \" AND NOT item_type:IDAHOImage AND item_type:DigitalGlobeProduct\"   \n",
    "results = gbdx.vectors.query(aoi, query)\n",
    "    \n",
    "for r in results:\n",
    "    props = r['properties']\n",
    "    print(props['attributes']['catalogID'], props['item_date'][0:10], props['item_type'][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gbdxtools import CatalogImage\n",
    "\n",
    "vnir_id = '104001001DB7BA00'\n",
    "swir_id = '104A01001DB7BA00'\n",
    "\n",
    "vnir = CatalogImage(vnir_id, proj=\"EPSG:3857\", bbox=[-73.9651, 40.7897, -73.9458, 40.8025])\n",
    "swir = CatalogImage(swir_id, proj=\"EPSG:3857\", bbox=[-73.9651, 40.7897, -73.9458, 40.8025])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we will plot our data so we can see the difference between the VNIR data and the SWIR data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from shapely.geometry import box\n",
    "import numpy as np\n",
    "\n",
    "plt.figure(figsize=(15,5))\n",
    "\n",
    "titles = ['Visbile Near-Infrared', 'Shortwave Infrared']\n",
    "\n",
    "bbox = box(-73.9651, 40.7897, -73.9458, 40.8025).bounds\n",
    "\n",
    "for idx, img in enumerate([vnir.aoi(bbox=bbox).rgb(), swir.aoi(bbox=bbox)[7,:,:]]):\n",
    "    plt.subplot(1,2,idx+1)\n",
    "    plt.axis('off')\n",
    "    plt.imshow(img, cmap=\"Greys_r\")\n",
    "    plt.title(titles[idx])\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the VNIR data is part of the visibile spectrum imagery product, which is in color and at a higher resolution that the black and white SWIR image.  The SWIR band in the imagery is not part of the visible spectrum of light - it falls between about 1400 and 3000 nanometers in wavelength.  This is why it doesn't show up in color and it accentuates shadow and contrast in the image.  This is also what makes SWIR particularily good for feature identification.  The second visible difference is the lower resolution (fuzziness) of the SWIR image.  This is because the SWIR imagery is captured with a different sensors on Worldview-3 that operate with a resolution of 3.7/7.5 meters.  The sensors that capture VNIR operate up to 1.2 meter in resolution.  \n",
    "\n",
    "For our spectral index example we'd like to calcualte a Normalized Difference Built-Up Index (NDBI).\n",
    "This NDBI index highlights urban areas where there is a higher reflectance in the shortwave-infrared (SWIR) region, compared to the near-infrared (NIR) region. Since our proposed index uses data from the two different imagery products we've highlighted; it presents a data problem for our analysis.  The VNIR and SWIR data are at different resolutions, so we need to normalize them before we can perform band math to create our spectral index.  To begin our normalization process we need to see what the pixel size difference is between our two images.  We can use the Python .shape command to get the relative sizes of our VNIR and SWIR images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "swir.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vnir.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from our .shape command that both images have eight bands, but the VNIR image has quite a few more pixels than the SWIR image. In order to create our NDBI index we are going to transform our VNIR image to be the same size as our SWIR image.  To do so we first need to decide which bands we want to use for our index and declare them.  Remember that for DigitalGlobe's visible spectrum imagery product the NIR band is six, so we'll declare that.  For the the SWIR image we'll use three."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "swir_1 = swir[3,:,:]\n",
    "nir = vnir[6,:,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've declared our bands we'll import the excellent Scikit Image library to resize our image by downsampling the pixel resolution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.transform import resize\n",
    "\n",
    "nir_rs = resize(nir.astype(np.float), swir_1.shape, order=1, mode='constant', clip=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check our image resize effort by using the .shape command again for swir_1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "swir_1.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Huzzah! Our SWIR image and our VNIR image now both have the same pixel dimensions and we can write our equation for calculating NDBI.  We can see below that our NDBI index follows a similar struture to NDVI and NHFD, where we are creating a normalized difference between two bands.  Once we have the equation in place we'll plot our results and see the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbdi = (swir_1 - nir_rs)/(swir_1 + nir_rs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def plot(image, w=15, h=12, cmap='jet'):\n",
    "    f, ax = plt.subplots(1, figsize=(w,h))\n",
    "    ax.axis('off')\n",
    "    plt.imshow(image, cmap=cmap)\n",
    "    plt.colorbar()\n",
    "    plt.show()\n",
    "\n",
    "plot(nbdi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NDBI index lives up to its billing and cleanly delineates central park and the heavily built up construction that surrounds it.  Further, the index clearly delineates the two water bodies in the image \"Harlem Meer\" and \"The Pool\".  SWIR data also is brilliant for detecting minerals, snow/ice, man made material and fires.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hub": {
   "id": "5a4d1b8ce7f21410429a150e",
   "published": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}