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).

Please refer to this page on using scripts for creation of DEMS as it is a fully scripted approach. This page is designed more as a tutorial for making a DEM in GRASS.

---

Using GRASS

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.

Import the LiDAR

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:

r.in.xyz -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 ​http://grass.fbk.eu/gdp/html_grass63/r.in.xyz.html

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 r.in.xyz 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 r.in.xyz command as above but without the -s flag:

r.in.xyz 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 ​http://grass.fbk.eu/gdp/html_grass63/r.patch.html.

Dealing with small holes

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
r.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.

This can then be patched onto the original higher resolution LiDAR data such that the high resolution takes priority (as it is listed first in the 'in' option list).

r.patch in=lidar_mosaic,lowres_mosaic  out=high_low_patched

Dealing with larger holes - using 3rd party data

If there are large areas of no data then it is recommended that 3rd party DEMs are used to fill these. ASTER and SRTM DEMs have near global coverage and are free to download. This assumes using the ASTER data.

As this process is very similar to extending the region of the DEM please see how to do it [here].

Dealing with larger holes - interpolating existing data

If you do not wish to use 3rd party data or none exists for your area then it is possible it interpolate and extrapolate.

This can be done using the following command. Note that this can take a long time to complete depending on the amount of null data that is present in the raster and the resolution and size of the raster.

r.fillnulls input=raster_to_fill output=filled_raster tension=40. smooth=0.1 

which will fill in the null values from 'raster_to_fill' and store the new data in 'filled_raster'. The tension and smooth values are used to define the spline interpolation. Note that this method will spam the screen with a lot of "Warnings" which GRASS says can be ignored. See ​http://grass.fbk.eu/gdp/html_grass63/r.fillnulls.html for more information.

Alternative Method

If you wished to use an alternative method of interpolation, and not extrapolate outside the edges of the data, then you could follow a procedure similar to this.

First create a mask that covers the region of the data. 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 create lower resolution data, say at 20m, use:

g.region res=20
g.resamp.stats input=raster_name output=mask_name

To assign the mask and reset the region resolution (in this example to 2m) use:

g.region res=2
r.mask input=mask_name

To interpolate the data we will use 'r.surf.idw' which requires the data to be scaled first to retain precision.

This command will scale the data by 100:

r.mapcalculator amap=raster_to_scale formula=100*A outfile=rasterx100 help=-

Run the interpolation command:

r.surf.idw input=rasterx100 output=interpolatedx100

And rescale the data back

r.mapcalculator amap=interpolatedx100 formula=0.01*A outfile=interpolated help=-

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

Filtering the data

It may be wise to filter / smooth the data to remove localised spikes. This can be done using the following command:

r.neighbors input=raster_to_filter output=filtered_raster method=median size=3 --overwrite

which will apply a median filter of kernel size 3 to the raster_to_filter, resulting in filtered_raster. More information on the filters that can be applied can be found at ​http://grass.fbk.eu/gdp/html_grass63/r.neighbors.html.

Extending the coverage using ASTER data

If you wish to use the DEM for processing your hyperspectral data then it may be that you need to extend the coverage such that it covers the hyperspectral data. This can be done by patching on ASTER data. Information on this, and how to access ASTER data, can be found here, along with how to offset the ASTER data from geoidal heights to spheroidal heights.

Writing out the DEM

The DEM can be written out from GRASS in various formats (see ​http://grass.fbk.eu/gdp/html_grass63/raster.html functions starting 'r.out').

To write the data out as an ASCII grid (for example to use in azgcorr) you could use:

r.out.ascii input=raster_name output=output_filename null=0

which will write out the raster 'raster_name' to the file 'output_filename' and set any null values to 0. To use this DEM the header will need to be converted into a style azgcorr understands, such as here.

To write out the data as an ENVI BIL file (for example to use in APL) you could use:

r.out.gdal format=ENVI type=Float32 input=raster_name output=output_filename nodata=9999

which will write out raster 'raster_name' to a 32-bit floating point ENVI file named 'output_filename', using 9999 to represent null data.

IMPORTANT NOTE. If you wish to use the DEM in APL it MUST be in WGS84 Latitude / Longitude. If you wish to do this and you are not working in a WGS84 Latitude/Longitude region then you will need to reproject the data into that region before outputting.

---

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

where:

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