4-thresholding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thresholding\n",
    "\n",
    "In this notebook we will focus in change detection by the using of thresholds.\n",
    "We will see different methods to define the thresholds and will apply them to some datasets previously introduced"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import chi2,norm,normaltest\n",
    "from skimage.filters import threshold_otsu, threshold_yen, threshold_li\n",
    "from helpers import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to detect changes by using thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize temporal serie and see differences\n",
    "img_data1=load_image('data/backscatterNiamey20180130.tif')\n",
    "img_data2=load_image('data/backscatterNiamey20180611.tif')\n",
    "img_data3=load_image('data/backscatterNiamey20180927.tif')\n",
    "layers = []\n",
    "layers.append(img_data1)\n",
    "layers.append(img_data2)\n",
    "layers.append(img_data3)\n",
    "image = np.dstack(layers)\n",
    "plot_bands_row(image, colormap=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute of difference image\n",
    "Backscatter1 - Backscatter2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#How to visualize changes\n",
    "diff=img_data2[:,:,0].astype(np.int32)-img_data1[:,:,0].astype(np.int32)\n",
    "diff2=img_data3[:,:,0].astype(np.int32)-img_data2[:,:,0].astype(np.int32)\n",
    "layersDiff = []\n",
    "layersDiff.append(diff)\n",
    "layersDiff.append(diff2)\n",
    "imageDiff = np.dstack(layersDiff)\n",
    "plot_bands_row(imageDiff,vmin=-255,vmax=255,cmap=plt.cm.bwr, colormap=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#RGB generation\n",
    "plot_rgb(img_data1,img_data2,img_data2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thresholds \n",
    "The pixels whose differences is lower than a certain min_threshold or higher than the upper_threshold can be considered change.\n",
    "\n",
    "How can this threshold be computed?\n",
    "Is the same threshold valid for all the areas?\n",
    "Are the thresholds the same for a time serie of images?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User defined threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Positive changes\n",
    "#maskedDiff50 = ma.masked_where(diff>50, diff)\n",
    "maskedDiff100 = ma.masked_where(diff>100, diff)\n",
    "maskedDiff150 = ma.masked_where(diff>150, diff)\n",
    "maskedDiff200 = ma.masked_where(diff>200, diff)\n",
    "maskedDiff250 = ma.masked_where(diff>250, diff)\n",
    "layers = []\n",
    "#layers.append(maskedDiff50.mask)\n",
    "layers.append(maskedDiff100.mask)\n",
    "layers.append(maskedDiff150.mask)\n",
    "layers.append(maskedDiff200.mask)\n",
    "layers.append(maskedDiff250.mask)\n",
    "image = np.dstack(layers)\n",
    "plot_bands_row(image,vmin=0,vmax=1)\n",
    "\n",
    "#Negative changes\n",
    "#maskedDiff50 = ma.masked_where(diff<-50, diff)\n",
    "maskedDiff100 = ma.masked_where(diff<-100, diff)\n",
    "maskedDiff150 = ma.masked_where(diff<-150, diff)\n",
    "maskedDiff200 = ma.masked_where(diff<-200, diff)\n",
    "maskedDiff250 = ma.masked_where(diff<-250, diff)\n",
    "layers = []\n",
    "#layers.append(maskedDiff50.mask)\n",
    "layers.append(maskedDiff100.mask)\n",
    "layers.append(maskedDiff150.mask)\n",
    "layers.append(maskedDiff200.mask)\n",
    "layers.append(maskedDiff250.mask)\n",
    "image = np.dstack(layers)\n",
    "plot_bands_row(image,vmin=0,vmax=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Another way to visualize changes with both min and max thresholds\n",
    "mask = np.zeros(diff.shape).astype(np.uint8)\n",
    "mask = np.logical_or(np.logical_or(diff>200,diff<(0-200)),mask)\n",
    "plot_bands_row(mask,vmin=0,vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: find some thresholds to be able to distinguish:\n",
    "* New constructions in Charpiya\n",
    "* Flooding areas in Charpiya\n",
    "* Burned Areas in Amazonas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Automatic threshold computation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lets have a look to the histogram of the diff images in Niamey\n",
    "\n",
    "plt.hist(diff.flatten()) \n",
    "plt.title(\"Diff histogram\") \n",
    "plt.show()\n",
    "\n",
    "plt.hist(diff2.flatten()) \n",
    "plt.title(\"Diff2 histogram\") \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can consider that the diff value, if no changes, would be similar to a gaussian distribution. Ideally it would have mean = 0, but it depends on:\n",
    "* acquisition conditions\n",
    "* preprocessing of images (e.g. data conversion...)\n",
    "\n",
    "With this approach we can compute with the level of statistical significance that we want, the values that are 'anomalies'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# statistical threshold (normal, chi2...)\n",
    "#Normalize diff\n",
    "std=np.std(diff)\n",
    "mean=np.mean(diff)\n",
    "diffNorm=(diff-mean)/std\n",
    "print('Mean diff: ' +str(mean))\n",
    "print('Std diff: ' +str(std))\n",
    "#print(diffNorm)\n",
    "\n",
    "error1=norm.isf(0.000001)\n",
    "print('Normalized Error: ' +str(error1))\n",
    "print('Error: ' +str(error1*std+mean))\n",
    "\n",
    "mask = np.zeros(diff.shape).astype(np.uint8)\n",
    "mask = np.logical_or(np.logical_or(diffNorm>error1,diffNorm<(0-error1)),mask)\n",
    "plot_bands_row(mask,vmin=0,vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposed exercise\n",
    "* Apply same method to difference at date 2 in Niamey\n",
    " * Is it OK to use the same threshold?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#other standards thresholds\n",
    "yen=threshold_yen(np.absolute(diff))\n",
    "print(yen)\n",
    "mask = np.zeros(diff.shape).astype(np.uint8)\n",
    "mask = np.logical_or(np.absolute(diff)>yen,mask)\n",
    "plot_bands_row(mask,vmin=0,vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposed exercise \n",
    "* Apply and visualize results using other thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary of methodology\n",
    "\n",
    "Till now, we have computed the difference image and then we have selected a threshold to apply. We have seen that are different methods for the threshold calculation but what about the difference image?\n",
    "\n",
    "There are different 'difference' images that are used:\n",
    "- Backscatter1 - Backscatter2\n",
    "- Ratio  of backscatter\n",
    "- Ratio of log of backscatters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "charp_data1=load_image('data/backscatterCharpiya20171220.tif')\n",
    "charp_data2=load_image('data/backscatterCharpiya20180910.tif')\n",
    "dlr=np.log(charp_data1/charp_data2)\n",
    "layers = []\n",
    "layers.append(charp_data1)\n",
    "layers.append(charp_data2)\n",
    "#layers.append(dlr)\n",
    "image = np.dstack(layers)\n",
    "plot_bands_row(image,vmin=0,vmax=1)\n",
    "plot_bands_row(dlr,vmin=-4,vmax=4,cmap=plt.cm.bwr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proposed exercise \n",
    "* Compute the binary map of new constructions/floodings in Charpiya using the time series using one or more of the methods proposed\n",
    "* Compute the binary map of Burned Area using the time series using one or more of the methods proposed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "better",
   "language": "python",
   "name": "better"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}