notebooks/DEACoastlines_generation_raster.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DEA Coastlines raster generation <img align=\"right\" src=\"https://github.com/GeoscienceAustralia/dea-notebooks/raw/develop/Supplementary_data/dea_logo.jpg\">\n",
    "\n",
    "This code conducts raster generation for DEA Coastlines:\n",
    "\n",
    "* Load stack of all available Landsat 5, 7, 8 and 9 satellite imagery for a location using [ODC Virtual Products](https://docs.dea.ga.gov.au/notebooks/Frequently_used_code/Virtual_products.html)\n",
    "* Convert each satellite image into a remote sensing water index (MNDWI)\n",
    "* For each satellite image, model ocean tides into a grid based on exact time of image acquisition\n",
    "* Interpolate tide heights into spatial extent of image stack using the [FES2014 global tide model](https://github.com/GeoscienceAustralia/dea-coastlines/wiki/Setting-up-tidal-models-for-DEA-Coastlines)\n",
    "* Mask out high and low tide pixels by removing all observations acquired outside of 50 percent of the observed tidal range centered over mean sea level\n",
    "* Combine tidally-masked data into annual median composites representing the most representative position of the coastline at approximately mean sea level each year\n",
    "\n",
    "This is an interactive version of the code intended for prototyping; to run this analysis at scale, use the [command line tools](DEACoastlines_generation_CLI.ipynb).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "Set working directory to top level of repository to ensure links work correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load packages\n",
    "\n",
    "First we import the required Python packages, then we connect to the database, and load the catalog of virtual products."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pip install -r requirements.in --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext line_profiler\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "import sys\n",
    "import dask\n",
    "import xarray as xr\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import partial\n",
    "from datacube.utils.geometry import Geometry\n",
    "\n",
    "# Load DEA Coastlines and DEA tools code\n",
    "import coastlines.raster\n",
    "import coastlines.utils\n",
    "from dea_tools.dask import create_local_dask_cluster\n",
    "from dea_tools.coastal import model_tides, pixel_tides\n",
    "from dea_tools.datahandling import parallel_apply\n",
    "\n",
    "# Connect to datacube\n",
    "import datacube\n",
    "import odc.geo.xr\n",
    "\n",
    "dc = datacube.Datacube(app=\"DEACoastlines\")\n",
    "\n",
    "# Create local dask client for parallelisation\n",
    "client = create_local_dask_cluster(return_client=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Setup\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set analysis parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# study_area = 8\n",
    "# raster_version = 'development'\n",
    "# start_year = 1988\n",
    "# end_year = 2021\n",
    "# config_path = 'configs/dea_coastlines_config_development.yaml'\n",
    "\n",
    "study_area = 283\n",
    "raster_version = 'testing2022'\n",
    "start_year = 1988\n",
    "end_year = 2022\n",
    "config_path = 'configs/dea_coastlines_config.yaml'\n",
    "\n",
    "# Load analysis params from config file\n",
    "config = coastlines.utils.load_config(\n",
    "    config_path=config_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load supplementary data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Albers grid cells used to process the analysis\n",
    "gridcell_gdf = (\n",
    "    gpd.read_file(config['Input files']['grid_path']).to_crs(\n",
    "        epsg=4326).set_index('id'))\n",
    "gridcell_gdf.index = gridcell_gdf.index.astype(int).astype(str)\n",
    "gridcell_gdf = gridcell_gdf.loc[[str(study_area)]]\n",
    "gridcell_gdf.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading data\n",
    "### Create spatiotemporal query\n",
    "This establishes the spatial and temporal extent used to search for Landsat satellite data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create query based on analysis gridcell extent\n",
    "geopoly = Geometry(gridcell_gdf.iloc[0].geometry, crs=gridcell_gdf.crs)\n",
    "query = {\n",
    "    'geopolygon': geopoly.buffer(0.05),\n",
    "    'time': (str(start_year - 1), str(end_year + 1)),\n",
    "    'dask_chunks': {\n",
    "        'time': 1,\n",
    "        'x': 2048,\n",
    "        'y': 2048\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load satellite data as MNDWI\n",
    "This step loads satellite data from Landsat 5, 7, 8 and 9, and returns the data as a cloud-masked array converted to the Modified Normalised Difference Water Index (MNDWI).\n",
    "For Digital Earth Australia Coastlines, this is achieved using [ODC Virtual Products](https://docs.dea.ga.gov.au/notebooks/Frequently_used_code/Virtual_products.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load MNDWI virtual product\n",
    "ds = coastlines.raster.load_water_index(\n",
    "    dc,\n",
    "    query,\n",
    "    yaml_path=config[\"Virtual product\"][\"virtual_product_path\"],\n",
    "    product_name=config[\"Virtual product\"][\"virtual_product_name\"],\n",
    "    mask_terrain_shadow=False,\n",
    ")\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tidal modelling\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Interpolate tides into each satellite timestep\n",
    "For each satellite timestep, model tide heights into a low-resolution 5 x 5 km grid (matching resolution of the FES2014 tidal model), then reproject modelled tides into the spatial extent of our satellite image. Add  this new data as a new variable in our satellite dataset to allow each satellite pixel to be analysed and filtered/masked based on the tide height at the exact moment of satellite image acquisition. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds[\"tide_m\"], tides_lowres = pixel_tides(ds, resample=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot example interpolated tide surface for a single timestep:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot\n",
    "timestep = 15\n",
    "ds_i = ds[\"tide_m\"].isel(time=timestep)\n",
    "ds_lowres_i = tides_lowres.isel(time=timestep)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(15, 8))\n",
    "ds_lowres_i.plot.imshow(\n",
    "    ax=axes[0],\n",
    "    robust=True,\n",
    "    cmap=\"viridis\",\n",
    "    vmin=ds_i.min().item(),\n",
    "    vmax=ds_i.max().item(),\n",
    ")\n",
    "ds_i.plot.imshow(\n",
    "    ax=axes[1],\n",
    "    robust=True,\n",
    "    cmap=\"viridis\",\n",
    "    vmin=ds_i.min().item(),\n",
    "    vmax=ds_i.max().item(),\n",
    ")\n",
    "for ax in axes:\n",
    "    gridcell_gdf.to_crs(ds.odc.geobox.crs).plot(\n",
    "        ax=ax, facecolor=\"none\", edgecolor=\"black\"\n",
    "    )\n",
    "axes[0].set_title(\"Low resolution (5 x 5 km) modelled tides\")\n",
    "axes[1].set_title(\"Modelled tides reprojected into input satellite grid\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate per-pixel tide cutoffs\n",
    "Based on the entire time-series of tide heights, compute the max and min satellite-observed tide height for each pixel, then calculate tide cutoffs used to restrict our data to satellite observations centred over mid-tide (0 m Above Mean Sea Level)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Determine tide cutoff\n",
    "tide_cutoff_min, tide_cutoff_max = coastlines.raster.tide_cutoffs(ds, tides_lowres, tide_centre=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate yearly composites\n",
    "Export tidally-masked MNDWI median composites for each year, and three-yearly composites used to gapfill poor data coverage areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# If output folder doesn't exist, create it\n",
    "output_dir = f\"data/interim/raster/{raster_version}/{study_area}_{raster_version}\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "\n",
    "# Iterate through each year and export annual and 3-year gapfill composites\n",
    "coastlines.raster.export_annual_gapfill(\n",
    "    ds, output_dir, tide_cutoff_min, tide_cutoff_max, start_year, end_year\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Close Dask client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "***\n",
    "\n",
    "## Additional information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**License:** The code in this notebook is licensed under the [Apache License, Version 2.0](https://www.apache.org/licenses/LICENSE-2.0). \n",
    "Digital Earth Australia data is licensed under the [Creative Commons by Attribution 4.0](https://creativecommons.org/licenses/by/4.0/) license.\n",
    "\n",
    "**Contact:** For assistance with any of the Python code or Jupyter Notebooks in this repository, please post a [Github issue](https://github.com/GeoscienceAustralia/dea-coastlines/issues/new).\n",
    "\n",
    "**Last modified:** November 2022\n",
    "\n",
    "**To cite:**\n",
    "\n",
    "> Bishop-Taylor, R., Nanson, R., Sagar, S., Lymburner, L. (2021). Mapping Australia's dynamic coastline at mean sea level using three decades of Landsat imagery. Remote Sensing of Environment, 267, 112734. Available: https://doi.org/10.1016/j.rse.2021.112734\n",
    ">\n",
    "> Nanson, R., Bishop-Taylor, R., Sagar, S., Lymburner, L., (2022). Geomorphic insights into Australia's coastal change using a national dataset derived from the multi-decadal Landsat archive. Estuarine, Coastal and Shelf Science, 265, p.107712. Available: https://doi.org/10.1016/j.ecss.2021.107712\n",
    ">\n",
    "> Bishop-Taylor, R., Sagar, S., Lymburner, L., Alam, I., Sixsmith, J. (2019). Sub-pixel waterline extraction: characterising accuracy and sensitivity to indices and spectra. Remote Sensing, 11 (24):2984. Available: https://doi.org/10.3390/rs11242984"
   ]
  }
 ],
 "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.8.10"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}