docs/SpatialFilteringGEE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GEOG0027 (2022/2023) Spatial Filtering with Google Earth Engine (GEE)\n",
    "\n",
    "This notebook uses the image exploration scripts from the [Intro Google Earth Engine](docs/Intro_to_GEE.ipynb) practical, and introduces spatial filtering functions to the images. First complete the introduction practical before attempting this one.\n",
    "\n",
    "Again we recommend using UCL's [Jupyter Hub](https://jupyter.data-science.rc.ucl.ac.uk/) (VPN required) to work with this chapter. To download this notebook, to execute and modify the code, run the following command in a `Terminal`:\n",
    "\n",
    "`wget https://github.com/UCL-EO/GEOG0027/raw/main/docs/SpatialFilteringGEE.ipynb`\n",
    "\n",
    "As with the introduction practical, modify the local lon/lat coordinates to get the location you desire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ca2b921f98964872b31f93276bd7d76e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[51.49163903397572, -0.08313179193795397], controls=(WidgetControl(options=['position', 'transparen…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import geemap, ee, os, numpy\n",
    "import ipyleaflet\n",
    "\n",
    "loc_coords = [51.49163903397572, -0.08313179193795397] ### EDIT THIS LINE WITH YOUR COORDS, [Lat,Lon] format\n",
    "\n",
    "Map = geemap.Map(center=loc_coords, zoom=9)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A basic Google Map-like interface should be displayed here now. If you can't see anything, please ensure that the ipyleaflet nbextension has been enabled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the Area Of Interest (AOI)\n",
    "Now that we have a location to centre the Map on, we need to determine the size of the map. This can be done using a range in latitude/longitude using the function below. This cell just defines how a function works, and does not perform any actual calculations. To do this we then call the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def expand_coords(centre,lon_expand = 1.0,lat_expand = 0.5):\n",
    "    lat1 = centre[0]-lat_expand\n",
    "    lat2 = centre[0]+lat_expand\n",
    "    lon1 = centre[1]-lon_expand\n",
    "    lon2 = centre[1]+lon_expand\n",
    "    return [lon1,lat1,lon2,lat2,]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This next piece of code will then calculate the map size. If the map isn't the correct size for your liking, add into the function alternative settings for 'lon_expand' and 'lat_expand' from those set in the previous cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec_coords = expand_coords(loc_coords)\n",
    "# rec_coords = expand_coords(loc_coords,lat_expand=0.7) ## example adjusted rectangle size\n",
    "local_rec = ee.Geometry.Rectangle(rec_coords) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Landsat data collections from GEE\n",
    "Now we can see the the buffered AOI displayed on the Map. Next, let's load some Landsat images for the elected area. I've defined here a python function called `display_landsat_collection` to do so. It automatically loads both the [surface reflectance](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR) and [annual NDVI](https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_ANNUAL_NDVI) image collections from GEE's data catalog and also calculates the annual means for each band. \n",
    "\n",
    "You can skip most of the details of what's inside the code cell, but only to look at the first (and last) line of code. In order to run such function, you will need to choose a year (any year since 1984) and an AOI. In the following example, I choose year 2019 and the Shenzhen buffer to demonstrate the use of the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "landsat_vis_param = {\n",
    "            'min': 0,\n",
    "            'max': 3000,\n",
    "            'bands': ['NIR', 'Red', 'Green']  # False Colour Composit bands to be visualised \n",
    "}\n",
    "RGB_vis_param = {\n",
    "            'min': 0,\n",
    "            'max': 3000,\n",
    "            'bands': [ 'Red', 'Green', 'Blue']  # False Colour Composit bands to be visualised \n",
    "}\n",
    "ndvi_colorized_vis = {\n",
    "            'min': 0.0,\n",
    "            'max': 1.0,\n",
    "            'palette': [\n",
    "            'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',\n",
    "            '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',\n",
    "            '012E01', '011D01', '011301']\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_landsat_collection(inmap,year, aoi, cloud_tolerance = 3.0, \n",
    "                            DISPLAY_ON_MAP = False, MEDIAN_ONLY = False):\n",
    "    '''This function allows GEE to display a Landsat data collection \n",
    "    from any year between 1984 and present year\n",
    "    that fall within the AOI and cloud tolerance, e.g. 3.0%.\n",
    "    There are two optional flag:\n",
    "    When DISPLAY_ON_MAP is TRUE, display this layer onto Map;\n",
    "    When return_series = 'MEDIAN_ONLY', only median SR is loaded into landsat_ts, and\n",
    "    Setting this option to MEDIAN_ONLY would be faster than loading other collections. \n",
    "    '''\n",
    "    assert year >= 1984\n",
    "    \n",
    "    def renameBandsETM(image):\n",
    "        if year >=2013: #LS8\n",
    "            bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7'] #, 'pixel_qa'\n",
    "            new_bands = ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'] #, 'pixel_qa'\n",
    "        elif year <=1984:\n",
    "            bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa']\n",
    "            new_bands = ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']\n",
    "        return image.select(bands).rename(new_bands)\n",
    "        \n",
    "    if not(MEDIAN_ONLY):\n",
    "        if year >= 2013:\n",
    "            layer_name = 'LC08' # LS8: 2013-now        \n",
    "        elif year == 2012: # # LS7: 1999- , however SLC error >= 1999:\n",
    "            layer_name = 'LE07' \n",
    "        elif year >=1984:\n",
    "            layer_name = 'LT05' # LS5: 1984-2012\n",
    "       \n",
    "        collection_name_sr = f\"LANDSAT/{layer_name}/C01/T1_SR\" \n",
    "        # You can also use the following line, if interested in incorperating ndvi\n",
    "        collection_name_ndvi = f\"LANDSAT/{layer_name}/C01/T1_ANNUAL_NDVI\" \n",
    "\n",
    "        all_sr_image = ee.ImageCollection(collection_name_sr) \\\n",
    "            .filterBounds(aoi) \\\n",
    "            .filterDate(f'{year}-01-01', f'{year}-12-31') \\\n",
    "            .filter(ee.Filter.lt('CLOUD_COVER', cloud_tolerance))\\\n",
    "            .sort('system:time_start') \\\n",
    "            .select('B[1-7]') \\\n",
    "            .sort('CLOUD_COVER')\n",
    "\n",
    "        all_sr_image = all_sr_image.map(renameBandsETM) # rename bands with 'renameBandsETM' function\n",
    "        \n",
    "        # reduce all_sr_image to annual average per pixel\n",
    "        mean_image = all_sr_image.mean()\n",
    "        mean_image = mean_image.clip(aoi).unmask()\n",
    "\n",
    "        ndvi_image = ee.ImageCollection(collection_name_ndvi)\\\n",
    "            .filterBounds(aoi) \\\n",
    "            .filterDate(f'{year}-01-01', f'{year}-12-31')\\\n",
    "            .select('NDVI')\\\n",
    "            .median()\n",
    "        ndvi_image = ndvi_image.clip(aoi).unmask()\n",
    "\n",
    "        #mean_image.addBands(ndvi_image, 'NDVI')\n",
    "    \n",
    "    # This line loads all annual median surface ref\n",
    "    landsat_ts = geemap.landsat_timeseries(roi=aoi, start_year=year, end_year=year, \\\n",
    "                                       start_date='01-01', end_date='12-31')\n",
    "\n",
    "    median_image = landsat_ts.first().clip(aoi).unmask()\n",
    "    \n",
    "    if DISPLAY_ON_MAP == True:\n",
    "        \n",
    "        if not(MEDIAN_ONLY):\n",
    "            inmap.addLayer(ndvi_image, ndvi_colorized_vis, 'NDVI '+str(year),  opacity=0.9)\n",
    "            inmap.addLayer(mean_image, landsat_vis_param, \"Mean Ref \"+str(year))\n",
    "        inmap.addLayer(median_image, landsat_vis_param, \"Median Ref \"+str(year))\n",
    "\n",
    "    if MEDIAN_ONLY:\n",
    "        return median_image\n",
    "    else:\n",
    "        return all_sr_image, mean_image, median_image, ndvi_image "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the `load_landsat_collection` function has been defined, and we will run/execute it by calling the function name with appropriate input parameters (or 'arguments). The output of such function will be returned to the variables on the LHS of the equal sign, i.e. all_image_2019, mean_2019, median_2019, and ndvi_2019 in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2c4ea71ada24895adefda441145bd66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[51.49163903397572, -0.08313179193795397], controls=(WidgetControl(options=['position', 'transparen…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map = geemap.Map(center=loc_coords, zoom=10)\n",
    "\n",
    "# All you need to modify here is the YEAR below:\n",
    "all_image_2019, mean_2019, median_2019, ndvi_2019 = load_landsat_collection(Map,2019,\\\n",
    "                                        local_rec, cloud_tolerance = 2,\\\n",
    "                                        DISPLAY_ON_MAP = True)\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should examine the functions under the `toolbar` and `layer` buttons on the top-right corner of the Map, e.g. use the `inspector` and `plotting` tools to check the data values, or use `layers` control to switch on/off layers or to adjust opacity.\n",
    "\n",
    "We can also check the metadata from the Landsat image collection we just loaded from the Cloud. Have a look of the output. Any useful information?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'CLOUD_COVER': 0.01, 'CLOUD_COVER_LAND': 0.03, 'EARTH_SUN_DISTANCE': 0.989582, 'ESPA_VERSION': '2_23_0_1b', 'GEOMETRIC_RMSE_MODEL': 8.908, 'GEOMETRIC_RMSE_MODEL_X': 6.443, 'GEOMETRIC_RMSE_MODEL_Y': 6.152, 'IMAGE_DATE': '2019-02-24', 'IMAGE_QUALITY_OLI': 9, 'IMAGE_QUALITY_TIRS': 9, 'LANDSAT_ID': 'LC08_L1TP_201025_20190224_20190308_01_T1', 'LEVEL1_PRODUCTION_DATE': 1552088520000, 'NOMINAL_SCALE': 30, 'PIXEL_QA_VERSION': 'generate_pixel_qa_1.6.0', 'SATELLITE': 'LANDSAT_8', 'SENSING_TIME': '2019-02-24T10:52:30.0138839Z', 'SOLAR_AZIMUTH_ANGLE': 157.231583, 'SOLAR_ZENITH_ANGLE': 62.337173, 'SR_APP_VERSION': 'LaSRC_1.3.0', 'WRS_PATH': 201, 'WRS_ROW': 25, 'system:asset_size': '392.213325 MB', 'system:band_names': ['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'], 'system:id': 'LANDSAT/LC08/C01/T1_SR/LC08_201025_20190224', 'system:index': 'LC08_201025_20190224', 'system:time_end': '2019-02-24 10:52:30', 'system:time_start': '2019-02-24 10:52:30', 'system:version': 1563396284578517}\n"
     ]
    }
   ],
   "source": [
    "first_image = all_image_2019.first() \n",
    "\n",
    "props = geemap.image_props(first_image)\n",
    "print( props.getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, examine the mean, median surface reflectance and/or NDVI layers we've visualized. Which one is better? How do these relate to the individual images?\n",
    "\n",
    "One particular issue that can arrise is taking the 'mean' of individiual NDVI images. Instead we can take the 'NDVI' of mean images (or medians).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to switching layers on and off, adjusting opacity, we can also use python code to some simple mathmatical operations, e.g. calculating the differences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2c4ea71ada24895adefda441145bd66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(bottom=87463.0, center=[51.49163903397572, -0.08313179193795397], controls=(WidgetControl(options=['positi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "median_ndvi = median_2019.normalizedDifference(['NIR', 'Red'])\n",
    "median_ndwi = median_2019.normalizedDifference(['Green','NIR'])\n",
    "\n",
    "Map.addLayer(median_ndvi, ndvi_colorized_vis, 'NDVI from Median LS')\n",
    "Map.addLayer(median_ndwi, ndvi_colorized_vis, 'NDWI from Median LS')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying kernels to the data\n",
    "Here we apply the spatial filtering techniques previously performed using ENVI (see [here](docs/SpatialFiltering.ipynb)) on the GEE images. The same method of convolving a kernel over an image is available, though we have to use the cloud based functional methods of Google Earth Engine. Some example kernels are shown below. What do these kernels do to the displayed data?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2c4ea71ada24895adefda441145bd66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(bottom=87463.0, center=[51.49163903397572, -0.08313179193795397], controls=(WidgetControl(options=['positi…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### high and low pass 3x3 kernels as we used in ENVI. \n",
    "low_pass_K=ee.Kernel.fixed(3,3,[[1/9,1/9,1/9],[1/9,1/9,1/9],[1/9,1/9,1/9]],1,1,False)\n",
    "high_pass_K=ee.Kernel.fixed(3,3,[[-1,-1,-1],[-1,8,-1],[-1,-1,-1]],1,1,False)\n",
    "#### circle averaging kernel - in a functional form, but more flexible (try adjusting the radius)\n",
    "low_pass_C=ee.Kernel.circle(radius = 5,units = 'pixels')\n",
    "\n",
    "Map.addLayer(median_2019.convolve(high_pass_K),landsat_vis_param,'FCC high passed')\n",
    "Map.addLayer(median_2019.convolve(low_pass_K),landsat_vis_param,'FCC low passed')\n",
    "Map.addLayer(median_2019.convolve(low_pass_C),landsat_vis_param,'FCC low passed (circle)')\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">## Advanced Challenge\n",
    ">Can you copy the cell above, and adjust the code to add a high passed version of a visual coloured median image? You'll need the `RGB_vis_param` display specifications, and the `median_2019` image.\n",
    ">\n",
    "> Can you create a low passed version of the NDVI image?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping regions with high Vegetation or Water Index\n",
    "When we used ENVI we used ranges of pixel values to select certain areas. This can be done with GEE too, although it is a littel more tricky. The code below will create maps that split the pixels depending upon:\n",
    "- NDVI > 0.5\n",
    "- NDWI > 0.0\n",
    "\n",
    "And then by combining the opposite of the previous two maps\n",
    "- NDVI < 0.5 AND NDWI < 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2c4ea71ada24895adefda441145bd66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(bottom=1395239.0, center=[51.493995984702856, -0.15440940856933594], controls=(WidgetControl(options=['pos…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### NDVI > 0.5\n",
    "ndvi_gt = median_ndvi.gt(median_ndvi.constant(0.5))\n",
    "Map.addLayer(ndvi_gt,{},'NDVI > 0.5')\n",
    "\n",
    "### NDWI > 0.0\n",
    "ndwi_gt=median_ndwi.gt(median_ndwi.constant(0.0))\n",
    "Map.addLayer(ndwi_gt,{},'NDWI > 0.0')\n",
    "\n",
    "### Opposite of BOTH\n",
    "both_gt=ndwi_gt.Not().And(ndvi_gt.Not())\n",
    "Map.addLayer(both_gt,{},'NDVI < 0.5, AND NDWI < 0.0')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Majority Analysis\n",
    "For a highly variable finely detailed classification, such as we performed above, the resulting maps can appear noisey, and contain data that is not useful for large scale mean metrics. We can remove this detail by applying a majority kernel. The code below will first apply this to the NDVI > 0.5 map. \n",
    ">Can you adjust the code (add and remove `#` characters) to apply the different kernels?\n",
    ">\n",
    ">Can you apply the majority kernel to the BOTH map created previously?\n",
    "\n",
    "Compare the maps using the layer selection tool, how do they compare at different zoom levels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d2c4ea71ada24895adefda441145bd66",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(bottom=1395239.0, center=[51.493995984702856, -0.15440940856933594], controls=(WidgetControl(options=['pos…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "majority_ndvi = ndvi_gt.reduceNeighborhood(\n",
    "                reducer=ee.Reducer.mode(), \n",
    "                kernel=ee.Kernel.fixed(3,3,[[1]*3]*3,1,1,False) #3x3\n",
    "#                 kernel=ee.Kernel.fixed(9,9,[[1]*9]*9,5,5,False) #9x9\n",
    "#                 kernel=ee.Kernel.circle(4)\n",
    "                )\n",
    "Map.addLayer(majority_ndvi,{},'NDVI > 0.5, 3x3 majority')\n",
    "Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}