Fixing spikes in DSM from LiDAR?

Question

I have a DSM derived from drone LiDAR that I want to use as a base to drape imagery and, paired with a DTM, as a basis for estimating canopy height/volume.

Unfortunately, the DSM has has many "spikes" along boundaries between flight lines. I inverted the DSM by multiplying by -1, making the spikes sinks. I determined that the sinks >50 feet deep were errors that I want to fill while sinks <50 feet deep were trees/shrubs that I want to retain. But the FILL tool only allows setting a z limit for the maximum elevation difference between a sink and its pour point to be filled.

I want to set the minimum elevation difference to be filled (e.g. > 50 feet).

Answer 1

I've given an answer below that is based on the WhiteboxTools geoprocessing library (for which I am the developer), but I would guess that you could use a similar workflow in just about any GIS. The only fairly unusual tool is the 'Depth-In-Sink' tool.

from whitebox_tools import WhiteboxTools

wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"

in_dem = "DEM.tif"

# Invert the DEM
wbt.multiply(in_dem, -1, "iDEM.tif")

# Find the depressions, giving each one a unique ID
wbt.sink("iDEM.tif", "depressions.tif")

# Calcuate the depth-in-sink
wbt.depth_in_sink("iDEM.tif", "depth_in_sink.tif")

# Assign each depression its maximum depth-in-sink value
wbt.extract_raster_statistics("depth_in_sink.tif", "depressions.tif", "max_dis.tif", stat="maximum")

# Threshold max_dis and create Booleans for deep deps vs everything else
wbt.greater_than("max_dis.tif", 50.0, "deep_depressions.tif", incl_equals=True)
wbt.less_than("max_dis.tif", 50.0, "not_deep_depressions.tif")

# Now fill the whole inverted DEM
wbt.fill_depressions("iDEM.tif", "iDEM_filled.tif", fix_flats=False)

# Multiply iDEM_filled by deep_depressions
wbt.multiply("iDEM_filled.tif", "deep_depressions.tif", "temp1.tif")

# Multiply iDEM by not_deep_depressions
wbt.multiply("iDEM.tif", "not_deep_depressions.tif", "temp2.tif")

# Add the two temp rasters together
wbt.add("temp1.tif", "temp2.tif", "temp3.tif")

# Invert the result back into the upward direction
wbt.multiply("temp3.tif", -1, "DEM_spikes_removed.tif")

If the features are true spikes (i.e. single grid cells with anomalously high values compared with neighbours), as is often the case with LiDAR DEMs, then you might consider using an image filtering approach instead of the depression-filling method. Filters that might be of use for this purpose (also available in WhiteboxTools) include the AdaptiveFilter, ConservativeSmoothingFilter, LeeFilter, MedianFilter, or the OlympicFilter. If your DEM is quite noisy, you might consider applying the FeaturePreservingDenoise tool. All of these tools may help to improve your DEM with a single operation compared to the above approach. If however your spikes are 'fat' (multiple cells in size) then the depression-filling approach may be best.

Lastly, you might consider instead filtering these points out of the LiDAR point cloud from which the DSM was created. The LidarRemoveOutliers tool will remove LiDAR points that exceed a specified height difference from the average elevation within a neighbourhood. This could be a very effective approach.