This tutorial is part of a series that includes:
This tutorial is intended to demonstrate how to use the WbW geospatial analysis library to map building footprints from a LiDAR point cloud. It requires a valid license for WbW-Pro. For more information about the WbW library you may refer to the user manual. You may download a copy of the raw Jupyter Notebooks file (*.ipynb) from here.
If you haven't already done so, install the Whitebox Worfklows for Python pip package.
pip install whitebox-workflows
We need to import the whitebox_workflows library into our script and set up our floating license ID, which will be used by WbW-Pro. Once you register your WbW-Pro license, you will be emailed your unique floating license ID, which will likely be a randomly generated three-word phrase, involving an animal. The license below was used for the Geomorphometry 2023 conference in Iasi, Romania and will be valid until mid-August 2023. After that point, you will need to purchase a license for WbW-Pro (about 350USD) to use the scripts below. Licenses can be purchased from Whitebox Geospatial Inc.
import whitebox_workflows
license_id = 'geomorphometry-2023' # Update this value with your own license ID when this one expires.
wbe = whitebox_workflows.WbEnvironment(license_id)
wbe.verbose = True # Let each of the function calls output to stdout.
print(wbe.version()) # Let's see what version of WbW we're working with
Now let's set the script parameters.
# Script parameters
resolution = 0.5 # in meters; determines cell size of DEM/DTM
filter_size = 151 # In grid cells; will need to adjust filter size for largest building
slope_threshold = 15.0 # In degrees; 15 works well but may have to adjust if applied on steeper terrain
min_height1 = 1.5 # affects the definition of edges of features
min_height2 = 3.0 # will need to set this for minimum building height
min_area = 200 # in grid cells; at 200 grid cells, with a resoltion of 0.5 m, it would mean a building has to be at least 50 m^2.
smoothing_factor = 5 # size of smoothing filter; must be odd integer, higher applies more smoothing and set to zero for none
building_footprint_filename = 'building_footprints.shp'
Download the Kitchener_lidar sample LiDAR file. Notice that it may take a few minutes to download the data. In the event that the download takes more than a few minutes, the connection may timeout and you will receive an error. If this should happen, you may download the dataset directly from here but you will need to update the wbe.working_directory to your download folder.
# Download the Kitchener_lidar sample lidar tile
wbe.working_directory = whitebox_workflows.download_sample_data('Kitchener_lidar')
print(f'Data have been stored in: {wbe.working_directory}')
lidar = wbe.read_lidar('Kitchener_lidar.laz') # read in the lidar data set
Finally, let's map the building footprints...
# Interpolate a last-return DEM
dem = wbe.lidar_tin_gridding(lidar, returns_included='last', cell_size=resolution, excluded_classes=[3,4,5])
wbe.write_raster(dem, 'DEM.tif')
# Remove the off-terrain objects (OTOs)
dtm = wbe.remove_off_terrain_objects(dem, filter_size=filter_size, slope_threshold=slope_threshold)
oto_heights = dem - dtm # measure OTO height as a DEM of diff
# wbe.write_raster(oto_heights, 'oto_heights.tif') # uncomment for quality control
# Filter out features based on height and area
otos = oto_heights > min_height1
otos = wbe.clump(otos, zero_background=True)
otos_max_hgt, tmp = wbe.zonal_statistics(oto_heights, otos, stat_type='maximum')
otos = otos_max_hgt > min_height2
# wbe.write_raster(otos, 'otos.tif') # uncomment for quality control
otos = wbe.generalize_classified_raster(raster=otos, area_threshold=min_area, method = "largest")
# wbe.write_raster(otos, 'otos2.tif') # uncomment for quality control
building_footprints = wbe.raster_to_vector_polygons(otos)
if smoothing_factor > 0:
building_footprints = wbe.smooth_vectors(building_footprints, filter_size=smoothing_factor)
# Save the final map
wbe.write_vector(building_footprints, building_footprint_filename)
print('Done!')
You should end up with something that looks a little like this:
Here, we have overlaid the final building_footprints.shp as black lines overtop the DEM.tif output.
Note that the reason that this process requires a license for WbW-Pro, rather than a less expensive WbW license, is the use of the generalize_classified_raster function, which lives in WbW-Pro. This function is essential for simplifying footprint shapes (including removing 'donut holes') and removing small features.
Don't forget to check your license in after you're done using it. If you skip this step, your checked-out license won't be returned to the license pool for others to use for one hour. Be kind to others!
print(wbe.check_in_license(license_id))