notebooks/DEACoastlines_generation_continental.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# DEA Coastlines continental hotspots <img align=\"right\" src=\"https://github.com/GeoscienceAustralia/dea-notebooks/raw/develop/Supplementary_data/dea_logo.jpg\">\n",
    "\n",
    "This code combines individual datasets into continental DEA Coastlines layers:\n",
    "* Combines output shorelines and rates of change statistics point vectors into single continental datasets\n",
    "* Aggregates this data to produce moving window hotspot datasets that summarise coastal change at regional and continental scale.\n",
    "* Writes outputs to GeoPackage and zipped shapefiles\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",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "Set working directory to top level of repo to ensure links work correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/Robbi/dea-coastlines\n"
     ]
    }
   ],
   "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": {},
   "outputs": [],
   "source": [
    "pip install -r requirements.in --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "%load_ext line_profiler\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geohash as gh\n",
    "import geopandas as gpd\n",
    "from pathlib import Path\n",
    "\n",
    "from coastlines.utils import STYLES_FILE\n",
    "from coastlines.continental import wms_fields\n",
    "from coastlines.vector import points_on_line, change_regress, vector_schema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "vector_version = \"testing2023\"\n",
    "continental_version = \"testing2023\"\n",
    "baseline_year = 2023\n",
    "hotspots_radius = [10000, 5000, 1000]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make output directory and identify files to load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Make output directory \n",
    "output_dir = Path(f\"data/processed/{continental_version}\")\n",
    "output_dir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "# Setup input and output file paths\n",
    "shoreline_paths = (\n",
    "    f\"data/interim/vector/{vector_version}/*/\" f\"annualshorelines*.shp\"\n",
    ")\n",
    "ratesofchange_paths = (\n",
    "    f\"data/interim/vector/{vector_version}/*/\" f\"ratesofchange*.shp\"\n",
    ")\n",
    "\n",
    "# Output path for geopackage and zipped shapefiles\n",
    "OUTPUT_GPKG = output_dir / f\"coastlines_{continental_version}.gpkg\"\n",
    "OUTPUT_SHPS = output_dir / f\"coastlines_{continental_version}.shp.zip\"\n",
    "\n",
    "# If shapefile zip exists, delete it first\n",
    "if OUTPUT_SHPS.exists():\n",
    "    OUTPUT_SHPS.unlink()    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine data\n",
    "### Shorelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system(\n",
    "    f\"ogrmerge.py -o \"\n",
    "    f\"{OUTPUT_GPKG} {shoreline_paths} \"\n",
    "    f\"-single -overwrite_ds -t_srs epsg:3577 \"\n",
    "    f\"-nln shorelines_annual\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rate of change points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system(\n",
    "    f\"ogrmerge.py \"\n",
    "    f\"-o {OUTPUT_GPKG} {ratesofchange_paths} \"\n",
    "    f\"-single -update -t_srs epsg:3577 \"\n",
    "    f\"-nln rates_of_change\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Continental hotspots\n",
    "### Prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Load continental rates of change data\n",
    "ratesofchange_gdf = gpd.read_file(OUTPUT_GPKG, layer=\"rates_of_change\").set_index('uid')\n",
    "\n",
    "# Load continental shorelines data\n",
    "shorelines_gdf = gpd.read_file(OUTPUT_GPKG, layer=\"shorelines_annual\").set_index('year')\n",
    "shorelines_gdf = shorelines_gdf.loc[shorelines_gdf.geometry.is_valid]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Calculate hotspots\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating hotspots at 10000 m\n",
      "Calculating hotspots at 5000 m\n",
      "Calculating hotspots at 1000 m\n"
     ]
    }
   ],
   "source": [
    "# Convert radius to list if not already\n",
    "hotspots_radius = (\n",
    "    [hotspots_radius] if not isinstance(hotspots_radius, list) else hotspots_radius\n",
    ")\n",
    "\n",
    "for i, radius in enumerate(hotspots_radius):\n",
    "\n",
    "    # Extract hotspot points\n",
    "    print(f\"Calculating hotspots at {radius} m\")\n",
    "    hotspots_gdf = points_on_line(\n",
    "        shorelines_gdf,\n",
    "        index=baseline_year,\n",
    "        distance=int(radius / 2),\n",
    "    )\n",
    "\n",
    "    # Create polygon windows by buffering points\n",
    "    buffered_gdf = hotspots_gdf[[\"geometry\"]].copy()\n",
    "    buffered_gdf[\"geometry\"] = buffered_gdf.buffer(radius)\n",
    "\n",
    "    # Spatial join rate of change points to each polygon\n",
    "    hotspot_grouped = (\n",
    "        ratesofchange_gdf.loc[\n",
    "            ratesofchange_gdf.certainty == \"good\",\n",
    "            ratesofchange_gdf.columns.str.contains(\"dist_|geometry\"),\n",
    "        ]\n",
    "        .sjoin(buffered_gdf, predicate=\"within\")\n",
    "        .groupby(\"index_right\")\n",
    "    )\n",
    "\n",
    "    # Aggregate/summarise values by taking median of all points\n",
    "    # within each buffered polygon\n",
    "    hotspot_values = hotspot_grouped.median(numeric_only=True).round(2)\n",
    "\n",
    "    # Extract year from distance columns (remove \"dist_\")\n",
    "    x_years = hotspot_values.columns.str.replace(\"dist_\", \"\").astype(int)\n",
    "\n",
    "    # Compute coastal change rates by linearly regressing annual\n",
    "    # movements vs. time\n",
    "    rate_out = hotspot_values.apply(\n",
    "        lambda row: change_regress(\n",
    "            y_vals=row.values.astype(float), x_vals=x_years, x_labels=x_years\n",
    "        ),\n",
    "        axis=1,\n",
    "    )\n",
    "\n",
    "    # Add rates of change back into dataframe\n",
    "    hotspot_values[\n",
    "        [\"rate_time\", \"incpt_time\", \"sig_time\", \"se_time\", \"outl_time\"]\n",
    "    ] = rate_out\n",
    "\n",
    "    # Join aggregated values back to hotspot points after\n",
    "    # dropping unused columns (regression intercept)\n",
    "    hotspots_gdf = hotspots_gdf.join(hotspot_values.drop(\"incpt_time\", axis=1))\n",
    "\n",
    "    # Add hotspots radius attribute column\n",
    "    hotspots_gdf[\"radius_m\"] = radius\n",
    "\n",
    "    # Initialise certainty column with good values\n",
    "    hotspots_gdf[\"certainty\"] = \"good\"\n",
    "\n",
    "    # Identify any points with insufficient observations and flag these as\n",
    "    # uncertain. We can obtain a sensible threshold by dividing the\n",
    "    # hotspots radius by 30 m along-shore rates of change point distance)\n",
    "    hotspots_gdf[\"n\"] = hotspot_grouped.size()\n",
    "    hotspots_gdf[\"n\"] = hotspots_gdf[\"n\"].fillna(0)\n",
    "    hotspots_gdf.loc[\n",
    "        hotspots_gdf.n < (radius / 30), \"certainty\"\n",
    "    ] = \"insufficient points\"\n",
    "\n",
    "    # Generate a geohash UID for each point and set as index\n",
    "    uids = (\n",
    "        hotspots_gdf.geometry.to_crs(\"EPSG:4326\")\n",
    "        .apply(lambda x: gh.encode(x.y, x.x, precision=11))\n",
    "        .rename(\"uid\")\n",
    "    )\n",
    "    hotspots_gdf = hotspots_gdf.set_index(uids)\n",
    "\n",
    "    # Export hotspots to file, incrementing name for each layer\n",
    "    try:\n",
    "\n",
    "        # Export to geopackage\n",
    "        layer_name = f\"hotspots_zoom_{range(0, 10)[i + 1]}\"\n",
    "        hotspots_gdf.to_file(\n",
    "            OUTPUT_GPKG,\n",
    "            layer=layer_name,\n",
    "            schema={\"properties\": vector_schema(hotspots_gdf), \"geometry\": \"Point\"},\n",
    "        )\n",
    "\n",
    "        # Add additional WMS fields and add to shapefile\n",
    "        hotspots_gdf = pd.concat([hotspots_gdf, wms_fields(gdf=hotspots_gdf)], axis=1)\n",
    "        hotspots_gdf.to_file(\n",
    "            OUTPUT_SHPS,\n",
    "            layer=f\"coastlines_{continental_version}_{layer_name}\",\n",
    "            schema={\"properties\": vector_schema(hotspots_gdf), \"geometry\": \"Point\"},\n",
    "        )\n",
    "\n",
    "    except ValueError as e:\n",
    "\n",
    "        print(f\"Failed to generate hotspots with error: {e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write zipped shapefiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Add rates of change points to shapefile zip\n",
    "# Add additional WMS fields and add to shapefile\n",
    "ratesofchange_gdf = pd.concat(\n",
    "    [ratesofchange_gdf, wms_fields(gdf=ratesofchange_gdf)], axis=1\n",
    ")\n",
    "\n",
    "ratesofchange_gdf.to_file(\n",
    "    OUTPUT_SHPS,\n",
    "    layer=f\"coastlines_{continental_version}_rates_of_change\",\n",
    "    schema={\"properties\": vector_schema(ratesofchange_gdf), \"geometry\": \"Point\"},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Add annual shorelines to shapefile zip\n",
    "shorelines_gdf.to_file(\n",
    "    OUTPUT_SHPS,\n",
    "    layer=f\"coastlines_{continental_version}_shorelines_annual\",\n",
    "    schema={\n",
    "        \"properties\": vector_schema(shorelines_gdf),\n",
    "        \"geometry\": [\"MultiLineString\", \"LineString\"],\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write styles to Geopackage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Insert styles table into GeoPackage\n",
    "styles = gpd.read_file(STYLES_FILE)\n",
    "styles.to_file(OUTPUT_GPKG, layer=\"layer_styles\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.10.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}