03_Shoreline_Derivation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f43c3a3b-26a4-4bc9-962a-25b9cbffbe50",
   "metadata": {},
   "source": [
    "# Shoreline Derivation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edcb2fbb-2274-4892-96f1-d3830bc9a3f0",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "_**Caitlin Haedrich and Pratikshya Regmi**, North Carolina State University_\n",
    "\n",
    "\n",
    "A shoreline is the interface between the land and the ocean. The location of the shoreline is fundamentally important to various aspects of coastal science and engineering. However, the cross-shore and vertical location of the shoreline changes continuously. Thus, differing definitions of the shoreline exist for various applications. Shoreline location is often approximated by a constant elevation contour drawn in a tidal datum. The specific tidal datum and contour elevation may depend on the application. The following figure shows the shoreline migration range in the town of Nags Head, North Carolina (image from [Hardin et al. (2014)](https://link.springer.com/book/10.1007/978-1-4939-1835-5).\n",
    "\n",
    "![](img/shoreline.png)\n",
    "\n",
    "In this notebook we will:\n",
    "* [Derive shorelines from DEMs](#2.-Deriving-and-Smoothing-Shorelines)\n",
    "* [Visualize the extracted shorelines](#3.-Visualization)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e83f37d1-2a56-48fc-a59b-91339c812f57",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "## 1. Import Python Packages and Start GRASS Session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae316400-3f54-4434-8309-35a8def11dfd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True, shell=False).strip()\n",
    ")\n",
    "\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "gj.init(\"./nags_head/PERMANENT\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3081240d-5c25-4f2b-beec-9e3f4cf74ed3",
   "metadata": {},
   "source": [
    "Before we get started, we check the computational region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "860deea3-9392-4b6c-a5c7-46b9f9de7cfa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!g.region -p region=jockeys_ridge"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "168c679f-3863-4ad0-bbd0-38fd2676f70f",
   "metadata": {},
   "source": [
    "We can double check the extent by using the `use_region=True` keyword in `gj.Map`. This will only areas that are within the current computational extent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2aa89c9-5de1-413e-9b2f-b118072b1e73",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img= gj.Map(use_region=True)\n",
    "image=img.d_rgb(red=\"naip_2020.1\", green=\"naip_2020.2\", blue=\"naip_2020.3\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7463c58-54d8-4831-97c4-83de878097b7",
   "metadata": {},
   "source": [
    "**Try it Yourself!**\n",
    "\n",
    "_Modify and save the computational region so it is just the eastern shore. This will be useful because we're going to focus on shoreline changes on the seaward side of the barrier island._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "401ddf28-f4bc-472a-827f-ba4ab4675ac5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "758e9c75-6a69-4236-8ab8-592c08ca45dd",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>👉 <b>click to example</b></summary>\n",
    "    \n",
    "```python\n",
    "!g.region --help\n",
    "!g.region region=jockeys_ridge w=w+1200 -p\n",
    "```\n",
    "</details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c29780cc-589b-441a-b2df-adaa881ea43d",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "## 2. Deriving and Smoothing Shorelines"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5e18edb-9eae-46d4-b2c4-7a13bae3029e",
   "metadata": {},
   "source": [
    "In this next part, we use [`r.contour`](https://grass.osgeo.org/grass83/manuals/r.contour.html) to extract the shoreline contour and [`v.generalize`](https://grass.osgeo.org/grass83/manuals/v.generalize.html) to smooth and clean the resulting line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54fc7f90-1ee4-49fe-abdb-3ff8f0b10b0d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shoreline_elev = 0.8 # shoreline proxy elevation in meters\n",
    "min_points = 400 # minimum number of points for a contour line (threshold to exclude small anomalies)\n",
    "\n",
    "gs.run_command('r.contour', input=\"JR_2014\", output=\"JR_2014_contour\", levels=shoreline_elev, cut=min_points)\n",
    "gs.run_command('v.generalize',input=\"JR_2014_contour\", output=\"JR_2014_shoreline\", method='sliding_averaging', threshold=0, look_ahead=51)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e80ddbaf-1fb5-475e-b321-ca4aa2d498ee",
   "metadata": {},
   "source": [
    "**Try it Yourself!**\n",
    "\n",
    "_Modify the parameters to see how the resulting shoreline changes._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00e2e9f2-cc77-41ae-84a5-37b3a90547d1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img= gj.Map(use_region=True)\n",
    "img.d_rgb(red=\"naip_2014.1\", green=\"naip_2014.2\", blue=\"naip_2014.3\")\n",
    "img.d_vect(map=\"JR_2014_shoreline\", color=\"red\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "828c184e-8532-4b2a-b382-985b4d91f143",
   "metadata": {},
   "source": [
    "### Batch Shoreline Derivation\n",
    "\n",
    "Let's use a for-loop to automate the process and repeat it for all of our elevation rasters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d00f2430-1568-4a5a-a556-16a72c8eee1c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "elev_rasters = gs.read_command(\"g.list\", type=\"raster\", pattern=\"JR_????\").split()\n",
    "elev_rasters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b06ad16e-98ce-4c55-94f1-89f32a6045c4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shoreline_elev = 0.8 # shoreline proxy elevation\n",
    "min_points = 400 # threshold to exclude small anomalies\n",
    "\n",
    "for raster in elev_rasters:\n",
    "    # Set output names\n",
    "    contour = f'{raster}_contour'\n",
    "    shoreline = f'{raster}_shoreline'\n",
    "    # derive shorelines\n",
    "    gs.run_command('r.contour', input=raster, output=contour, levels=shoreline_elev, cut=min_points)\n",
    "    gs.run_command('v.generalize',input=contour, output=shoreline, method='sliding_averaging', threshold=0, look_ahead=51)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe090e84-2294-4e6b-a7bf-895f16d458ff",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "## 3. Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8189679c-ecd0-4f23-bf9a-2906781b3255",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shoreline_series = gs.read_command(\"g.list\", type=\"vector\", pattern=\"*shoreline\").split()\n",
    "shoreline_series"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23c85cb1-968c-4c02-860e-21a000770f98",
   "metadata": {},
   "source": [
    "Here's a quick little bit of code that makes unique colors for each shoreline so we can see change over time better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99c8b1e0-9514-4a3c-8ff6-abbcfa7812b3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n = len(shoreline_series)\n",
    "colors = []\n",
    "\n",
    "for i in range(n):\n",
    "    # Generate shades of blue\n",
    "    blue = int((i / (n)) * 255)\n",
    "    # Convert RGB to hexadecimal format\n",
    "    color = \"#{:02x}{:02x}{:02x}\".format(0, 0, blue)  # Pure blue with varying intensity\n",
    "    colors.append(color)\n",
    "\n",
    "print(colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a51efd39-db4d-47d8-b75f-dd190cf28c18",
   "metadata": {},
   "source": [
    "Let's also create a composite rgb image from our NAIP rgb channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cf0e7fe-ea9a-42c3-8ab2-2c2dfbe92d63",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!r.composite blue=naip_2020.3 green=naip_2020.2 red=naip_2020.1 output=naip_rgb"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc09da4b-0312-4064-9fc8-6dc8ba7eeb12",
   "metadata": {},
   "source": [
    "With a list of our shoreline vectors, hex colors and our NAIP image, we can make an interactive map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2461ed20-42e0-4f45-8878-6e6a5e489b33",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "img=gj.InteractiveMap(width=800, tiles=\"OpenStreetMap\")\n",
    "img.add_raster(\"naip_rgb\")\n",
    "for i, shoreline in enumerate(shoreline_series):\n",
    "    img.add_vector(shoreline, color=colors[i])\n",
    "img.add_layer_control()\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "667c82cd-63c9-4474-818e-ff4bdf967923",
   "metadata": {},
   "source": [
    "## Take it further\n",
    "\n",
    "What to get more advanced? Check out the [additional tutorials](./Additional_Tutorials/Dune_Toe_and_Ridge.ipynb) where we identify the foredune toe and ridge!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a27a27d-f393-4040-87eb-bb52d69266b8",
   "metadata": {},
   "source": [
    "## Further Reading\n",
    "\n",
    "[Paris P.J., Hardin E., Mitasova H., Starek M.J., Kurum O.M., Overton M.F., 2013, Lines in the Sand: Geomorphic and Geospatial Characterization and Interpretation of Sandy Shorelines and Beaches, Geography Compass 7(5), p. 315-343.](https://compass.onlinelibrary.wiley.com/doi/abs/10.1111/gec3.12041)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e7217ba-1df1-4b92-82fe-a93dbcf6fbff",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "GRASS GIS",
   "language": "python",
   "name": "grass"
  },
  "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}