2. Stratified Sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stratified sampling, filter and binarize\n",
    "The following notebook does a few extra things. We use a global permanent water dataset and create points for stratified sampling over permanent water and land. The notebook then introduces Normalized Difference Water Index (MNDWI) and the stratified sampling value to get a estimated threshold to MNDWI. This is translated as a binary land water image and added back. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bqplot\n",
    "import datetime\n",
    "import dateutil.parser\n",
    "import ee\n",
    "import ipywidgets\n",
    "import ipyleaflet\n",
    "import IPython.display\n",
    "import numpy as np\n",
    "import pprint\n",
    "import pandas as pd\n",
    "import traitlets\n",
    "\n",
    "# Configure the pretty printing output & initialize earthengine.\n",
    "pp = pprint.PrettyPrinter(depth=4)\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to get sizes in Human readable format\n",
    "suffixes = ['B', 'KB', 'MB', 'GB', 'TB', 'PB']\n",
    "def humansize(nbytes):\n",
    "    i = 0\n",
    "    while nbytes >= 1024 and i < len(suffixes)-1:\n",
    "        nbytes /= 1024.\n",
    "        i += 1\n",
    "    f = ('%.2f' % nbytes).rstrip('0').rstrip('.')\n",
    "    return '%s %s' % (f, suffixes[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*While single images are great to do quick analytics, the true power of the Earth Engine environment comes with the possibility of looking at really large and heavy image collections. In GEE environment image collections have their own characteristic setup and are often composed of multiple single images. They often have the similar band structure and generally share a similar metadata structure for filtering and querying.*\n",
    "\n",
    "In this step we query an image collection and get the size and numbe of items in a collection. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of JRC Yearly Water History assets with filters: 1\n"
     ]
    }
   ],
   "source": [
    "#Public Image Collections\n",
    "l8sr = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').filterDate('2018-08-01','2018-12-01').filterMetadata('CLOUD_COVER','less_than',2)\n",
    "global_water = (ee.ImageCollection('JRC/GSW1_0/YearlyHistory')\n",
    "                .filterDate('2015-01-01','2018-01-01'))\n",
    "\n",
    "#Select a single grid to be used for the analysis\n",
    "grid=ee.FeatureCollection('users/samapriya/vector/la-subset').filterMetadata('Grids','equals',\"Grid_1\")\n",
    "                \n",
    "# Get collection size\n",
    "print('Total number of JRC Yearly Water History assets with filters: '+str(global_water.size().getInfo()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'bands': [{'crs': 'EPSG:4326',\n",
      "            'crs_transform': [0.00026949458523585647,\n",
      "                              0.0,\n",
      "                              -180.00001488697754,\n",
      "                              0.0,\n",
      "                              -0.00026949458523585647,\n",
      "                              84.00011474509029],\n",
      "            'data_type': {'max': 255,\n",
      "                          'min': 0,\n",
      "                          'precision': 'int',\n",
      "                          'type': 'PixelType'},\n",
      "            'id': 'waterClass'}],\n",
      " 'id': 'JRC/GSW1_0/YearlyHistory/31',\n",
      " 'properties': {'system:footprint': {'coordinates': [[...]],\n",
      "                                     'geodesic': False,\n",
      "                                     'type': 'Polygon'},\n",
      "                'system:index': '31',\n",
      "                'system:time_start': 1420070400000,\n",
      "                'year': 2015},\n",
      " 'type': 'Image',\n",
      " 'version': 1550150111879087}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['waterClass']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get sample image from collection\n",
    "pp.pprint(global_water.first().getInfo())\n",
    "band_names_original = global_water.first().bandNames()\n",
    "band_names_original.getInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to get tilelayer url from earthengine server\n",
    "def GetTileLayerUrl(ee_image_object):\n",
    "  map_id = ee.Image(ee_image_object).getMapId()\n",
    "  tile_url_template = \"https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}\"\n",
    "  return tile_url_template.format(**map_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "508c25d2c9e2479793e74841c37c3c38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Create a slider widget to add both Landsat 8 and PlanetScope 4B SR imagery\n",
    "map1 = ipyleaflet.Map(\n",
    "    center=(29.2963,-90.3446), zoom=12,\n",
    "    layout={'height':'500px'},\n",
    ")\n",
    "\n",
    "#Definte mask for not water\n",
    "global_nw=global_water.median().updateMask(global_water.median().eq(0))\n",
    "jrc_tile_url=GetTileLayerUrl(global_water.median().visualize(min=0, max=3, bands=['waterClass'],palette='grey,brown,FFFFFF,blue'))\n",
    "l8sr_tile_url = GetTileLayerUrl(l8sr.median().visualize(min=100, max=3500, gamma=1.5, bands= ['B5', 'B3', 'B2']))  #Landsat 8 SR\n",
    "left = ipyleaflet.TileLayer(url=jrc_tile_url)\n",
    "right=ipyleaflet.TileLayer(url=l8sr_tile_url)\n",
    "control = ipyleaflet.SplitMapControl(left_layer=left, right_layer=right)\n",
    "map1.add_control(control)\n",
    "map1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stratified Sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The JRC land water mask has a single band called **waterClass** which has the following classification\n",
    "\n",
    "* 0 nodata\n",
    "* 1 not water\n",
    "* 2 seasonal water\n",
    "* 3 permanent water\n",
    "\n",
    "In the following setup we are going to do use startified sampling to select permanent land and water. These points are then used to collect the reflectance values for each of these points as a table. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Land Mean for Bands:\n",
      "['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11', 'sr_aerosol', 'pixel_qa', 'radsat_qa']\n",
      "[230.12, 2962.62, 2942.94, 267.94, 402.78, 397.42, 1023.04, 838.32, 471.16, 322.0, 0.0, 97.2]\n",
      "\n",
      "water Mean for Bands:\n",
      "['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11', 'sr_aerosol', 'pixel_qa', 'radsat_qa']\n",
      "[174.96, 2953.84, 2934.94, 197.46, 298.24, 242.1, 63.08, 25.86, 17.48, 324.0, 0.0, 68.0]\n"
     ]
    }
   ],
   "source": [
    "#Select permanent water and not water\n",
    "water=global_water.first().clip(grid).gte(3)\n",
    "land=global_water.first().clip(grid).eq(1)\n",
    "\n",
    "def geomft(ft):\n",
    "    ftgeom=ft.setGeometry(ee.Geometry.Point([ft.get('longitude'), ft.get('latitude')]))\n",
    "    return ftgeom\n",
    "\n",
    "def strat(image,numpoints):\n",
    "    img = image.updateMask(image).addBands(ee.Image.pixelLonLat()).stratifiedSample(\n",
    "    numPoints=numpoints,\n",
    "    classBand='waterClass',\n",
    "    seed=1,\n",
    "    projection='EPSG:4326',\n",
    "    scale=30,\n",
    "    region= grid)\n",
    "    return img.map(geomft)\n",
    "\n",
    "stratified_water=strat(image=water,numpoints=50)\n",
    "stratified_land=strat(image=land,numpoints=50)\n",
    "\n",
    "# print(stratified_water.getInfo())\n",
    "# print(stratified_land.getInfo())\n",
    "\n",
    "# Get Landsat 8 SR band values for each class type from Sampling\n",
    "landMean = l8sr.filterBounds(grid).median().reduceRegion(\n",
    "  reducer= ee.Reducer.mean(), \n",
    "  geometry= stratified_land, \n",
    "  scale= 30\n",
    "  ).values()\n",
    "\n",
    "waterMean = l8sr.filterBounds(grid).median().reduceRegion(\n",
    "  reducer= ee.Reducer.mean(), \n",
    "  geometry= stratified_water, \n",
    "  scale= 30\n",
    "  ).values()\n",
    "\n",
    "print('Land Mean for Bands:'+'\\n'+str(l8sr.first().bandNames().getInfo()))\n",
    "print(landMean.getInfo())\n",
    "\n",
    "print('\\n'+'water Mean for Bands:'+'\\n'+ str(l8sr.first().bandNames().getInfo()))\n",
    "print(waterMean.getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating a Binary using an Index Approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The modified normalized difference water index ($MNDWI$) is a band ratio that is related to the amount of percentage water that is present in a pixel. For Landsat 8 the setup would be\n",
    "\n",
    "\\begin{equation*}\n",
    "MNDWI = \\frac{GREEN — SWIR}{GREEN + SWIR} = \\frac{Band3 — Band6}{Band3 + Band6}\n",
    "\\end{equation*}\n",
    "\n",
    "where $SWIR$ is the short wave infrared band and $Green$ is green band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "35bc27a1fc714e878d0d9c03664b4e25",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Create a slider widget to add both Landsat 8 and Landsat NDVI imagery\n",
    "img = l8sr.median()\n",
    "ndwi = img.normalizedDifference(['B3','B6']).rename('ndwi')\n",
    "\n",
    "map2 = ipyleaflet.Map(\n",
    "    center=(29.2963,-90.3446), zoom=12,\n",
    "    layout={'height':'500px'},\n",
    ")\n",
    "l8sr_tile_url = GetTileLayerUrl(l8sr.median().visualize(min=100, max=3500, gamma=1.5, bands= ['B5', 'B3', 'B2']))  #Landsat 8 SR\n",
    "ndwi_tile_url = GetTileLayerUrl(ndwi.visualize(bands=['ndwi'], min=-0.2, max=0.6, palette='brown,9999ff,blue'))  #Landsat 8 SR\n",
    "left = ipyleaflet.TileLayer(url=l8sr_tile_url)\n",
    "right=ipyleaflet.TileLayer(url=ndwi_tile_url)\n",
    "control = ipyleaflet.SplitMapControl(left_layer=left, right_layer=right)\n",
    "map2.add_control(control)\n",
    "map2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max Min NDVI for Stratified Sampling of Land [0.10585305105853052, -0.5666973321067157]\n",
      "Max Min NDVI for Stratified Sampling of Water [0.9204545454545454, 0.5058139534883721]\n"
     ]
    }
   ],
   "source": [
    "#Get min max region for the stratified sampling points for ndwi over land\n",
    "ndviland = ee.Image(ndwi).reduceRegion(\n",
    "  reducer= ee.Reducer.minMax(), \n",
    "  geometry= stratified_land, \n",
    "  scale= 30\n",
    "  ).values()\n",
    "\n",
    "ndviwater = ee.Image(ndwi).reduceRegion(\n",
    "  reducer= ee.Reducer.minMax(), \n",
    "  geometry= stratified_water, \n",
    "  scale= 30\n",
    "  ).values()\n",
    "\n",
    "#pp.pprint(zonal.getInfo())\n",
    "print('Max Min NDVI for Stratified Sampling of Land '+str(ndviland.getInfo()))\n",
    "print('Max Min NDVI for Stratified Sampling of Water '+str(ndviwater.getInfo()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'type': 'Image', 'bands': [{'id': 'ndwi', 'data_type': {'type': 'PixelType', 'precision': 'int', 'min': 0, 'max': 1}, 'crs': 'EPSG:4326', 'crs_transform': [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]}]}\n"
     ]
    }
   ],
   "source": [
    "# From the earlier min max we find that we can safely assume that NDVI >0 equals land (assumption)\n",
    "constrained= ee.Image(ndwi).clip(grid).gt(0.2)\n",
    "print(constrained.getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d9124bef2ffe40dabd7e63f2f92d71ea",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(basemap={'url': 'https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', 'max_zoom': 19, 'attribution': 'Map …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Create a slider widget to add both Landsat 8 and Landsat NDVI imagery\n",
    "img = l8sr.median()\n",
    "\n",
    "map3 = ipyleaflet.Map(\n",
    "    center=(29.2963,-90.3446), zoom=12,\n",
    "    layout={'height':'500px'},\n",
    ")\n",
    "l8sr_tile_url = GetTileLayerUrl(l8sr.median().clip(grid).visualize(min=100, max=3500, gamma=1.5, bands= ['B5', 'B3', 'B2']))  #Landsat 8 SR\n",
    "constrained_tile_url = GetTileLayerUrl(constrained.visualize(bands=['ndwi'],palette='brown,blue'))  #NDWI constrained\n",
    "left = ipyleaflet.TileLayer(url=l8sr_tile_url)\n",
    "right=ipyleaflet.TileLayer(url=constrained_tile_url)\n",
    "control = ipyleaflet.SplitMapControl(left_layer=left, right_layer=right)\n",
    "map3.add_control(control)\n",
    "map3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We learned how to \n",
    "* Filter and mask a collection\n",
    "* Use a binary mask to create stratified sampling points\n",
    "* Create NDWI to get percentage water and get min max values using the Stratified sampling points\n",
    "* Generate a binary land water classification using Normalized Difference Water Index (NDWI) thresholds"
   ]
  }
 ],
 "metadata": {
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}