6_0_Remote Sensing Image analysis in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading satellite raster files with Rasterio\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the required module\n",
    "import os\n",
    "import sys\n",
    "import rasterio\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sys.path.append(os.path.abspath('..')) # set default directory where this notebook is running from"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rasterio is a very useful module for raster processing for reading and writing several different raster formats in Python. \n",
    "\n",
    "Rasterio is based on GDAL and Python automatically registers all GDAL drivers. \n",
    "Commonly used file formats include e.g., TIFF and GeoTIFF, ASCII Grid and Erdas Imagine .img -files.\n",
    "\n",
    "Landsat 8 bands are stored as separate GeoTIFF -files in the original package. Each band contains information of surface reflectance from different ranges of the electromagnetic spectrum.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "installation of GDAL at https://anaconda.org/conda-forge/gdal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read raster file\n",
    "fp= \"Data/CPH_masked.tif\"\n",
    "\n",
    "# Open the file:\n",
    "raster = rasterio.open(fp)\n",
    "\n",
    "# Check type of the variable 'raster'\n",
    "print (\"The object type of the raster image is: \", type(raster))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that our raster variable is a rasterio._io.RasterReader type, which means that we have opened the raster file ready for reading.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading raster file properties\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore the metadata of the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Projection\n",
    "raster.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check what projection it contains\n",
    "https://epsg.io/?q=32632 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dimensions --> width\n",
    "raster.width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dimensions --> height\n",
    "raster.height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of bands\n",
    "raster.count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bounds of the file\n",
    "raster.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Driver (data format)\n",
    "raster.driver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# No data values for all channels\n",
    "raster.nodatavals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All Metadata for the whole raster dataset\n",
    "raster.meta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(raster.meta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write a code to print each part of metadata in a single line\n",
    "meta_dic= raster.meta\n",
    "\n",
    "for k, v in meta_dic.items():\n",
    "    print (\"The \", k, \"is \", v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get raster bands\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different bands of a satellite images are often stacked together in one raster dataset. In our case, all 11 bands of the Landsat 8 scene are included in our GeoTIFF and the count is hence 11.\n",
    "Remember landsat 8 bands (https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites?qt-news_science_products=0#qt-news_science_products)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read the raster band as separate variable\n",
    "band2 = raster.read(2)\n",
    "band3 = raster.read(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check type of the variable 'band'\n",
    "type(band2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data type of the values\n",
    "band2.dtype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "have a look at raster data types at https://scikit-image.org/docs/dev/user_guide/data_types.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Band statistics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOW, let’s check the values that are stored in the band. \n",
    "The values of the bands are stored as numpy arrays, so it is extremely easy to calculate basic statistics by using basic numpy functions.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read all bands\n",
    "array = raster.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate statistics for each band\n",
    "stats = []\n",
    "for band in array:\n",
    "    stats.append({\n",
    "        'min': band.min(),\n",
    "        'mean': band.mean(),\n",
    "        'median': np.median(band),\n",
    "        'max': band.max()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing raster layers\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, it is always highly useful to take a look how the data looks like. This is easy with the \"rasterio.plot.show()\" function of rasterio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import rasterio\n",
    "from rasterio.plot import show\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "fp= \"Data/CPH_masked.tif\"\n",
    "\n",
    "# Open the file:\n",
    "raster = rasterio.open(fp)\n",
    "\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "show((raster, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image looks dark and we can not see clear content, so we should do image enhancement in order to strecth the diginal numbers (DN) mkaing a good contrast between grey values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the grid values into numpy arrays\n",
    "red = raster.read(4)\n",
    "green = raster.read(3)\n",
    "blue = raster.read(2)\n",
    "\n",
    "# Function to normalize the grid values\n",
    "def normalize(array):\n",
    "    \"\"\"Normalizes numpy arrays into scale 0.0 - 1.0\"\"\"\n",
    "    array_min, array_max = array.min(), array.max()\n",
    "    return ((array - array_min)/(array_max - array_min))\n",
    "\n",
    "# Normalize the bands\n",
    "redn = normalize(red)\n",
    "greenn = normalize(green)\n",
    "bluen = normalize(blue)\n",
    "\n",
    "print(\"Normalized bands\")\n",
    "print(redn.min(), '-', redn.max(), 'mean:', redn.mean())\n",
    "print(greenn.min(), '-', greenn.max(), 'mean:', greenn.mean())\n",
    "print(bluen.min(), '-', bluen.max(), 'mean:', bluen.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as you can see, the grey values are very close to the first quantile of values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create RGB natural color composite\n",
    "rgb = np.dstack((redn, greenn, bluen))\n",
    "\n",
    "# Let's see how our color composite looks like\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.imshow(rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the grid values into numpy arrays\n",
    "swir2 = raster.read(7)\n",
    "swir = raster.read(6)\n",
    "nir = raster.read(5)\n",
    "red = raster.read(4)\n",
    "green = raster.read(3)\n",
    "\n",
    "# Normalize the values using the function that we defined earlier\n",
    "nirn = normalize(nir)\n",
    "redn = normalize(red)\n",
    "greenn = normalize(green)\n",
    "\n",
    "# Create the composite by stacking\n",
    "nrg = np.dstack((nirn, redn, greenn))\n",
    "\n",
    "# Let's see how our color composite looks like\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.imshow(nrg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GDAL\n",
    "\n",
    "GDAL is a translator library for raster and vector geospatial data formats that is released under an MIT style Open Source License by the Open Source Geospatial Foundation. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing.\n",
    "\n",
    "[MORE](https://gdal.org)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's repeat some of these steps with GDAL\n",
    "from osgeo import gdal\n",
    "image_data = gdal.Open(\"Data/CPH_masked.tif\")\n",
    "\n",
    "\n",
    "total_number_of_bands = image_data.RasterCount\n",
    "print (\"total number of bands in this image is %s \" %total_number_of_bands)\n",
    "rows = image_data.RasterYSize\n",
    "columns = image_data.RasterXSize\n",
    "print (\"rows: %s pixels\" %columns)\n",
    "print (\"columns: %s pixels\"  %rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projection = image_data.GetProjection()\n",
    "print (projection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's read some of the bands\n",
    "band_1 = image_data.GetRasterBand(1).ReadAsArray()\n",
    "band_2 = image_data.GetRasterBand(2).ReadAsArray()\n",
    "band_3 = image_data.GetRasterBand(3).ReadAsArray()\n",
    "band_4 = image_data.GetRasterBand(4).ReadAsArray()\n",
    "band_5 = image_data.GetRasterBand(5).ReadAsArray()\n",
    "band_6 = image_data.GetRasterBand(6).ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_color = np.dstack([band_3, band_2, band_1])\n",
    "color_infrared = np.dstack([band_4, band_3, band_2])\n",
    "false_natural_color = np.dstack([band_5, band_4, band_3])\n",
    "false_color_II = np.dstack([band_6, band_5, band_3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show((raster, 1), cmap='Reds', title ='landsat image')\n",
    "show((raster, 2), cmap='Greens')\n",
    "show((raster, 3), cmap='Blues')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Histogram of the raster data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is recommended to always look at the histogram of your data. It is as easy as using \"rasterio.plot.show_hist()\" function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.plot import show_hist\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "\n",
    "show_hist(raster, bins=50, lw=0.0, stacked=False, alpha=0.3,\n",
    "          histtype='stepfilled' ,title=\"Histogram\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clipping or Masking a raster\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One common task in raster processing is clipping raster files based on a Polygon. We are going to clip the raster data from a polygon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import required packages\n",
    "from rasterio.plot import show\n",
    "from rasterio.plot import show_hist\n",
    "from rasterio.mask import mask\n",
    "from shapely.geometry import box\n",
    "import geopandas as gpd\n",
    "from fiona.crs import from_epsg\n",
    "import pycrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to specify input and output filepaths\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the output path or filename\n",
    "out_tif = r\"CPH_masked_CLIP.tif\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "WHY using \"r\" in the line above? (Since some of you faced this problem)\n",
    "\n",
    "In Python, backslash is used to signify special characters.\n",
    "\n",
    "e.g. \"hello\\nworld\" -- the \"\\n\" means a newline. Remember we used it for printing in a new line.\n",
    "\n",
    "Path names on Windows tend to have backslashes in them. But we want them to mean actual backslashes, not special characters.\n",
    "\n",
    "r stands for \"raw\" and will cause backslashes in the string to be interpreted as actual backslashes rather than special characters.\n",
    "\n",
    "e.g. r\"hello\\nworld\" literally means the characters \"hello\\nworld\". Again, try printing it.\n",
    "\n",
    "More info is in the Python docs, it's a good idea to search them for questions like these.\n",
    "\n",
    "https://docs.python.org/3/tutorial/introduction.html#strings\n",
    "\n",
    "Takeaway: ALWAYS put an \"r\" before a path/directory"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the raster in read mode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = rasterio.open(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the data\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "show((data, 2), cmap='terrain')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOW, let's create a bounding box with Shapely for clipping\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# UTM coordinates EPSG: 32632\n",
    "# http://bboxfinder.com/\n",
    "from shapely.geometry import box\n",
    "\n",
    "bbox = box(717621.8079, 6172470.4282, 731060.0044, 6180333.8423)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "# Insert the bbox into a GeoDataFrame\n",
    "geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=\"EPSG:32632\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to get the coordinates of the geometry in such a format that rasterio wants them. This can be conducted easily with following function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getFeatures(gdf):\n",
    "    \"\"\"Function to parse features from GeoDataFrame in such a manner that rasterio wants them\"\"\"\n",
    "    import json\n",
    "    return [json.loads(gdf.to_json())['features'][0]['geometry']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the geometry coordinates by using the function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = getFeatures(geo)\n",
    "\n",
    "print(coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so rasterio wants to have the coordinates of the Polygon in this kind of format as above!\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to clip the raster with the polygon using the coords variable that we just created. Clipping the raster can be done easily with the mask function that we imported in the beginning from rasterio, and specifying clip=True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting the output file and its ttransformation \n",
    "out_img , out_transform= mask(dataset=data, shapes=coords, crop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to modify the metadata. Let’s start by copying the metadata from the original data file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy the metadata\n",
    "out_meta = data.meta.copy()\n",
    "\n",
    "print(out_meta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to parse the EPSG value from the CRS so that we can create a Proj4 string using PyCRS library (to ensure that the projection information is saved correctly).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.crs.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parse EPSG code\n",
    "epsg_code = int(data.crs.data['init'][5:])\n",
    "print(epsg_code)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to update the metadata with new dimensions, transform (affine) and CRS (as Proj4 text)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_img.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_meta.update({\"driver\": \"GTiff\",\n",
    "                 \"height\": out_img.shape[1],\n",
    "                 \"width\": out_img.shape[2],\n",
    "                 \"transform\": out_transform,\n",
    "                 \"crs\": epsg_code})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can save the clipped raster to disk with following command.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(out_tif, \"w\", **out_meta) as dest:\n",
    "    dest.write(out_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s still check that the result is correct by plotting our new clipped raster.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clipped = rasterio.open(out_tif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 10))\n",
    "show((clipped, 5), cmap='terrain')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clipped.meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raster calculations\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculations between bands of a raster file is a common GIS task. \n",
    "We will calculate NDVI (Normalized difference vegetation index) of Copenhagen based on our Landsat dataset. Conducting calculations with rasterio is fairly straightforward if the extent etc. matches because the values of the rasters are stored as numpy arrays (similar to the columns stored in Geo/Pandas, i.e. Series).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let’s start by importing the necessary modules rasterio and numpy\n",
    "from rasterio.plot import show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let’s read the masked file.\n",
    "fp = r\"CPH_masked_CLIP.tif\"\n",
    "\n",
    "raster = rasterio.open(fp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For calculating the NDVI (Normalized difference vegetation index) you need two bands:::\n",
    "band-4 which is the Red channel and band-5 which is the Near Infrared (NIR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let’s read those bands from our raster source\n",
    "red = raster.read(4)\n",
    "nir = raster.read(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "imgplot = plt.imshow(nir, norm=LogNorm(),  cmap='summer')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "imgplot = plt.imshow(red, norm=LogNorm(), cmap='summer')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we are ready to calculate the NDVI. \n",
    "# First, we can create a filter where we calculate the values on such pixels that have a value larger than 0.\n",
    "\n",
    "check = np.logical_or ( red > 0, nir > 0 )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now we can apply that filter and calculate the ndvi index.\n",
    "ndvi = np.where ( check,  (nir - red ) / ( nir + red ), -999 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's calcuate mean of ndvi\n",
    "ndvi.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's calcuate SD of ndvi\n",
    "ndvi.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "you can look for colormaps at https://matplotlib.org/3.1.1/tutorials/colors/colormaps.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_pure = np.logical_and ( -1 <= ndvi, ndvi <= 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "imgplot = plt.imshow(ndvi, norm=LogNorm(), cmap='PiYG')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density slicing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap, BoundaryNorm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Create a colormap from a list of colors\n",
    "colors = ['white','blue','grey','green']\n",
    "\n",
    "cmap = ListedColormap(colors)\n",
    "\n",
    "class_bins = [-1, -0.5, 0.02, 0.3, 0.8]\n",
    "norm = BoundaryNorm(class_bins,len(colors))\n",
    "\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "imgplot = plt.imshow(ndvi, norm=norm,  cmap=cmap)\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear stretch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "an alternative way of linear streching using Scikit-image library "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.figure(figsize=(20, 10))\n",
    "\n",
    "plt.imshow(rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import exposure\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "rgb = np.dstack((redn, greenn, bluen))\n",
    "\n",
    "def linearStretch(input, percent):\n",
    "    pLow, pHigh = np.percentile(input[~np.isnan(input)], (percent, 100 - percent))\n",
    "    img_rescale = exposure.rescale_intensity(input, in_range=(pLow, pHigh))\n",
    "    return img_rescale\n",
    "\n",
    "img_rescaled = linearStretch(rgb, 1.15)\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.imshow(img_rescaled)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a histogram equalization solution that is working better than linear stretch.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Some of my images have NaNs for NoData so I have to remove them when creating the histogram\n",
    "shape = rgb.shape\n",
    "rgb_vector = rgb.reshape([rgb.shape[0] * rgb.shape[1], rgb.shape[2]])\n",
    "rgb_vector = rgb_vector[~np.isnan(rgb_vector).any(axis=1)]\n",
    "\n",
    "# View histogram of RGB values\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "fig.set_facecolor('white')\n",
    "for color, channel in zip('rgb_vector', np.rollaxis(rgb, axis=-1)):\n",
    "    counts, centers = exposure.histogram(channel)\n",
    "    plt.plot(centers[1::], counts[1::], color=color)\n",
    "plt.show()\n",
    "\n",
    "# Get cutoff values based on standard deviations. Ideally these would be on either side of each histogram peak and cutoff the tail. \n",
    "lims = []\n",
    "for i in range(3):\n",
    "    x = np.mean(rgb_vector[:, i])\n",
    "    sd = np.std(rgb_vector[:, i])\n",
    "    low = x-(1.75*sd)  # Adjust the coefficient here if the image doesn't look right\n",
    "    high = x + (1.75 * sd)  # Adjust the coefficient here if the image doesn't look right\n",
    "    if low < 0:\n",
    "        low = 0\n",
    "    if high > 1:\n",
    "        high = 1\n",
    "    lims.append((low, high))\n",
    "\n",
    "r = exposure.rescale_intensity(rgb[:, :, 0], in_range=lims[0])\n",
    "g = exposure.rescale_intensity(rgb[:, :, 1], in_range=lims[1])\n",
    "b = exposure.rescale_intensity(rgb[:, :, 2], in_range=lims[2])\n",
    "rgb_enhanced = np.dstack((r, g, b))\n",
    "plt.figure(figsize=(20, 10))\n",
    "plt.imshow(rgb_enhanced)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}