Nick M Cummins

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.
   

Recent versions of SAGA GIS no longer support importing files in the .las format. However, a .las file can be converted using PDAL’s writers.text module into a text file format, which in turn can be imported in SAGA. For this example, we will begin with a .las file projected in a state plane CRS that is the output of my post Reprojecting LIDAR .las files to Latitude/Longitude to State Plane Example.

1. Open a new text file named las2xyz.json in the same directory as your .las file and paste the following into it:

   {
     "pipeline": [
       {
         "type": "readers.las",
         "filename": "20131013_usgs_olympic_wa_10TDT310990.las"
       },
       {
         "type": "writers.text",
         "order": "X,Y,Z",
         "keep_unspecified": "false",
         "quote_header": "false",
         "delimiter": ",",
         "filename": "20131013_usgs_olympic_wa_10TDT310990.txt"
       }
     ]
   }

2. Now run the pipeline. This may take a while. Once it completes, you should now see the converted .txt file in the working directory.

   pdal pipeline las2xyz.json
   ls -hs -1 
   total 1.1G
   227M 20131013_usgs_olympic_wa_10TDT310990.las
   148M 20131013_usgs_olympic_wa_10TDT310990.laz
   745M 20131013_usgs_olympic_wa_10TDT310990.txt
   4.0K las2xyz.json
   4.0K lat_lng_WA_S_FIPS_ft_reprojection.pipeline.json
   

3. Let’s perform a quick sanity check on this new file:

   less 20131013_usgs_olympic_wa_10TDT310990.txt
   X,Y,Z
   803169.892,932562.842,2093.730
   803170.060,932566.158,2093.070
   803178.839,932573.593,2054.360
   803170.851,932564.554,2093.700
   803225.508,932562.372,2043.070
   803232.723,932560.928,2040.090
   803232.361,932562.222,2041.790
   803232.443,932564.117,2041.470
   803234.863,932568.358,2030.210
   803235.004,932568.206,2034.440
   803236.294,932567.895,2028.310
   
   

4. Open SAGA GIS. From the Geoprocessing menu, select File -> Tables -> Import -> Import Text Table with Numbers only
   
5. Select the .txt point cloud file that was created by PDAL, with the appropriate options as shown in the screenshot. Click Okay. Since the example text file used is nearly 1.0 GB in size, importing may take some time.
   
6. Once the text file has been imported as a table into SAGA, select Geoprocessing -> Shapes -> Point Clouds -> Conversion -> Point Cloud from Table.
   
7. Select the table that was created in the previous step, and then set the values for X/Y/Z to these columns from the table.
   
8. Afterwards, you will get the point cloud generated in your data sources.
   
9. From the Geoprocessing menu, select Visualization -> 3D Viewer -> Point Cloud Viewer. Select Z as the colored attribute.
   
   
10. Click Okay. The point cloud will be displayed in SAGA’s Point Cloud Viewer tool:
    

1. Download one or more lidar tiles:

   wget https://coast.noaa.gov/htdata/lidar1_z/geoid12b/data/5008/20131013_usgs_olympic_wa_10TDT310990.laz

2. Note the source crs (coordinate reference system, i.e. projection). This 2013 USGS Lidar: Olympic Peninsula (WA) Point Cloud dataset from NOAA uses EPSG 4269. This means that the point cloud coordinates are using latitude/longitude, which is not well supported by many downstream GIS tools. Running lasinfo https://coast.noaa.gov/htdata/lidar1_z/geoid12b/data/5008/20131013_usgs_olympic_wa_10TDT310990.laz confirms that the x and y coordinates are in lat/lng:

    
   lasinfo (170203) report for 20131013_usgs_olympic_wa_10TDT310990.laz
   reporting all LAS header entries:
     file signature:             'LASF'
     file source ID:             0
     global_encoding:            1
     project ID GUID data 1-4:   00000000-0000-0000-6C4F-6D7900636970
     version major.minor:        1.2
     system identifier:          'LAStools (c) by rapidlasso GmbH'
     generating software:        'lasduplicate (160119) commercia'
     file creation day/year:     231/2014
     header size:                227
     offset to point data:       331
     number var. length records: 1
     point data format:          1
     point data record length:   28
     number of point records:    25204541
     number of points by return: 16034492 6849061 2147980 173008 0
     scale factor x y z:         0.0000001 0.0000001 0.001
     offset x y z:               0 0 0
     min x y z:                  -123.9223149 47.8406412 248.190
     max x y z:                  -123.9087939 47.8497435 827.800
   

   Thus, re-projecting the .las file is necessary before any further processing can be done. Common CRS systems include UTM based projections, as well as state plane. For this example, I’ve chosen to reproject from lat/lng to NAD 1983 StatePlane Washington South FIPS 4602 Feet as it is what is used by data available through the Puget Sound Lidar Consortium.

Re-projection can be done with PDAL’s filters.reprojection. Create a new JSON pipeline file named lat_lng_WA_S_FIPS_ft_reprojection.pipeline.json with the following contents:

{
  "pipeline": [
    {
      "type": "readers.las",
      "compression": "laszip",
      "filename": "20131013_usgs_olympic_wa_10TDT310990.laz"
    },
    {
      "type": "filters.reprojection",
      "in_srs": "EPSG:4326",
      "out_srs": "EPSG:102749"
    },
    {
      "type": "writers.las",
      "compression": "laszip",
      "scale_x": "0.01",
      "scale_y": "0.01",
      "scale_z": "0.01",
      "offset_x": "auto",
      "offset_y": "auto",
      "offset_z": "auto",
      "filename": "20131013_usgs_olympic_wa_10TDT310990.las"
    }
  ]
}

Now save the file and run the pipeline. This will probably take a while as it involves de-compressing the input .laz file, iterating through and reprojecting every point in the point cloud, and writing the output file to disk. Now, examining the outputted, reprojected .las file with lasinfo 20131013_usgs_olympic_wa_10TDT310990.las will indicate the min/max in the new projection, as well as the z (elevation) values converted into feet:

lasinfo 20131013_usgs_olympic_wa_10TDT310990.las 
lasinfo (170203) report for 20131013_usgs_olympic_wa_10TDT310990.las
reporting all LAS header entries:
  file signature:             'LASF'
  file source ID:             0
  global_encoding:            0
  project ID GUID data 1-4:   00000000-0000-0000-0000-000000000000
  version major.minor:        1.2
  system identifier:          'PDAL'
  generating software:        'PDAL 1.5.0 (Releas)'
  file creation day/year:     211/2017
  header size:                227
  offset to point data:       690
  number var. length records: 3
  point data format:          3
  point data record length:   34
  number of point records:    25204541
  number of points by return: 16034492 6849061 2147980 173008 0
  scale factor x y z:         0.01 0.01 0.01
  offset x y z:               800114.623206489603035 932554.57963061647024 814.270025000000487
  min x y z:                  800114.62 932554.58 814.27
  max x y z:                  803498.24 935938.24 2715.87