Version 28 (modified by mark1, 10 years ago) (diff)


Creating a DEM from a LIDAR point cloud

This page contains instructions for creating a simple DEM from the ASCII point clouds acquired from the Leica ALS50 II in 2009 onwards. Instructions are given for using your ARSF data in the GRASS GIS system.

There is also a section giving instructions on how to make your DEM suitable for use in the azgcorr software.

This assumes that you have removed any noisy points that you do not wish to be used in the making of your DEM (see Processing/LidarNoisyPoints).


The first step is to create a location of the required area. See the GRASS reference manual for help with this. Alternatively you can download a template grass database from here: GRASS template database. The template database contains both a UK National Grid and WGS84 Latitude/Longitude region. Start GRASS using a region suitable for the LiDAR data (creating one if it doesn't exist).

When this is done the ASCII point data can be loaded in. To import the ASCII data it must be within the GRASS region limits, any data outside the current region will not be imported.

Scan through the point cloud and find the min/max Eastings and Northings: -s input=<LDRfilename> output=<outputmapname> x=2 y=3 z=4 fs=' ' 

where <LDRfilename> is the filename you want to read in, <outputmapname> is the name you wish to call this within GRASS, x, y, and z are equal to the column numbers which contain the Easting, Northing and elevation values, fs is the field separator. For more info on this command see the GRASS manual

Set the region such that it contains all the point cloud data, and the resolution you wish to use – in this case 2.0m:

g.region  n=max_northing s=min_northing w=min_easting e=max_easting res=2.0

replacing the keywords max_northing, min_northing, min_easting, max_easting with the values from the command above. You may wish to add a small buffer on to ensure points on the boundary are imported.

Now we can import the laser point cloud data into GRASS. This uses the same command as above but without the -s flag: input=<LDRfilename> output=<outputmapname> x=2 y=3 z=4 fs=' '

The above 3 steps need to be repeated for each of the point cloud files you wish to use to make the DEM with.

When all point cloud files are loaded into GRASS, the region needs to be changed such that it covers the area of all the point clouds:

g.region rast=mapname1,mapname2,....,mapnameN

where mapname1, mapname2, etc are the map names of the imported point cloud data from the above steps.

All the separate LIDAR maps can now be patched together to make a single large raster data set. To do this, use the following command:

r.patch in=mapname1,mapname2,... out=lidar_mosaic

where the map names are as before, the imported point cloud rasters, and lidar_mosaic is the output name for the single concatenated raster set. For more details on this command see

Point cloud mosaic - gridded 2m

Figure: Mosaic of point cloud data gridded at 2m resolution. Holes (no data) are shown as white

If small holes are present in your data, say due to small lakes or removed noise, the easiest way to fill these in is to patch the data with a lower resolution version of the LiDAR. This can be done by changing the region to a lower resolution and resampling:

g.region res=5
g.resamp.stats input=lidar_mosaic output=lowres_mosaic

which in this example will create a new raster called lowres_mosaic at 5m resolution. Note that this is not suitable for larger holes in the data and other methods should be used in this case such as patching external data or interpolating.

6. r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=- input=lidar_mosaicx100 output=lidar_mosaic_idwx100

r.mapcalculator amap=lidar_mosaic_idwx100 formula=0.01*A outfile=lidar_mosaic_idw help=-

where the output raster, lidar_mosaic_idw, has been interpolated using an inverse distance weighted formula. As well as filling in holes within the lidar swath, this will also interpolate over empty parts of the GRASS region, resulting in possibly unrealistic data values.

interpolated full region

Figure: Interpolated point cloud data over full GRASS region.

This can be improved if you wish by using a mask when interpolating the data. To create a mask, the easiest way is to use your input data at a low resolution. This creates a raster which covers the LIDAR swath but contains no holes in the data (assuming the resolution is selected low enough). To import the data at a lower resolution , repeat steps 1-3 for each ASCII point cloud setting the res variable to a suitable value, e.g. 50.0 and outputting to a new map. Then repeat steps 4 and 5 to create a raster covering the combined LIDAR swath.

Then to use the mask:

  1. r.mask input=<maskmapname>

and perform the interpolation step (with precision maintaining commands before and after as above):

  1. Re-set the resolution to 2 for the LiDAR before any calculations.

g.region res=2

9. r.mapcalculator amap=lidar_mosaic formula=100*A outfile=lidar_mosaicx100 help=- input=lidar_mosaicx100 output=lidar_mosaic_idwx100

r.mapcalculator amap=lidar_mosaic_idwx100 formula=0.01*A outfile=lidar_mosaic_idw help=-

this results in a map where the internal holes have been filled, but the area outside the swath coverage remains unchanged.

interpolated with mask

Figure: Interpolated point cloud data after implementing a mask.

To extend the coverage of your DEM, if required, it is suggested to patch on external DEM data. If you have access to a good quality DEM then use that, else the SRTM 3 arc second DEM is freely available and covers most of the globe between +-60 degrees latitude. To see how to make a DEM from SRTM data see the SRTM DEM page. Make sure to select the projection the same as your LIDAR DEM. Also, be aware that the SRTM elevations are with respect to a geoid model and will need to be converted to the LIDAR vertical datum.

Once you have an SRTM DEM of sufficient coverage you can patch the LIDAR and SRTM DEMs together, such that the LIDAR takes precedence. This means the lidar_dem should be the first of the input maps on the r.patch command. If the mask is still applied then remove it before patching:

  1. r.mask input=<maskmapname> -r
  1. r.patch in=lidar_dem,srtm_dem out=combined_dem

Combined SRTM and Lidar DEM

Figure: Interpolated point cloud data after implementing a mask with SRTM DEM data to fill in the rest of the region.

To output the DEM (as ascii) and set null values as 0:

  1. r.out.ascii input=<DEMmapname> output=<outputDEMfile> null=0

Making the DEM suitable for azgcorr

To make a suitable ASCII DEM for use in the azgcorr software, the header information of the ASCII DEM file must be in a certain format. The required format is to have a header of one line (the first line of the file) with the DEM data following. The format is (as given by the azgcorr help):

or c r xm ym xx yx gi


or = row order of the data [0 if South to North, 1 if North to South]

c = number of columns

r = number of rows

xm = minimum easting

ym = minimum northing

xx = maximum easting

yx = maximum northing

gi = grid size (spacing)

An example header, for a DEM of 2000 rows x 2000 columns at 5m resolution, might be: 1 2000 2000 400000 850000 410000 860000 5

Attachments (4)

Download all attachments as: .zip