Drainage networks for wetlands

wetland_drainage.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/gist/giswqs/c4003e5ab541a464c5061fefaff84a37)\n",
    "\n",
    "**Delineating Drainage Networks for Wetlands**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Installation\n",
    "\n",
    "Uncomment and run the following cell to install necessary packages for this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install \"leafmap[raster]\" geopandas py3dep mapclassify"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import leafmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "## Create interactive maps\n",
    "\n",
    "Specify the map center, zoom level, and height."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[40, -100], zoom=4)\n",
    "m.add_basemap(\"HYBRID\")\n",
    "m.add_basemap(\"USGS 3DEP Elevation Index\", opacity=0.5)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Add wetland sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "742b8171",
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://github.com/opengeos/datasets/releases/download/hydrology/NC_wetland_site_data.csv\"\n",
    "df = leafmap.csv_to_df(url)\n",
    "print(f\"Number of stations: {len(df)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "834ab213",
   "metadata": {},
   "source": [
    "Add the gaging stations to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9276d2e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[34.9399, -77.3657], zoom=7)\n",
    "m.add_basemap(\"HYBRID\")\n",
    "m.add_basemap(\"USGS 3DEP Elevation Index\", opacity=0.5)\n",
    "m.add_basemap(\"FWS NWI Wetlands\")\n",
    "m.add_basemap(\"USGS Hydrography\", show=False)\n",
    "m.add_points_from_xy(df, x=\"Longitude\", y=\"Latitude\", layer_name=\"Wetland Sites\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d1f8fa0",
   "metadata": {},
   "source": [
    "You can check out the interactive map at https://apps.opengeos.org/NC_wetland_sites.html."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {},
   "source": [
    "## Get watershed data\n",
    "\n",
    "Use the drawing tool to draw a polygon or place a marker on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c163b409",
   "metadata": {},
   "outputs": [],
   "source": [
    "if m.user_roi is not None:\n",
    "    geometry = m.user_roi_bounds()\n",
    "else:\n",
    "    geometry = (-81.861739,35.580196)\n",
    "    m.set_center(geometry[0], geometry[1], zoom=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3072dde",
   "metadata": {},
   "source": [
    "Download NHD WBD HUC12 data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77879bd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "point_geometry = {\"x\": geometry[0], \"y\": geometry[1]}\n",
    "gdf = leafmap.get_wbd(point_geometry, digit=12, return_geometry=True)\n",
    "huc12_id = gdf[\"huc12\"].values[0]\n",
    "print(f\"HUC12 ID: {huc12_id}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be9cf24",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.explore()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d426afc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_gdf(gdf, layer_name=\"WBD12\", info_mode=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {},
   "source": [
    "## Download DEM\n",
    "\n",
    "Download a digital elevation model (DEM) for the watershed from the USGS 3DEP Elevation service. Convert the DEM to a Cloud Optimized GeoTIFF (COG)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "dem = f\"{huc12_id}.tif\"\n",
    "ds = leafmap.get_3dep_dem(\n",
    "    gdf, resolution=10, output=dem, dst_crs=\"EPSG:3857\", to_cog=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22",
   "metadata": {},
   "source": [
    "Display the DEM on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(dem, palette=\"terrain\", layer_name=\"DEM\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {},
   "source": [
    "## Get DEM metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "metadata = leafmap.image_metadata(dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e59d67ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin, vmax = leafmap.image_min_max(dem)\n",
    "print(f\"Min: {vmin}, Max: {vmax}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {},
   "source": [
    "## Add colorbar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_colormap(cmap=\"terrain\", vmin=int(vmin), vmax=int(vmax), label=\"Elevation (m)\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "## Initialize WhiteboxTools\n",
    "\n",
    "Initialize the WhiteboxTools class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "31",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt = leafmap.WhiteboxTools()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32",
   "metadata": {},
   "source": [
    "Check the WhiteboxTools version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.version()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36",
   "metadata": {},
   "source": [
    "## Set working directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "37",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.set_working_dir(os.getcwd())\n",
    "wbt.verbose = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38",
   "metadata": {},
   "source": [
    "## Smooth DEM\n",
    "\n",
    "All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.feature_preserving_smoothing(dem, \"smoothed.tif\", filter=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40",
   "metadata": {},
   "source": [
    "Display the smoothed DEM and watershed boundary on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_basemap(\"HYBRID\")\n",
    "m.add_basemap(\"FWS NWI Wetlands\")\n",
    "m.add_raster(\"smoothed.tif\", palette=\"terrain\", layer_name=\"Smoothed DEM\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42",
   "metadata": {},
   "source": [
    "## Create hillshade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.hillshade(\"smoothed.tif\", \"hillshade.tif\", azimuth=315, altitude=35)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44",
   "metadata": {},
   "source": [
    "Overlay the hillshade on the smoothed DEM with transparency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "45",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"hillshade.tif\", layer_name=\"Hillshade\")\n",
    "m.layers[-1].opacity = 0.6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50",
   "metadata": {},
   "source": [
    "## Fill depressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.fill_depressions(\"smoothed.tif\", \"filled.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56",
   "metadata": {},
   "source": [
    "## Delineate flow direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.d8_pointer(\"filled.tif\", \"flow_direction.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58",
   "metadata": {},
   "source": [
    "## Calculate flow accumulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.d8_flow_accumulation(\"filled.tif\", \"flow_accum.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"flow_accum.tif\", layer_name=\"Flow Accumulation\", nodata=-32768, visible=False)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61",
   "metadata": {},
   "source": [
    "## Extract streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.extract_streams(\"flow_accum.tif\", \"streams.tif\", threshold=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "63",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"streams.tif\", layer_name=\"Streams\", visible=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64",
   "metadata": {},
   "source": [
    "## Calculate distance to outlet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.distance_to_outlet(\n",
    "    \"flow_direction.tif\", streams=\"streams.tif\", output=\"distance_to_outlet.tif\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "66",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"distance_to_outlet.tif\", layer_name=\"Distance to Outlet\", visible=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67",
   "metadata": {},
   "source": [
    "## Vectorize streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.raster_streams_to_vector(\n",
    "    \"streams.tif\", d8_pntr=\"flow_direction.tif\", output=\"streams.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69",
   "metadata": {},
   "source": [
    "The raster_streams_to_vector tool has a bug. The output vector file is missing the coordinate system. Use leafmap.vector_set_crs() to set the coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.vector_set_crs(source=\"streams.shp\", output=\"streams.shp\", crs=\"EPSG:3857\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"streams.shp\", layer_name=\"Streams Vector\", style={\"color\": \"#ff0000\", \"weight\": 3}, info_mode=None\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75",
   "metadata": {},
   "source": [
    "## Delineate the longest flow path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4607022e",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.basins(\"flow_direction.tif\", \"basins.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.longest_flowpath(\n",
    "    dem=\"filled.tif\", basins=\"basins.tif\", output=\"longest_flowpath.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77",
   "metadata": {},
   "source": [
    "Select only the longest flow path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "78",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.select_largest(\n",
    "    \"longest_flowpath.shp\", column=\"LENGTH\", output=\"longest_flowpath.shp\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"longest_flowpath.shp\",\n",
    "    layer_name=\"Longest Flowpath\",\n",
    "    style={\"color\": \"#ff0000\", \"weight\": 3},\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80",
   "metadata": {},
   "source": [
    "## Generate a pour point\n",
    "\n",
    "Use the drawing tool to place a pour point on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "81",
   "metadata": {},
   "outputs": [],
   "source": [
    "if m.user_roi is not None:\n",
    "    m.save_draw_features(\"pour_point.shp\", crs=\"EPSG:3857\")\n",
    "else:\n",
    "    coords = geometry\n",
    "    leafmap.coords_to_vector(coords, output=\"pour_point.shp\", crs=\"EPSG:3857\")\n",
    "    m.add_marker(location=coords[::-1], name=\"Pour Point\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82",
   "metadata": {},
   "source": [
    "## Snap pour point to stream"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.snap_pour_points(\n",
    "    \"pour_point.shp\", \"flow_accum.tif\", \"pour_point_snapped.shp\", snap_dist=50\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "84",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\"pour_point_snapped.shp\", layer_name=\"Snapped Pour Point\", info_mode=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85",
   "metadata": {},
   "source": [
    "## Delineate watershed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.watershed(\"flow_direction.tif\", \"pour_point_snapped.shp\", \"watershed.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"watershed.tif\", layer_name=\"Watershed\", visible=False)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88",
   "metadata": {},
   "source": [
    "## Convert watershed raster to vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89",
   "metadata": {},
   "outputs": [],
   "source": [
    "wbt.raster_to_vector_polygons(\"watershed.tif\", \"watershed.shp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "90",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_shp(\n",
    "    \"watershed.shp\",\n",
    "    layer_name=\"Watershed Vector\",\n",
    "    style={\"color\": \"#ffff00\", \"weight\": 3},\n",
    "    info_mode=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "a29105d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_basemap(\"USGS Hydrography\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geo",
   "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}