WhiteboxTools

17. WhiteboxTools#

Open In Colab

17.1. Installation#

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

# %pip install "leafmap[raster]" geopandas

17.2. Import libraries#

import os
import leafmap

17.3. Create interactive maps#

Specify the map center, zoom level, and height.

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

17.4. Add basemaps#

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

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.

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.

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

17.5. Get watershed data#

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

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

Plot the watershed boundary on the map.

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

Save the watershed boundary to a GeoJSON or shapefile.

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

gdf.to_file("basin.shp")

17.6. 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).

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

Display the DEM on the map.

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

17.7. Get DEM metadata#

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

Get a summary statistics of the DEM.

metadata["bands"]

17.8. Add colorbar#

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

17.9. Initialize WhiteboxTools#

Initialize the WhiteboxTools class.

wbt = leafmap.WhiteboxTools()

Check the WhiteboxTools version.

wbt.version()

Display the WhiteboxTools interface.

leafmap.whiteboxgui()

17.10. Set working directory#

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

17.11. Smooth DEM#

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

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

Display the smoothed DEM and watershed boundary on the map.

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

17.12. Create hillshade#

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

Overlay the hillshade on the smoothed DEM with transparency.

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

17.13. Find no-flow cells#

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

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

Display the no-flow cells on the map.

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

17.14. Fill depressions#

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

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

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

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

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

17.15. Delineate flow direction#

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

17.16. Calculate flow accumulation#

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

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

17.17. Extract streams#

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

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

17.18. Calculate distance to outlet#

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

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

17.19. Vectorize streams#

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.

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

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

17.20. Delineate basins#

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

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

17.21. Delineate the longest flow path#

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

Select only the longest flow path.

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

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

17.22. Generate a pour point#

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")

17.23. Snap pour point to stream#

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

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

17.24. Delineate watershed#

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

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

17.25. Convert watershed raster to vector#

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

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