{
"cells": [
{
"cell_type": "markdown",
"id": "0",
"metadata": {},
"source": [
"[](https://demo.leafmap.org/lab/index.html?path=workshops/ICRW_2023.ipynb)\n",
"[](https://colab.research.google.com/github/opengeos/leafmap/blob/master/docs/workshops/ICRW_2023.ipynb)\n",
"\n",
"**An Introduction to Watershed Analysis with Leafmap and WhiteboxTools**\n",
"\n",
"This [notebook](https://leafmap.org/workshops/ICRW_2023) provides an introduction to watershed analysis with Leafmap and WhiteboxTools. It is designed for The Interagency Conference on Research in the Watersheds ([ICRW](https://icrwatersheds.org/)) 2023 workshop - [Working with Geospatial Hydrologic Data for Watershed Analyses in R and Python Using Web Services](https://icrwatersheds.org/icrw8/icrw8-schedule/june-5-workshops/).\n",
"\n",
"- Leafmap: <https://leafmap.org>\n",
"- WhiteboxTools: <https://www.whiteboxgeo.com>\n",
"- WhiteboxTools User Manual: <https://www.whiteboxgeo.com/manual/wbt_book>"
]
},
{
"cell_type": "markdown",
"id": "1",
"metadata": {},
"source": [
"## Installation\n",
"\n",
"Uncomment and run the following cell to install necessary packages for this notebook, including leafmap, geopandas, localtileserver, rio-cogeo, pynhd, py3dep."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2",
"metadata": {},
"outputs": [],
"source": [
"# %pip install leafmap[raster] geopandas"
]
},
{
"cell_type": "markdown",
"id": "3",
"metadata": {},
"source": [
"## Import libraries"
]
},
{
"cell_type": "code",
"execution_count": null,
"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, height=\"600px\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "7",
"metadata": {},
"source": [
"## Add basemaps\n",
"\n",
"Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8",
"metadata": {},
"outputs": [],
"source": [
"m = leafmap.Map()\n",
"m.add_basemap(\"OpenTopoMap\")\n",
"m.add_basemap(\"USGS 3DEP Elevation\")\n",
"m.add_basemap(\"USGS Hydrography\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "9",
"metadata": {},
"source": [
"Add NLCD land cover map and legend."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"metadata": {},
"outputs": [],
"source": [
"m = leafmap.Map(center=[40, -100], zoom=4)\n",
"m.add_basemap(\"HYBRID\")\n",
"m.add_basemap(\"NLCD 2019 CONUS Land Cover\")\n",
"m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "11",
"metadata": {},
"source": [
"Add WMS layers."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"metadata": {},
"outputs": [],
"source": [
"m = leafmap.Map(center=[40, -100], zoom=4)\n",
"m.add_basemap(\"Esri.WorldImagery\")\n",
"url = \"https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?\"\n",
"m.add_wms_layer(\n",
" url,\n",
" layers=\"NLCD_2019_Land_Cover_L48\",\n",
" name=\"NLCD 2019 CONUS Land Cover\",\n",
" format=\"image/png\",\n",
" transparent=True,\n",
")\n",
"m.add_legend(builtin_legend=\"NLCD\", title=\"NLCD Land Cover Type\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "13",
"metadata": {},
"source": [
"## Get watershed data\n",
"\n",
"Let's download watershed data for the Calapooia River basin in Oregon."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "14",
"metadata": {},
"outputs": [],
"source": [
"gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource=\"comid\", simplified=False)"
]
},
{
"cell_type": "markdown",
"id": "15",
"metadata": {},
"source": [
"Plot the watershed boundary on the map."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "16",
"metadata": {},
"outputs": [],
"source": [
"m = leafmap.Map()\n",
"m.add_gdf(gdf, layer_name=\"Catchment\", info_mode=None)\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "17",
"metadata": {},
"source": [
"Save the watershed boundary to a GeoJSON or shapefile."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "18",
"metadata": {},
"outputs": [],
"source": [
"gdf.to_file(\"basin.geojson\", driver=\"GeoJSON\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "19",
"metadata": {},
"outputs": [],
"source": [
"gdf.to_file(\"basin.shp\")"
]
},
{
"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": [
"leafmap.get_3dep_dem(\n",
" gdf, resolution=30, output=\"dem.tif\", 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.tif\", palette=\"terrain\", layer_name=\"DEM\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "24",
"metadata": {},
"source": [
"## Get DEM metadata"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25",
"metadata": {},
"outputs": [],
"source": [
"metadata = leafmap.image_metadata(\"dem.tif\")\n",
"metadata"
]
},
{
"cell_type": "markdown",
"id": "26",
"metadata": {},
"source": [
"Get a summary statistics of the DEM."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "27",
"metadata": {},
"outputs": [],
"source": [
"metadata[\"bands\"]"
]
},
{
"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=\"60\", vmax=1500, 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": null,
"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": "34",
"metadata": {},
"source": [
"Display the WhiteboxTools interface."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "35",
"metadata": {},
"outputs": [],
"source": [
"leafmap.whiteboxgui()"
]
},
{
"cell_type": "markdown",
"id": "36",
"metadata": {},
"source": [
"## Set working directory"
]
},
{
"cell_type": "code",
"execution_count": null,
"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.tif\", \"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_raster(\"smoothed.tif\", palette=\"terrain\", layer_name=\"Smoothed DEM\")\n",
"m.add_geojson(\"basin.geojson\", layer_name=\"Watershed\", info_mode=None)\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": null,
"id": "45",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"hillshade.tif\", layer_name=\"Hillshade\")\n",
"m.layers[-1].opacity = 0.6"
]
},
{
"cell_type": "markdown",
"id": "46",
"metadata": {},
"source": [
"## Find no-flow cells\n",
"\n",
"Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "47",
"metadata": {},
"outputs": [],
"source": [
"wbt.find_no_flow_cells(\"smoothed.tif\", \"noflow.tif\")"
]
},
{
"cell_type": "markdown",
"id": "48",
"metadata": {},
"source": [
"Display the no-flow cells on the map."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"noflow.tif\", layer_name=\"No Flow Cells\")\n",
"m"
]
},
{
"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": "52",
"metadata": {},
"source": [
"Alternatively, you can use depression breaching to fill the depressions."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "53",
"metadata": {},
"outputs": [],
"source": [
"wbt.breach_depressions(\"smoothed.tif\", \"breached.tif\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "54",
"metadata": {},
"outputs": [],
"source": [
"wbt.find_no_flow_cells(\"breached.tif\", \"noflow2.tif\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "55",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"noflow2.tif\", layer_name=\"No Flow Cells after Breaching\")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "56",
"metadata": {},
"source": [
"## Delineate flow direction"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "57",
"metadata": {},
"outputs": [],
"source": [
"wbt.d8_pointer(\"breached.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(\"breached.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\")\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=5000)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"streams.tif\", layer_name=\"Streams\")"
]
},
{
"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": null,
"id": "66",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"distance_to_outlet.tif\", layer_name=\"Distance to Outlet\")"
]
},
{
"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}\n",
")\n",
"m"
]
},
{
"cell_type": "markdown",
"id": "72",
"metadata": {},
"source": [
"## Delineate basins"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "73",
"metadata": {},
"outputs": [],
"source": [
"wbt.basins(\"flow_direction.tif\", \"basins.tif\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "74",
"metadata": {},
"outputs": [],
"source": [
"m.add_raster(\"basins.tif\", layer_name=\"Basins\")"
]
},
{
"cell_type": "markdown",
"id": "75",
"metadata": {},
"source": [
"## Delineate the longest flow path"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "76",
"metadata": {},
"outputs": [],
"source": [
"wbt.longest_flowpath(\n",
" dem=\"breached.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": null,
"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"
]
},
{
"cell_type": "code",
"execution_count": null,
"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 = [-122.613559, 44.284383]\n",
" leafmap.coords_to_vector(coords, output=\"pour_point.shp\", crs=\"EPSG:3857\")"
]
},
{
"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=300\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "84",
"metadata": {},
"outputs": [],
"source": [
"m.add_shp(\"pour_point_snapped.shp\", layer_name=\"Pour Point\")"
]
},
{
"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\")\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": null,
"id": "90",
"metadata": {},
"outputs": [],
"source": [
"m.add_shp(\n",
" \"watershed.shp\",\n",
" layer_name=\"Watershed Vector\",\n",
" style={\"color\": \"#ffff00\", \"weight\": 3},\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}