Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Geoprocessing with WhiteboxTools

Introduction

Learning Objectives

Why Whitebox?

What is Whitebox?

What can Whitebox do?

How is Whitebox different?

Useful Resources for Whitebox

Installing Whitebox

# %pip install whitebox "pygis[lidar]"
import os
import leafmap
import numpy as np
os.environ["NODATA"] = "-32768"

Watershed Analysis

Create Interactive Maps

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

Download Watershed Data

lat = 44.361169
lon = -122.821802

m = leafmap.Map(center=[lat, lon], zoom=10)
m.add_marker([lat, lon])
m
geometry = {"x": lon, "y": lat}
gdf = leafmap.get_wbd(geometry, digit=10, return_geometry=True)
gdf.explore()
gdf.to_file("basin.geojson")

Download and Display DEM

array = leafmap.get_3dep_dem(
    gdf,
    resolution=30,
    output="dem.tif",
    dst_crs="EPSG:3857",
    to_cog=True,
    overwrite=True,
)
array
m.add_raster("dem.tif", palette="terrain", nodata=np.nan, layer_name="DEM")
m

Get DEM metadata

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

Add colorbar

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

Initialize WhiteboxTools

wbt = leafmap.WhiteboxTools()
wbt.version()
leafmap.whiteboxgui()

Initialize WhiteboxTools

Set working directory

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

Smooth DEM

wbt.feature_preserving_smoothing("dem.tif", "smoothed.tif", filter=9)
m = leafmap.Map()
m.add_basemap("Satellite")
m.add_raster("smoothed.tif", colormap="terrain", layer_name="Smoothed DEM")
m.add_geojson("basin.geojson", layer_name="Watershed", info_mode=None)
m.add_basemap("USGS Hydrography", show=False)
m

Create hillshade

wbt.hillshade("smoothed.tif", "hillshade.tif", azimuth=315, altitude=35)
m.add_raster("hillshade.tif", layer_name="Hillshade")
m.layers[-1].opacity = 0.6

Find no-flow cells

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

Fill depressions

wbt.fill_depressions("smoothed.tif", "filled.tif")
wbt.breach_depressions("smoothed.tif", "breached.tif")
wbt.find_no_flow_cells("breached.tif", "noflow2.tif")
m.layers[-1].visible = False
m.add_raster("noflow2.tif", layer_name="No Flow Cells after Breaching")
m

Delineate flow direction

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

Calculate flow accumulation

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

Extract streams

wbt.extract_streams("flow_accum.tif", "streams.tif", threshold=5000)
m.layers[-1].visible = False
m.add_raster("streams.tif", layer_name="Streams")
m

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

Vectorize streams

wbt.raster_streams_to_vector(
    "streams.tif", d8_pntr="flow_direction.tif", output="streams.shp"
)
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},
    info_mode=None,
)
m

Delineate the longest flow path

wbt.basins("flow_direction.tif", "basins.tif")
wbt.longest_flowpath(
    dem="breached.tif", basins="basins.tif", output="longest_flowpath.shp"
)
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

Generate a pour point

if m.user_roi is not None:
    m.save_draw_features("pour_point.shp", crs="EPSG:3857")
else:
    lat = 44.284642
    lon = -122.611217
    leafmap.coords_to_vector([lon, lat], output="pour_point.shp", crs="EPSG:3857")
    m.add_marker([lat, lon])

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", info_mode=False)

Delineate watershed

wbt.watershed("flow_direction.tif", "pour_point_snapped.shp", "watershed.tif")
m.add_raster("watershed.tif", layer_name="Watershed")
m

Convert watershed raster to vector

wbt.raster_to_vector_polygons("watershed.tif", "watershed.shp")
m.layers[-1].visible = False
m.add_shp(
    "watershed.shp",
    layer_name="Watershed Vector",
    style={"color": "#ffff00", "weight": 3},
    info_mode=False,
)

LiDAR Data Analysis

Set up whitebox

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

Download a sample dataset

url = "https://github.com/opengeos/datasets/releases/download/lidar/madison.zip"
filename = "madison.las"
leafmap.download_file(url, "madison.zip", quiet=True)

Read LAS/LAZ data

laz = leafmap.read_lidar(filename)
laz
str(laz.header.version)

Upgrade file version

las = leafmap.convert_lidar(laz, file_version="1.4")
str(las.header.version)

Write LAS data

leafmap.write_lidar(las, "madison.las")

Histogram analysis

wbt.lidar_histogram("madison.las", "histogram.html")

Visualize LiDAR data

leafmap.view_lidar("madison.las")

Remove outliers

wbt.lidar_elevation_slice("madison.las", "madison_rm.las", minz=0, maxz=450)

Visualize LiDAR data after removing outliers

leafmap.view_lidar("madison_rm.las", cmap="terrain")

Create DSM

wbt.lidar_digital_surface_model(
    "madison_rm.las", "dsm.tif", resolution=1.0, minz=0, maxz=450
)
leafmap.add_crs("dsm.tif", epsg=2255)

Visualize DSM

m = leafmap.Map()
m.add_basemap("Satellite")
m.add_raster("dsm.tif", colormap="terrain", layer_name="DSM")
m

Create DEM

wbt.remove_off_terrain_objects("dsm.tif", "dem.tif", filter=25, slope=15.0)

Visualize DEM

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

Create CHM

chm = wbt.subtract("dsm.tif", "dem.tif", "chm.tif")
m.add_raster("chm.tif", palette="gist_earth", layer_name="CHM")
m.add_layer_manager()
m

Key Takeaways

Exercises

Exercise 1: Watershed Analysis for a Different Location

Exercise 2: LiDAR Analysis for Forest Structure Assessment

Geospatial Programming with Python
Interactive Visualization with Leafmap
Geospatial Programming with Python
3D Mapping with MapLibre