Interpolating LiDAR data
- How do I convert a LAS point cloud into a raster?
- How do I exclude points with certain classifications?
- I have many LAS files and want to interpolate all of them at once
- What if my data contains anomalously high/low points?
- My data are in LAZ format. How do I interpolate them?
- How do I interpolate an image from the intensity data?
- How do I decide on an appropriate grid resolution?
- My raster contains NoData gaps. How do I remove these?
- How do I combine many LiDAR tiles into a single raster?
- Is there a complete example LiDAR processing workflow available?
How do I convert a LAS point cloud into a raster?
Converting your LiDAR data into a raster requires an interpolation operation. There are many such interpolation methods. The following is an example of how to interpolate the last-return points of a LAS (LiDAR) file using an inverse distance weighted (IDW) interpolation scheme, with a search window radius of 2.5 m, an exponent of 2.0, and an output grid resolution of 1.5 m.
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
wbt.lidar_idw_interpolation(
i="myFile.las",
output="myRaster.tif",
parameter="elevation",
returns="last",
resolution=1.5,
weight=2.0,
radius=2.5
)Other methods for gridding a LAS file include nearest neighbour, Delaunay triangulation (TINing), block minimum, and block maximum gridding schemes.
How do I exclude points with certain classifications?
It is commonly the case that points with certain
class values should be excluded from the gridding of LiDAR data. For
example, you may wish to exclude points associated with vegetation,
buildings, bridges, utility lines, etc. The LidarIdwInterpolation
and LidarNearestNeighbourGridding tools allow for excluded
point classes using the exclude_cls parameter. The parameter
takes a numeric list as input, e.g. exclude_cls='3,4,5,6,7,18'.
Class values follow those of the LAS v.1.4 specifications:
| Classification Value | Meaning |
|---|---|
| 0 | Created, never classified |
| 1 | Unclassified3 |
| 2 | Ground |
| 3 | Low Vegetation |
| 4 | Medium Vegetation |
| 5 | High Vegetation |
| 6 | Building |
| 7 | Low Point (noise) |
| 8 | Reserved |
| 9 | Water |
| 10 | Rail |
| 11 | Road Surface |
| 12 | Reserved |
| 13 | Wire – Guard (Shield) |
| 14 | Wire – Conductor (Phase) |
| 15 | Transmission Tower |
| 16 | Wire-structure Connector (e.g. Insulator) |
| 17 | Bridge Deck |
| 18 | High Noise |
Of course, not all LAS files have had point classifications applied and stored. To determine whether your data contains point class data, you can run the LidarInfo tool before interpolation.
I have many LAS files and want to interpolate all of them at once
When you have hundreds, or even thousands, of LAS files you might be inclined to write a Python script that calls the above function for each input file contained within a folder. But that isn't the best way to handle this common situation. Instead, if the input (i) and output parameters are left unspecified, each of WhiteboxTool's LiDAR gridding methods will interpolate all of the LAS files in the working directory, e.g.
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
wbt.lidar_idw_interpolation(
parameter="elevation",
returns="last",
resolution=1.0,
weight=1.0,
radius=2.5
)Using this approach to folder-based interpolation has some advantages other than a greatly simplified script. WhiteboxTools will be able to parallelize the operation better, greatly improving the overall time required to interpolate the batch of files. Also, the gridding operations will be carried out with a strip of buffered data surrounding each LiDAR tile, i.e. there will be reduced edge-effects. This will reduce the potential for artifacts in the final mosaiced DEM.
What if my data contains anomalously high/low points?
This is a fairly common problem with LiDAR data.
if you're fortunate, these points, which often fall hundreds of meters above
or below the terrain surface, will be classified appropriately. When
this is the case, you may simply exclude the points with class values of
7 (low point) and 18 (high point). Alternatively, you may use the optional
minz and maxz interpolation parameters to
exclude unclassified outlier points. Lastly, you may remove these points
from the original point cloud data set using the LidarRemoveOutliers
tool.
My data are in LAZ format. How do I interpolate them?
WhiteboxTools does not currently support the compressed LiDAR format LAZ. To use these data, you will first need to decompress the files to a LAS format. You may wish to use LasTools for this purpose.
How do I interpolate an image from the intensity data?
The parameter argument of the IDW and nearest neighbour
interpolator tools allows you to interpolate intensity data (options include 'elevation',
'intensity', 'class', 'scan angle', and 'user data'). Here is an example:
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
wbt.lidar_nearest_neighbour_gridding(
"in.las", "out.tif", parameter="intensity")How do I decide on an appropriate grid resolution?
You want to choose a grid resolution where the vast majority of grid cells in the area covered by data have at least one return point. If you are interpolating with last-return points only, then this will necessarily reduce the potential resolution. Ultimately, there is not single appropriate value and the range of suitable resolutions will depend on the distribution of point density with the area of coverage. If the specified resolution is too high given the point density of the LiDAR data set, many of the grid cells will either be NoData holes, or represent interpolated values from relatively distant (up to the search radius) points. A higher than necessary grid resolution will also make working with the final mosaiced DEM raster more challenging, due to the computational effort needed to work with massive rasters and increase the storage and memory requirements. It is advisable to experiment with the LidarPointDensity and LidarPointStats tools before deciding upon a grid resolution for interpolation.
My raster contains NoData gaps. How do I remove these?
First, we need to distinguish between two common areas of NoData values in the interpolated rasters of LiDAR data sets. Because LiDAR data are often collected for irregularly shaped sites, it is frequently the case that LiDAR DEMs have large NoData areas beyond the area of LiDAR point coverage. These are generally acceptable void areas and should not be altered. The more problemmatic void areas are interior data gaps (so called doughnut holes). These generally arise because the point density in an area of LiDAR coverage is lower than the grid resolution (and search radius) dictate in an area. Sometimes these NoData areas are associated with specific non-reflective surfaces, such as water, or areas of dense vegetation (and therefore the last return point density is far lower than in other areas). If the NoData gaps are extensive and spread throughout the area of coverage, that is a sign that you likely need to interpolate either with a coarser grid resolution or a larger search radius, or quite probably both. If your LiDAR DEM has a small number these void areas, and they are not extensive, then you may interpolate to remove the gaps using the FillMissingData tool:
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
wbt.fill_missing_data("dem.tif", "new_dem.tif", filter=11)The choice of a filter size will depend on the extent of the largest interior void area.
How do I combine many LiDAR tiles into a single raster?
Often you have many hundred LAS files, which you've intepolated into an equally large number of raster files. To combine these rasters into a single large DEM, use the Mosaic tool.
from os import listdir
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"
# find all GeoTIFFs in the path
input_files = ""
for f in listdir(wbt.work_dir):
if f.endswith(".tif"):
input_files += ";" + f
input_files = input_files[1:] # strips the first ';'
wbt.mosaic(inputs=input_files, output="big_DEM.tif")Is there a complete example LiDAR processing workflow available?
Yes! The following code is an example of some of the common tasks required in processing large LiDAR datasets of many hundreds of LAS files.
import os
from os import path
from WBT.whitebox_tools import WhiteboxTools
def main():
las_files_dir = "/Users/johnlindsay/Documents/data/test_lidar/"
filtered_las_dir = "/Users/johnlindsay/Documents/data/test_lidar/filtered_LAS/"
raster_data_dir = "/Users/johnlindsay/Documents/data/test_lidar/interpolated_grids/"
wbt = WhiteboxTools()
wbt.work_dir = las_files_dir #set working directory
wbt.verbose = False
if not os.path.exists(filtered_las_dir):
os.makedirs(filtered_las_dir)
# Sometimes, I like to extract all of the LAS tiles that overlap with a
# particular area. For example, I might want to interpolate all the files
# overlapping with a watershed. For this, you can use select_tiles_by_polygon.
# Uncomment the four lines below if you want to do this.
# outdir = "/Users/johnlindsay/Documents/data/LAS_files_in_watershed/"
# polygons = "/Users/johnlindsay/Documents/data/LAS_files_in_watershed/watershed.shp"
# wbt.select_tiles_by_polygon(las_files_dir, outdir, polygons)
# las_files_dir = outdir # this way the analysis below works only on the selected tiles.
##################################################################################
# Filter the ground points in the LAS files using lidar_ground_point_filter tool #
##################################################################################
# This one is the SLOWEST part of the workflow and can be avoided if you are
# confident that you have good point classification data, i.e. that the
# vegetation and building classes have been properly populated.
processed_files = []
num_filtered = 1
flag = True
while flag:
file_names = find_las_files(las_files_dir, processed_files)
if len(file_names) > 0: # and len(processed_files) < 3000:
for i in range(len(file_names)):
in_file = las_files_dir + file_names[i]
out_file = filtered_las_dir + file_names[i].replace(".las", "_filtered.las")
print("Processing LAS {} of {} (total filtered={}) {}".format(i+1, len(file_names), num_filtered, file_names[i]))
processed_files.append(file_names[i])
wbt.lidar_ground_point_filter(in_file, out_file, radius=2.0, slope_threshold=45, height_threshold=1.0)
num_filtered += 1
else:
flag = False
##############################
# Interpolate the LAS files. #
##############################
# If you don't use the above ground point filter the directory below must be
# updated to point to the original LAS files.
wbt.work_dir = filtered_las_dir
# You can use either IDW, nearest neighbour, or TINing (Delaunay triangulation)
# for the gridding step. TINing option is available as of WhiteboxTools v0.11.
wbt.lidar_idw_interpolation(parameter="elevation", returns="all", resolution=2.0, weight=1.0, radius=5.0, exclude_cls='3,4,5,6,7,18')
# wbt.lidar_nearest_neighbour_gridding(returns="last", resolution=1.5, radius=2.5, exclude_cls='3,4,5,6,7,18')
# wbt.lidar_tin_gridding(parameter="elevation", returns="all", resolution=2.0, exclude_cls='3,4,5,6,7,18')
###############################################################
# Now mosaic the tiles; this is done using intermediate steps #
###############################################################
if not os.path.exists(raster_data_dir):
os.makedirs(raster_data_dir)
wbt.work_dir = filtered_las_dir #set working directory
wbt.verbose = False
processed_files = []
num_mosaiced = 1
flag = True
round = 1
while flag:
# This will mosaic a maximum of 250 tiles together; these sub-files
# will subsequently be merged. Mosaicing many hundreds of tiles
# together at one time is otherwise too intensive.
file_names = find_tiff_files(filtered_las_dir, processed_files, 250)
if len(file_names) > 1:
in_files = ""
for i in range(len(file_names)):
if i < len(file_names)-1:
in_files += f"{file_names[i]};"
else:
in_files += f"{file_names[i]}"
processed_files.append(file_names[i])
num_mosaiced += 1
out_file = raster_data_dir + f"mosaic{round}.tif"
wbt.mosaic(inputs=in_files, output=out_file, method="nn")
print(f"Processing mosaic {round}; num. files = {num_mosaiced}")
# now clean up the individual tiles
for i in range(len(file_names)):
os.remove(filtered_las_dir + file_names[i])
else:
flag = False
round += 1
wbt.work_dir = raster_data_dir #set working directory
mosaic_file = raster_data_dir + f"final_mosaic.tif"
file_names = find_mosaic_files(raster_data_dir)
if len(file_names) > 1:
in_files = ""
for i in range(len(file_names)):
if i < len(file_names)-1:
in_files += f"{file_names[i]};"
else:
in_files += f"{file_names[i]}"
num_mosaiced += 1
wbt.mosaic(inputs=in_files, output=mosaic_file, method="nn")
# now clean up the intermediate mosaics
for i in range(len(file_names)):
os.remove(raster_data_dir + file_names[i])
##############################################
# Would you like to fill in the NoData gaps? #
##############################################
dem_nodata_filled = raster_data_dir + f"DEM_gaps_filled.tif"
wbt.fill_missing_data(mosaic_file, dem_nodata_filled, filter=11)
######################################################################
# I usually remove off-terrain objects, like any remaining buildings #
######################################################################
dem_no_otos = raster_data_dir + f"DEM_no_OTOs.tif"
wbt.remove_off_terrain_objects(dem_nodata_filled, dem_no_otos, filter=11, slope=15.0)
#####################################
# Would you like to smooth the DEM? #
#####################################
dem_smoothed = raster_data_dir + f"DEM_smoothed.tif"
wbt.feature_preserving_denoise(dem_no_otos, dem_smoothed, filter=11, norm_diff=8.0)
################################
# Want to fix the depressions? #
################################
dem_breached = raster_data_dir + f"DEM_breached.tif"
# Set the maximum breach depth appropriate for the terrain. You can
# also restrict breaching based on a maximum breach channel length.
wbt.breach_depressions(dem_smoothed, dem_breached, max_depth=5.0)
# because we restricted the use of very deep breach channels, there
# may still be depressions in the DEM. To get rid of these, we can
# perform a subsequent depression filling operation.
dem_filled = raster_data_dir + f"DEM_filled.tif"
wbt.fill_depressions(dem_breached, dem_filled)
####################################################################
# Okay, now we have a good base DEM from which we can extract #
# various land-surface parameters. There are really a large #
# number of these parameters available, but I'll just showcase #
# a few common ones here. See the User Manual for a complete list. #
####################################################################
# slope
slope_file = raster_data_dir + f"slope.tif"
wbt.slope(dem_filled, slope_file)
# plan curvature
plan_curv_file = raster_data_dir + f"plan_curv.tif"
wbt.plan_curvature(dem_filled, plan_curv_file)
# profile curvature; other curvatures are available too.
profile_curv_file = raster_data_dir + f"profile_curv.tif"
wbt.profile_curvature(dem_filled, profile_curv_file)
# hillshade (shaded relief raster)
hillshade_file = raster_data_dir + f"hillshade.tif"
wbt.hillshade(dem_filled, hillshade_file)
# relative topographic position (RTP) index
rtp_file = raster_data_dir + f"relative_topographic_position.tif"
wbt.relative_topographic_position(dem_filled, rtp_file, filterx=11, filtery=11)
# or even better, multiscale topographic position
dev_max_mag = raster_data_dir + f"multiscale_topo_position_mag.tif"
dev_max_scale = raster_data_dir + f"multiscale_topo_position_scale.tif"
wbt.max_elevation_deviation(dem_filled, dev_max_mag, dev_max_scale, min_scale=1, max_scale=100, step=2)
# ruggedness index
ruggedness_index_file = raster_data_dir + f"ruggedness_index.tif"
wbt.ruggedness_index(dem_filled, ruggedness_index_file)
# or even better, multiscale roughness
roughness_mag = raster_data_dir + f"multiscale_roughness_mag.tif"
roughness_scale = raster_data_dir + f"multiscale_roughness_scale.tif"
wbt.multiscale_roughness(dem_filled, roughness_mag, roughness_scale, min_scale=1, max_scale=100, step=2)
# D-infinity flow accumulation
flow_accum_file = raster_data_dir + f"dinf_flow_accum.tif"
wbt.d_inf_flow_accumulation(dem_filled, flow_accum_file, log=True)
# There literally hundreds of other useful parameters that could be
# extracted from our DEM using WhiteboxTools. Take a look at the User Manual.
print("Done!")
def find_las_files(input_dir, processed_files):
files = os.listdir(input_dir)
file_names = []
for f in files:
if f.endswith(".las") and f not in processed_files:
file_names.append(f)
return file_names
def find_tiff_files(input_dir, processed_files, max_num=10):
files = os.listdir(input_dir)
file_names = []
for f in files:
if f.endswith(".tif") and f not in processed_files:
if len(file_names) < max_num:
file_names.append(f)
else:
break
return file_names
def find_mosaic_files(input_dir):
files = os.listdir(input_dir)
file_names = []
for f in files:
if "mosaic" in f and (f.endswith(".tif") or f.endswith(".dep")):
file_names.append(f)
return file_names
main()
