21. Geoprocessing with WhiteboxTools#
21.1. Introduction#
21.2. Learning Objectives#
21.3. Why Whitebox?#
21.3.1. What is Whitebox?#
21.3.2. What can Whitebox do?#
21.3.3. How is Whitebox different?#
21.4. Useful Resources for Whitebox#
21.5. Installing Whitebox#
# %pip install whitebox "pygis[lidar]"
import os
import leafmap
import numpy as np
os.environ["NODATA"] = "-32768"
21.6. Watershed Analysis#
21.6.1. Create Interactive Maps#
m = leafmap.Map()
m.add_basemap("OpenTopoMap")
m.add_basemap("USGS 3DEP Elevation")
m.add_basemap("USGS Hydrography")
m
21.6.2. 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")
21.6.3. 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
21.6.4. Get DEM metadata#
metadata = leafmap.image_metadata("dem.tif")
metadata
21.6.5. Add colorbar#
leafmap.image_min_max("dem.tif")
m.add_colormap(cmap="terrain", vmin="60", vmax=1500, label="Elevation (m)")
m
21.6.6. Initialize WhiteboxTools#
wbt = leafmap.WhiteboxTools()
wbt.version()
leafmap.whiteboxgui()
21.6.7. Initialize WhiteboxTools#
21.6.8. Set working directory#
wbt.set_working_dir(os.getcwd())
wbt.verbose = False
21.6.9. 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
21.6.10. 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
21.6.11. Find no-flow cells#
wbt.find_no_flow_cells("smoothed.tif", "noflow.tif")
m.add_raster("noflow.tif", layer_name="No Flow Cells")
21.6.12. 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
21.6.13. Delineate flow direction#
wbt.d8_pointer("breached.tif", "flow_direction.tif")
21.6.14. Calculate flow accumulation#
wbt.d8_flow_accumulation("breached.tif", "flow_accum.tif")
m.add_raster("flow_accum.tif", layer_name="Flow Accumulation")
21.6.15. 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
21.6.16. 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")
21.6.17. 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
21.6.18. 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
21.6.19. 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])
21.6.20. 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)
21.6.21. Delineate watershed#
wbt.watershed("flow_direction.tif", "pour_point_snapped.shp", "watershed.tif")
m.add_raster("watershed.tif", layer_name="Watershed")
m
21.6.22. 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,
)
21.7. LiDAR Data Analysis#
21.7.1. Set up whitebox#
wbt = leafmap.WhiteboxTools()
wbt.set_working_dir(os.getcwd())
wbt.verbose = False
21.7.2. 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)
21.7.3. Read LAS/LAZ data#
laz = leafmap.read_lidar(filename)
laz
str(laz.header.version)
21.7.4. Upgrade file version#
las = leafmap.convert_lidar(laz, file_version="1.4")
str(las.header.version)
21.7.5. Write LAS data#
leafmap.write_lidar(las, "madison.las")
21.7.6. Histogram analysis#
wbt.lidar_histogram("madison.las", "histogram.html")
21.7.7. Visualize LiDAR data#
leafmap.view_lidar("madison.las")
21.7.8. Remove outliers#
wbt.lidar_elevation_slice("madison.las", "madison_rm.las", minz=0, maxz=450)
21.7.9. Visualize LiDAR data after removing outliers#
leafmap.view_lidar("madison_rm.las", cmap="terrain")
21.7.10. 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)
21.7.11. Visualize DSM#
m = leafmap.Map()
m.add_basemap("Satellite")
m.add_raster("dsm.tif", colormap="terrain", layer_name="DSM")
m
21.7.12. Create DEM#
wbt.remove_off_terrain_objects("dsm.tif", "dem.tif", filter=25, slope=15.0)
21.7.13. Visualize DEM#
m.add_raster("dem.tif", palette="terrain", layer_name="DEM")
m
21.7.14. 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