DSM to DEM Conversion

Often you may have some Digital Surface Model (DSM) acquired from a LiDAR data, and you may desire to get rid of buildings and trees in order to get a DEM. There can be found some articles on a subject that would provide more accurate result then that described here. But I will show here my quick-n-dirty way when one don't need extreme accuracy. It is meant to be used for a relatively flat urban areas.

There will be 3 steps:

1. Subtraction of the elevated objects like buildings and trees from DSM.
2. Filling the gaps of the raster cleaned of the elevated objects.
3. Polishing the DEM.

FOSS software you will need: SAGA GIS, QGIS.

Initial DSM (LiDAR)

Subtraction of the elevated objects like buildings and trees from DSM.

In SAGA GIS go Modules -> Grid-Filter -> DTM filter (slope based). This tool will divide your DSM in two rasters: Bare Earth and Removed Objects. You will have to set some options:

• Search radius: this will determine the range of pixels that are going to be processed when the very steep slope will be encountered. A value to small will leave central parts of the wide buildings in Bare Earth raster while value too big will reap almost all pixels to Removed Objects. Try to set a value that will subtract buildings from the raster completely even if a lot of actual surface will be gone with them. A value about 20 would do for a 2 meters resolution DSM.
• Approx. Terrain Slope: if you have a clue about general slope of the area - set it. You may just play with it to look for better result.
• Use Confidence interval: check this box, otherwise over 95% of your DSM will end up in Removed Objects.

Once the calculation completed examine the final rasters. If you will see that there are enough pixels left in Bare Earth (and there are no parts of buildings left) you may proceed with this raster to the next step. 

If too many actual surface pixels passed to Removed Objects - you should do some more additional work. Assess the the highest feasible value for a an actual surface caught in Removed Objects (in SAGA you may set in map layer options to show the cell values above the pixels). Then go Modules -> Grid-Tools -> Change Grid Values and set all the pixels that are below the assessed threshold to 0. Then load this grid into QGIS.

In QGIS you will need a RasterCalc plugin. Use it to assign NULL data pixels of a raster a value of 0 with the following statement: eq([your_raster], -9999, 0). I assume that -9999 is the value assigned to the NULL data pixels. Now subtract this raster from the original LiDAR raster and you will get the raster where all the elevated objects like buildings and trees will have a value of 0  (If your original data contains pixels with values of 0 you will need to create its mask beforehand and then change 0 values outside the mask to a different number). Now go to Raster -> Projection -> Wrap, chose resulted layer to process (and set it as an output), check the "No data values" box and set its value to 0 (or an appropriate value if you had 0 elevations in original LiDAR data and changed 0 values of subtracted raster outside the mask). Now you are finally ready to proceed further.

Removed elevated objects

Filling the gaps of the raster cleaned of the elevated objects

In order to get the surface values of the areas that were under elevated objects we will use SAGA. Modules -> Grid-Spline Interpolation ->Multilevel B-spline Interpolation (from Grid). The default parameters of this tool should be left as it is unless you know better. If some parts of the buildings remained on processed raster you will get some undesired hills on the resulted DEM so inspect the result carefully. Now you have a DEM, but you might wont to "polish" it further.

Polishing the DEM

By polishing I mean smoothing actually. Definitely there will be a lot of relatively rough terrain left after subtraction of elevated areas. To remove needless junk and to further smooth borders of interpolated areas I would suggest to use Modules -> Grid-Filter -> Gaussian Filter. The higher the value of Standard Deviation and Search Radius you will set - the smoother final DEM will be.

Here is the result:

Final DEM (the difference between red and yellow pixels on this sample is less then 1 metre)