Using LIDAR ground points to generate a DEM (Digital Elevation Model) and Contours using PDAL and SAGA

Global elevation data is available through NASA’s Shuttle Radar Topography Mission (SRTM) program. For those interested, this blog post compares the resolution of elevation slope maps derived from the SRTM project vs. LiDAR data, and provides a brief background behind the two methods. LiDAR, in short, provides for much higher resolution Digital Elevation Models (DEMs). Here I provide an extensive example of generating contour lines from LIDAR point clouds with classified ground points.

1. I have selected 4 LIDAR adjacent tiles of an area with Olympic National Park, Washington, downloaded from the Puget Sound LIDAR Consortium:
   • http://pugetsoundlidar.ess.washington.edu/lidardata/restricted/las/pslc2012/hoh/delivery2/q47123g8207.laz
   • http://pugetsoundlidar.ess.washington.edu/lidardata/restricted/las/pslc2012/hoh/delivery2/q47123g8208.laz
   • http://pugetsoundlidar.ess.washington.edu/lidardata/restricted/las/pslc2012/hoh/delivery2/q47123g8212.laz
   • http://pugetsoundlidar.ess.washington.edu/lidardata/restricted/las/pslc2012/hoh/delivery2/q47123g8213.laz
2. In order to download data from PSLC, you’ll need to register (free). Once you’ve downloaded your tiles of interest, the first step is filtering out only the points classified as ground. In forested settings, generating a DEM on the unfiltered .las tile would generate a Digital Surface Model, or DSM, in contrast to a Digital Terrain Model, or DTM. This can be done using PDAL’s filters.range module, because ground points in a classified point cloud are given the classification of 2. Note that not all point cloud .las files have classified points; methodologies for classifying unclassified points is beyond the scope of this article. Save a new text file with the following contents and name it filter_ground.pipeline.json:

   {
     "pipeline":[
       {
         "type":"readers.las",
         "compression": "laszip", 
         "filename":"q47123g8207.laz"
       },
       {
         "type":"filters.range",
         "limits": "Classification[2:2]"
       },
       {
         "type":"writers.text",
         "order":"X,Y,Z",
         "keep_unspecified": "false",
         "quote_header": "false",
         "delimiter": ",",
         "filename":"q47123g8207_Ground.txt"
       }
     ]
   }
   

3. Run this file by running pdal pipeline filter_ground.pipeline.json. If you have multiple tiles as I do in this example, you’ll need to run this filter for each file. You can simplify modify the input and output file arguments and save over the same file.
4. Import the outputted .txt point cloud files into SAGA GIS by following steps 4-8 from this post.
5. From the Geoprocessing menu, select Shapes -> Point Clouds -> Tools -> Merge Point Clouds to merge the point clouds you have imported into SAGA.
   
6. From the Geoprocessing menu, select Grid -> Gridding -> Interpolation from Points -> Triangulation. Fill out the dialog with the information and point cloud you have created from the merge.
   • Set Attribute to Z.
   • Set Fit to cells.

   

7. Click Okay. Running the gridding algorithm will likely take some time, potentially several hours. This is because the implementation does not leverage multiple processor cores. Once completed, you should see your generated grid similar in the screenshot below.
   
8. Now, from the Geoprocessing menu, select Shapes -> Grid -> Vectorization -> Contour Lines from Grid. Set Equidistance to the interval to which you would like contour lines rendered.
9. After your contour lines are created, you can optionally have them labeled via SAGA’s object properties dialog, and/or export the contour lines for use in outside programs.