Coastal_erosion_monitoring.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monitoring coastal change and erosion with Landsat 8\n",
    "\n",
    "### Background\n",
    "Over 40% of the world’s population lives within 100 km of the coastline. However, coastal environments are constantly changing, with erosion and coastal change presenting a major challenge to valuable coastal infrastructure and important ecological habitats. Up-to-date data on coastal change and erosion is essential for coastal managers to be able to identify and minimise the impacts of coastal change and erosion.\n",
    "\n",
    "### The Problem\n",
    "Monitoring coastlines and rivers using field surveys can be challenging and hazardous, particularly at regional or national scale. Aerial photography and LiDAR can be used to monitor coastal change, but this is often expensive and requires many repeated flights over the same areas of coastline to build up an accurate history of how the coastline has changed across time.\n",
    "\n",
    "### Digital Earth Australia use case\n",
    "Imagery from satellites such as the NASA/USGS Landsat program is available for free for the entire planet, making satellite imagery a powerful and cost-effective tool for monitoring coastlines and rivers at regional or national scale. By identifying and extracting the precise boundary between water and land based on satellite data, it is possible to extract accurate shorelines that can be compared across time to reveal hotspots of erosion and coastal change.\n",
    "\n",
    "The usefulness of satellite imagery in the coastal zone can be affected by the presence of clouds, sun-glint over water, poor water quality (e.g. sediment) and the influence of tides. The effect of these factors can be reduced by combining individual noisy images into cleaner \"summary\" or composite layers, and filtering the data to focus only on images taken at certain tidal conditions (e.g. high tide). \n",
    "\n",
    "In this example, we combine data from the Landsat 8 satellite with image compositing and tide filtering techniques to accurately map shorelines across time, and identify areas that have changed significantly between 2013 and 2018. The worked example takes users through the code required to:\n",
    "\n",
    "1.  Load in a cloud-free Landsat time series\n",
    "2.  Compute a water index (MNDWI)\n",
    "3.  Filter images by tide height\n",
    "4.  Create \"summary\" or composite images for given time periods\n",
    "5.  Extract and visualise shorelines across time\n",
    "\n",
    "### Technical details\n",
    "\n",
    "* Products used: `ls8_nbar_scene`\n",
    "* Bands used: `red`, `green`, `blue` and `swir1`. You can read more about the bands available for Landsat 8 [here](https://landsat.usgs.gov/what-are-best-spectral-bands-use-my-study)\n",
    "* Analyses used: MNDWI water index, tide filtering, compositing, waterline extraction\n",
    "\n",
    "**To run this analysis, run all the cells in the notebook, starting with the \"Load modules\" cell. When you finished the analysis, you can return to the \"Set up analysis\" cell, modify some values (e.g. choose a different location, time, tide height range or time period to analyse) and re-run the analysis. There are additional instructions on modifying the notebook at the end.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datacube\n",
    "import sys\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import waterline_funcs\n",
    "from utils.display import display_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up analysis\n",
    "\n",
    "The following lines allow us to set up our analysis by selecting the location (`lat_range` and `lon_range`) and time (`time_range`) we want to investigate. We can also set two other important parameters: `tide_range` allows us to extract shorelines based on satellite images taken only at certain part of the tide range (e.g. at low or high tide). This can be useful for seperating the effect of erosion from the influence of tides. Finally, the `time_step` parameter allows us to choose the length of the time periods we want to compare: e.g. shorelines for each year, or shorelines for each six months etc:\n",
    "\n",
    "* **lat_range:** The latitude range to analyse (e.g. `(-12.20, -12.27)`). For fast load times, keep this to a range of ~0.1 degrees or less.\n",
    "\n",
    "* **lon_range:** The longitude range to analyse (e.g. `(131.80, 131.91)`). For fast load times, keep this to a range of ~0.1 degrees or less.\n",
    "\n",
    "* **time_range:** The date range to analyse (e.g. `('2013-01-01', '2018-12-31') `)\n",
    "\n",
    "* **time_step:** The time period to compare coastline change. `1Y` will generate one coastline for every year in the dataset; `6M` will produce a coastline for every six months, etc.\n",
    "\n",
    "* **tide_range:** The lowest and highest tide heights to analyse (e.g. `(0.00, 4.00)` in metre units). This will remove all images with tides outside this range from the dataset, allowing you to produce shorelines for specific tidal conditions (e.g. low tide, average tide, high tide etc)\n",
    "\n",
    "* **tide_file:** A text file containing a list of times and tide heights for a location near to the study area. These tides are used to filter the satellite dataset to focus on images taken at certain tide conditions. Tide height files are stored in the `tide_heights` folder and are named using the following convention: `locationname_latitude_longitude_tides.csv`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Edit these lines to change any of the analysis parameters\n",
    "\n",
    "If you're running this notebook for the first time, we recommend you keep the default settings below. This will allow you to understand how the analysis works.\n",
    "\n",
    "The example we've selected looks at a section of the Daly River in the Northen Territory, Australia.\n",
    "\n",
    "Run the following cell to set values for the various parameters discussed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat_range = (-13.39, -13.47) \n",
    "lon_range = (130.27, 130.35) \n",
    "time_range = ('2013-01-01', '2018-12-31') \n",
    "time_step = '1Y'\n",
    "\n",
    "tide_range = (0.00, 3.50)\n",
    "tide_file = 'tide_heights/dalyriver_-13.32_130.23_tides.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View the selected location\n",
    "\n",
    "The next cell will display the selected area on an interactive map. Feel free to zoom in and out to get a better understanding of the area you'll be analysing. Clicking on any point of the map will reveal the latitude and longitude coordinates of that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "display_map(latitude=lat_range, longitude=lon_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load cloud-masked Landsat data\n",
    "The first step in this analysis is to load in Landsat data for the `lat_range`, `lon_range` and `time_range` we provided above. The code below first connects to the datacube database, and then uses the `load_cloudmaskedlandsat` function to load in data from the Landsat 8 satellite for the area and time included in `lat_range`, `lon_range` and `time_range`. The function will also automatically mask out clouds from the dataset, allowing us to focus on pixels that contain useful data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connect to the datacube\n",
    "dc = datacube.Datacube(app='Coastal erosion')\n",
    "\n",
    "# Create the 'query' dictionary object, which contains the longitudes, latitudes and time provided above\n",
    "query = {\n",
    "    'y': lat_range,\n",
    "    'x': lon_range,\n",
    "    'time': time_range,\n",
    "    'output_crs': 'EPSG:28352',\n",
    "    'resolution': (-25, 25)\n",
    "}\n",
    "\n",
    "# Load Landsat 8 data for the time and area in the query. This may take several minutes, so please be patient.\n",
    "landsat_ds = waterline_funcs.load_cloudmaskedlandsat(\n",
    "    dc=dc,\n",
    "    query=query,\n",
    "    platform='ls8',\n",
    "    bands=['red', 'green', 'blue', 'swir1']\n",
    ")\n",
    "\n",
    "print(landsat_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot example timestep in true colour\n",
    "To verify that we extracted data for the area we're interested in analysing, we can plot up the first image in the dataset in true colour. This makes use of the `red`, `green` and `blue` Landsat 8 bands that we loaded in the previous cell.\n",
    "\n",
    "Change `index=0` to another value (e.g. `index=1`) to plot up different timesteps. You should be able to clearly see where clouds have been masked out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "waterline_funcs.rgb(landsat_ds, index=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Modified Normalised Difference Water Index\n",
    "To extract shoreline locations, we need to be able to seperate water from land in our study area. To do this, we can use our Landsat data to calculate a water index called the `Modified Normalised Difference Water Index`, or MNDWI. This index uses the ratio of green and mid-infrared radiation to identify the presence of water. You can find more information on it in the [original research article](http://www.aari.ru/docs/pub/060804/xuh06.pdf). The formula is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\text{MNDWI} &= \\frac{(Green - MIR)}{(Green + MIR)},\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where `Green` is the green band and `MIR` is the mid-infrared band. For Landsat 8, we can use the Short-wave Infrared (SWIR) 1 band as our measure for `MIR`.\n",
    "\n",
    "When it comes to interpreting the index, <font color='blue'>**High values (greater than 0, blue colours) typically represent water pixels**</font>, while <font color='red'>**low values (less than 0, red colours) represent land**</font>. You can use the cell below to calculate and plot one of the images after calculating the index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "landsat_ds['MNDWI'] = (landsat_ds.green - landsat_ds.swir1)/(landsat_ds.green + landsat_ds.swir1)\n",
    "landsat_ds.MNDWI.isel(time=0).plot(cmap='RdBu', size=6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does the plot of the index compare to the optical image from earlier? Was there water or land anywhere you weren't expecting?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import tide height data\n",
    "The location of the shoreline can vary greatly from low to high tide. In the code below, we aim to reduce the effect of tides by loading in tide height data, and keeping only the satellite images that were taken at specific tidal conditions. For example, if `tide_range = (0.00, 2.00)`, we are telling the analysis to focus only on satellite images taken when the tide was between `0.00 m` and `2.00 m` relative to Mean Sea Level.\n",
    "\n",
    "For this example notebook, we will use tide data for a single location off the coast of Darwin, Northern Territory. You can load the existing tide height data in the following cell using the `pandas` module, which we imported earlier. The code will read the `tide_file` variable provided during the set up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tide_data = pd.read_csv(tide_file, parse_dates=['time'], index_col='time')\n",
    "tide_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the tide height data, we need to estimate the tide height for each of our Landsat 8 images. We can do this by interpolating between the data points we do have to get the tide height when each image was taken."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, we convert the data to an xarray dataset so we can analyse it in the same way as our Landsat data\n",
    "tide_data_xr = tide_data.to_xarray()\n",
    "\n",
    "# We want to convert our hourly tide heights to estimates of exactly how high the tide was at the time that\n",
    "# each satellite image was taken. To do this, we can use `.interp` to 'interpolate' a tide height for each\n",
    "# Landsat timestamp:\n",
    "landsat_tideheights = tide_data_xr.interp(time=landsat_ds.time)\n",
    "\n",
    "# We then want to put these values back into the Landsat dataset so that each image has an estimated tide height:\n",
    "landsat_ds['tide_height'] = landsat_tideheights.tide_height\n",
    "\n",
    "# Plot the resulting tide heights for each Landsat image:\n",
    "landsat_ds.tide_height.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot above should make it clear that limiting the range of the tides for the analysis should give you more consistent results. Consistency is critical as you want to compare the change in the shoreline from year to year, and large variance in the tide height could obscure your results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Filter Landsat images by tide height\n",
    "Here we take the Landsat dataset and only keep the images with tide heights we want to analyse (i.e. tides within the heights given by `tide_range`). This will result in a smaller number of images (e.g. 67 images compared to 129):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "landsat_hightide = landsat_ds.where((landsat_ds.tide_height > tide_range[0]) & \n",
    "                                    (landsat_ds.tide_height < tide_range[1]), drop=True)\n",
    "print(landsat_hightide)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create water index summary images\n",
    "Individual remote sensing images can be affected by noisy data, including clouds, sunglint and poor water quality conditions (e.g. sediment). To produce cleaner images that can be compared more easily across time, we can create 'summary' images or composites that combine multiple images into one image to reveal the 'typical' or median appearance of the landscape for a certain time period. In this case, we use the median as the summary statistic because it prevents strong outliers (like masked cloud values) from skewing the data, which would not be the case if we were to use the mean.\n",
    "\n",
    "In the code below, we take the time series of images and combine them into single images for each `time_step`. For example, if `time_step = '1Y'`, the code will produce one new image for each year of data in the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine by time\n",
    "landsat_resampled = landsat_hightide.MNDWI.resample(time=time_step).median('time')\n",
    "landsat_resampled.plot(col='time', cmap='RdBu', col_wrap=3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract shorelines from imagery\n",
    "We now want to extract an accurate shoreline for each each of the summary images above (e.g. 2013, 2014 etc. summaries). The code below identifies the boundary between land and water by tracing a line along pixels with a water index value of `0` (halfway between land and water water index values). It returns a shapefile with one line for each time step:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up attributes to assign to each waterline\n",
    "attribute_data = {'time': [str(i)[0:10] for i in landsat_resampled.time.values]}\n",
    "attribute_dtypes = {'time': 'str'}\n",
    "\n",
    "# Extract waterline contours for the '0' water index threshold:\n",
    "contour_gdf = waterline_funcs.contour_extract(\n",
    "    z_values=[0],\n",
    "    ds_array=landsat_resampled,\n",
    "    ds_crs=landsat_ds.crs,\n",
    "    ds_affine=landsat_ds.geobox.transform,\n",
    "    output_shp=f'output_waterlines.shp',\n",
    "    attribute_data=attribute_data,\n",
    "    attribute_dtypes=attribute_dtypes,\n",
    "    min_vertices=100\n",
    ")\n",
    "\n",
    "# Plot output shapefile\n",
    "contour_gdf.plot(cmap='YlOrRd', figsize=(15, 15))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above plot is a basic visualisation of the contours returned by the `contour_extract` function. Given we now have the shapfile, we can use a more complex function to make an interactive plot for viewing the change in shoreline over time below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot interactive map of output shorelines coloured by time\n",
    "The next cell provides an interactive map with an overlay of the shorelines identified in the previous cell. Run it to view the map.\n",
    "\n",
    "Zoom in to the map below to explore the resulting set of shorelines. Older shorelines are coloured in yellow; more recent shorelines in red (click on the shorelines to return their time period, which will appear below the map). Using this data, we can easily identify areas of rivers or coastline that have changed significantly over time, or areas that have remained stable over the past 5 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "waterline_funcs.map_shapefile(gdf=contour_gdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drawing conclusions\n",
    "Here are some questions to think about:\n",
    "* What can you conclude about the change in the shoreline? \n",
    "* Which sections of the shoreline have seen the most change? \n",
    "* Is the change consistent with erosion? \n",
    "* What other information might you need to draw additional conclusions about the cause of the change?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "When you are done, return to the \"Set up analysis\" cell, modify some values and rerun the analysis. Some suggestions for things you can change are:\n",
    "1. Time range\n",
    "1. Tide height range\n",
    "1. Time step\n",
    "1. Location\n",
    "\n",
    "If you're going to change the location, you'll need to make sure Landsat 8 data is available for the new location, which you can check at the [DEA Dashboard](https://dashboard.dea-sandbox.test.frontiersi.io/ls8_nbar_scene). You will  need to find tide heights for your region and store them in the same format as the `tide_heights/dalyriver_-13.32_130.23_tides.csv`. You will also need to update the `tide_file` variable in the set up."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further exploration\n",
    "While this notebook mostly utilised the functions available from the Open Data Cube library, some of the steps were conducted with functions from the `waterline_funcs` module, which you can find in the `utils` folder. If you're curious about how these functions work, feel free to explore the `utils/waterline_funcs.py` file. This option is for users who feel comfortable programming in Python. You'll be able to navigate to this file from the Jupyter file explorer page."
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}