ICRW 2023

An Introduction to Watershed Analysis with Leafmap and WhiteboxTools

This notebook provides an introduction to watershed analysis with Leafmap and WhiteboxTools. It is designed for The Interagency Conference on Research in the Watersheds (ICRW) 2023 workshop - Working with Geospatial Hydrologic Data for Watershed Analyses in R and Python Using Web Services.

• Leafmap: https://leafmap.org
• WhiteboxTools: https://www.whiteboxgeo.com
• WhiteboxTools User Manual: https://www.whiteboxgeo.com/manual/wbt_book

Installation¶

Uncomment and run the following cell to install necessary packages for this notebook, including leafmap, geopandas, localtileserver, rio-cogeo, pynhd, py3dep.

In [ ]:

# %pip install leafmap[raster] geopandas

# %pip install leafmap[raster] geopandas

Import libraries¶

In [ ]:

import os
import leafmap

import os import leafmap

Create interactive maps¶

Specify the map center, zoom level, and height.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4, height="600px")
m

m = leafmap.Map(center=[40, -100], zoom=4, height="600px") m

Add basemaps¶

Add OpenTopoMap, USGS 3DEP Elevation, and USGS Hydrography basemaps.

In [ ]:

m = leafmap.Map()
m.add_basemap("OpenTopoMap")
m.add_basemap("USGS 3DEP Elevation")
m.add_basemap("USGS Hydrography")
m

m = leafmap.Map() m.add_basemap("OpenTopoMap") m.add_basemap("USGS 3DEP Elevation") m.add_basemap("USGS Hydrography") m

Add NLCD land cover map and legend.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4)
m.add_basemap("HYBRID")
m.add_basemap("NLCD 2019 CONUS Land Cover")
m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type")
m

m = leafmap.Map(center=[40, -100], zoom=4) m.add_basemap("HYBRID") m.add_basemap("NLCD 2019 CONUS Land Cover") m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") m

Add WMS layers.

In [ ]:

m = leafmap.Map(center=[40, -100], zoom=4)
m.add_basemap("Esri.WorldImagery")
url = "https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?"
m.add_wms_layer(
    url,
    layers="NLCD_2019_Land_Cover_L48",
    name="NLCD 2019 CONUS Land Cover",
    format="image/png",
    transparent=True,
)
m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type")
m

m = leafmap.Map(center=[40, -100], zoom=4) m.add_basemap("Esri.WorldImagery") url = "https://www.mrlc.gov/geoserver/mrlc_display/NLCD_2019_Land_Cover_L48/wms?" m.add_wms_layer( url, layers="NLCD_2019_Land_Cover_L48", name="NLCD 2019 CONUS Land Cover", format="image/png", transparent=True, ) m.add_legend(builtin_legend="NLCD", title="NLCD Land Cover Type") m

Get watershed data¶

Let's download watershed data for the Calapooia River basin in Oregon.

In [ ]:

gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource="comid", simplified=False)

gdf = leafmap.get_nhd_basins(feature_ids=23763529, fsource="comid", simplified=False)

Plot the watershed boundary on the map.

In [ ]:

m = leafmap.Map()
m.add_gdf(gdf, layer_name="Catchment", info_mode=None)
m

m = leafmap.Map() m.add_gdf(gdf, layer_name="Catchment", info_mode=None) m

Save the watershed boundary to a GeoJSON or shapefile.

In [ ]:

gdf.to_file("basin.geojson", driver="GeoJSON")

gdf.to_file("basin.geojson", driver="GeoJSON")

In [ ]:

gdf.to_file("basin.shp")

gdf.to_file("basin.shp")

Download DEM¶

Download a digital elevation model (DEM) for the watershed from the USGS 3DEP Elevation service. Convert the DEM to a Cloud Optimized GeoTIFF (COG).

In [ ]:

leafmap.get_3dep_dem(
    gdf, resolution=30, output="dem.tif", dst_crs="EPSG:3857", to_cog=True
)

leafmap.get_3dep_dem( gdf, resolution=30, output="dem.tif", dst_crs="EPSG:3857", to_cog=True )

Display the DEM on the map.

In [ ]:

m.add_raster("dem.tif", palette="terrain", layer_name="DEM")
m

m.add_raster("dem.tif", palette="terrain", layer_name="DEM") m

Get DEM metadata¶

In [ ]:

metadata = leafmap.image_metadata("dem.tif")
metadata

metadata = leafmap.image_metadata("dem.tif") metadata

Get a summary statistics of the DEM.

In [ ]:

metadata["bands"]

metadata["bands"]

Add colorbar¶

In [ ]:

m.add_colormap(cmap="terrain", vmin="60", vmax=1500, label="Elevation (m)")
m

m.add_colormap(cmap="terrain", vmin="60", vmax=1500, label="Elevation (m)") m

Initialize WhiteboxTools¶

Initialize the WhiteboxTools class.

In [ ]:

wbt = leafmap.WhiteboxTools()

wbt = leafmap.WhiteboxTools()

Check the WhiteboxTools version.

In [ ]:

wbt.version()

wbt.version()

Display the WhiteboxTools interface.

In [ ]:

leafmap.whiteboxgui()

leafmap.whiteboxgui()

Set working directory¶

In [ ]:

wbt.set_working_dir(os.getcwd())
wbt.verbose = False

wbt.set_working_dir(os.getcwd()) wbt.verbose = False

Smooth DEM¶

All WhiteboxTools functions will return 0 if they are successful, and 1 if they are not.

In [ ]:

wbt.feature_preserving_smoothing("dem.tif", "smoothed.tif", filter=9)

wbt.feature_preserving_smoothing("dem.tif", "smoothed.tif", filter=9)

Display the smoothed DEM and watershed boundary on the map.

In [ ]:

m = leafmap.Map()
m.add_raster("smoothed.tif", palette="terrain", layer_name="Smoothed DEM")
m.add_geojson("basin.geojson", layer_name="Watershed", info_mode=None)
m

m = leafmap.Map() m.add_raster("smoothed.tif", palette="terrain", layer_name="Smoothed DEM") m.add_geojson("basin.geojson", layer_name="Watershed", info_mode=None) m

Create hillshade¶

In [ ]:

wbt.hillshade("smoothed.tif", "hillshade.tif", azimuth=315, altitude=35)

wbt.hillshade("smoothed.tif", "hillshade.tif", azimuth=315, altitude=35)

Overlay the hillshade on the smoothed DEM with transparency.

In [ ]:

m.add_raster("hillshade.tif", layer_name="Hillshade")
m.layers[-1].opacity = 0.6

m.add_raster("hillshade.tif", layer_name="Hillshade") m.layers[-1].opacity = 0.6

Find no-flow cells¶

Find cells with undefined flow, i.e. no valid flow direction, based on the D8 flow direction algorithm

In [ ]:

wbt.find_no_flow_cells("smoothed.tif", "noflow.tif")

wbt.find_no_flow_cells("smoothed.tif", "noflow.tif")

Display the no-flow cells on the map.

In [ ]:

m.add_raster("noflow.tif", layer_name="No Flow Cells")
m

m.add_raster("noflow.tif", layer_name="No Flow Cells") m

Fill depressions¶

In [ ]:

wbt.fill_depressions("smoothed.tif", "filled.tif")

wbt.fill_depressions("smoothed.tif", "filled.tif")

Alternatively, you can use depression breaching to fill the depressions.

In [ ]:

wbt.breach_depressions("smoothed.tif", "breached.tif")

wbt.breach_depressions("smoothed.tif", "breached.tif")

In [ ]:

wbt.find_no_flow_cells("breached.tif", "noflow2.tif")

wbt.find_no_flow_cells("breached.tif", "noflow2.tif")

In [ ]:

m.add_raster("noflow2.tif", layer_name="No Flow Cells after Breaching")
m

m.add_raster("noflow2.tif", layer_name="No Flow Cells after Breaching") m

Delineate flow direction¶

In [ ]:

wbt.d8_pointer("breached.tif", "flow_direction.tif")

wbt.d8_pointer("breached.tif", "flow_direction.tif")

Calculate flow accumulation¶

In [ ]:

wbt.d8_flow_accumulation("breached.tif", "flow_accum.tif")

wbt.d8_flow_accumulation("breached.tif", "flow_accum.tif")

In [ ]:

m.add_raster("flow_accum.tif", layer_name="Flow Accumulation")
m

m.add_raster("flow_accum.tif", layer_name="Flow Accumulation") m

Extract streams¶

In [ ]:

wbt.extract_streams("flow_accum.tif", "streams.tif", threshold=5000)

wbt.extract_streams("flow_accum.tif", "streams.tif", threshold=5000)

In [ ]:

m.add_raster("streams.tif", layer_name="Streams")

m.add_raster("streams.tif", layer_name="Streams")

Calculate distance to outlet¶

In [ ]:

wbt.distance_to_outlet(
    "flow_direction.tif", streams="streams.tif", output="distance_to_outlet.tif"
)

wbt.distance_to_outlet( "flow_direction.tif", streams="streams.tif", output="distance_to_outlet.tif" )

In [ ]:

m.add_raster("distance_to_outlet.tif", layer_name="Distance to Outlet")

m.add_raster("distance_to_outlet.tif", layer_name="Distance to Outlet")

Vectorize streams¶

In [ ]:

wbt.raster_streams_to_vector(
    "streams.tif", d8_pntr="flow_direction.tif", output="streams.shp"
)

wbt.raster_streams_to_vector( "streams.tif", d8_pntr="flow_direction.tif", output="streams.shp" )

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.

In [ ]:

leafmap.vector_set_crs(source="streams.shp", output="streams.shp", crs="EPSG:3857")

leafmap.vector_set_crs(source="streams.shp", output="streams.shp", crs="EPSG:3857")

In [ ]:

m.add_shp(
    "streams.shp", layer_name="Streams Vector", style={"color": "#ff0000", "weight": 3}
)
m

m.add_shp( "streams.shp", layer_name="Streams Vector", style={"color": "#ff0000", "weight": 3} ) m

Delineate basins¶

In [ ]:

wbt.basins("flow_direction.tif", "basins.tif")

wbt.basins("flow_direction.tif", "basins.tif")

In [ ]:

m.add_raster("basins.tif", layer_name="Basins")

m.add_raster("basins.tif", layer_name="Basins")

Delineate the longest flow path¶

In [ ]:

wbt.longest_flowpath(
    dem="breached.tif", basins="basins.tif", output="longest_flowpath.shp"
)

wbt.longest_flowpath( dem="breached.tif", basins="basins.tif", output="longest_flowpath.shp" )

Select only the longest flow path.

In [ ]:

leafmap.select_largest(
    "longest_flowpath.shp", column="LENGTH", output="longest_flowpath.shp"
)

leafmap.select_largest( "longest_flowpath.shp", column="LENGTH", output="longest_flowpath.shp" )

In [ ]:

m.add_shp(
    "longest_flowpath.shp",
    layer_name="Longest Flowpath",
    style={"color": "#ff0000", "weight": 3},
)
m

m.add_shp( "longest_flowpath.shp", layer_name="Longest Flowpath", style={"color": "#ff0000", "weight": 3}, ) m

Generate a pour point¶

In [ ]:

if m.user_roi is not None:
    m.save_draw_features("pour_point.shp", crs="EPSG:3857")
else:
    coords = [-122.613559, 44.284383]
    leafmap.coords_to_vector(coords, output="pour_point.shp", crs="EPSG:3857")

if m.user_roi is not None: m.save_draw_features("pour_point.shp", crs="EPSG:3857") else: coords = [-122.613559, 44.284383] leafmap.coords_to_vector(coords, output="pour_point.shp", crs="EPSG:3857")

Snap pour point to stream¶

In [ ]:

wbt.snap_pour_points(
    "pour_point.shp", "flow_accum.tif", "pour_point_snapped.shp", snap_dist=300
)

wbt.snap_pour_points( "pour_point.shp", "flow_accum.tif", "pour_point_snapped.shp", snap_dist=300 )

In [ ]:

m.add_shp("pour_point_snapped.shp", layer_name="Pour Point")

m.add_shp("pour_point_snapped.shp", layer_name="Pour Point")

Delineate watershed¶

In [ ]:

wbt.watershed("flow_direction.tif", "pour_point_snapped.shp", "watershed.tif")

wbt.watershed("flow_direction.tif", "pour_point_snapped.shp", "watershed.tif")

In [ ]:

m.add_raster("watershed.tif", layer_name="Watershed")
m

m.add_raster("watershed.tif", layer_name="Watershed") m

Convert watershed raster to vector¶

In [ ]:

wbt.raster_to_vector_polygons("watershed.tif", "watershed.shp")

wbt.raster_to_vector_polygons("watershed.tif", "watershed.shp")

In [ ]:

m.add_shp(
    "watershed.shp",
    layer_name="Watershed Vector",
    style={"color": "#ffff00", "weight": 3},
)

m.add_shp( "watershed.shp", layer_name="Watershed Vector", style={"color": "#ffff00", "weight": 3}, )