Creation of DEMs from LIDAR point clouds

This guide assumes you have a LIDAR dataset comprising a number of files containing LIDAR lines of the following format:

217755.001130 30563577.97 4513423.12 1293.23   178 30563577.96 4513423.12 1293.28   178
time(gps?)(last pulse) ZZeasting northing height(m) pulse_time? (first pulse) ZZeasting northing height(m) pulse_time? 

1. Merge files and strip off unnecessary UTM zone number by running trim_lidar.sh in the directory containing the lidar files:

Usage: trim_lidar.sh output_file

2. Generate a histogram for the lidar data to determine min/max reasonable data values by running lidar_histogram.py:

Usage: python lidar_histogram.py input_file

This generates a 5x log histogram (ie it's a histogram of the base-10 log of the lidar values binned into 10m bins, multiplied by 5 to scale). If there are a small number of low or high values in the histogram then decide on a scaling cutoff. eg:

20: 
40: #################
50: #############################
60: #############################
70: ##############################
80: ##############################
90: ###############################
100: ###############################
110: ###############################
120: ################################
130: ################################
140: ###############################
150: #############################
160: ############################
170: ############################
180: ###########################
190: ###########################
200: ##########################
210: #######################
220: ###################
230: ##############
590: #####
600: ######
610: ######
620: ####

In the given example, there are a very few values between 20-30m at the bottom end, and a big gap at the top end before a small number of values above 590m (nothing between 240m and 590m). So it would be sensible to pick 30m as a minimum cutoff and 250m as a maximum. Note that because the histogram is on a log scale the difference in numbers of points at both ends is much larger than it looks at first glance.

3. Check the LIDAR data projection in the header file. Open GRASS, select (or create) an a location in the appropriate projection and set the region appropriately for the dataset - easiest way is to look at the start of the trimmed lidar data file for the starting co-ordinates of the first line, set an approximate region from there, view the map in GRASS and then adjust the region appropriately.

4. Read in the trimmed lidar file using r.in.xyz:

   r.in.xyz input=<trimmed_lidar_file> output=<lidar_basemap> x=2 y=3 z=4 fs=" " zrange=<min_cutoff>,<max_cutoff>
   

5. Interpolate gaps in the map using r.fillnulls:

   r.fillnulls input=<lidar_basemap> output=<lidar_interp_map>
   

6. Reproject the dataset if necessary to match the projection of the flight data (should be the same, but conceivably won't be). See Creation of DEMs from SRTM 90m data for how to do this.

7. Generate an ASCII DEM as per Creation of DEMs from NextMap data:

   r.out.ascii input=<lidar_interp_map> output=lidar.dem null=0