Unit 18 - Lidar, DTM interpolation

LiDAR data used in this unit can be found in geodata/lidar directory. Data is stored in LAZ format (compressed LAS format).

Note

There is open Lidar data for Thuringia region available, see related blog by rapidlasso GmbH. Go to Download Höhendaten, choose region in map viewer and download selected tiles.

Data import

For importing LAS/LAZ data there are available two GRASS modules:

  • r.in.lidar for importing point cloud data as raster map
  • v.in.lidar for importing point cloud data as vector point map

Note

GRASS requires libLAS library for reading LAS data. For compressed LAS (LAZ) there is another dependency, a laszip library.

Basic metadata can be obtained by command line utility lasinfo which is part of libLAS library.

lasinfo las_681_5644_1_th_2014-2019.laz

...
Min X Y Z:                   681000.000 5644000.000 -43.382
Max X Y Z:                   681999.990 5644999.990 968.621
...

Tip

There are also new modules available r.in.pdal and v.in.pdal which are based on PDAL library.

Vector import

Let’s create a new GRASS location defined by EPSG:25832 and import selected tile into vector point map first. By Flag -t creation of attribute table is skipped. The import process will be significantly faster. Another flag which will speed up the import process is -b. It skips building topology (which is not needed for point features anyway). Information about data SRS can be missing, projection check can be skipped by -o flag.

v.in.lidar -otb input=/home/user/geodata/lidar/las_681_5644_1_th_2014-2019.laz output=las_681_5644

Note that computational region is ignored when importing data using v.in.lidar.

Overall point density can be checked by v.outlier. Note that v.outlier is working in the current computation region(!). Computational region must be set based on imported data. Since building topology was skipped by -b flag, the map extent must be computed by scanning all points in input vector map. That is reason why g.region can run for a while.

g.region vector=las_681_5644
v.outlier -e input=las_681_5644
Estimated point density: 8.937
Estimated mean distance between points: 0.3345

Note

Point density is calculated for map (square) unit.

Basic metadata can be printed by v.info (since no topology is built, the module must scan all features in the map):

v.info map=las_681_5644
...
|   Number of points:       8936470         Number of centroids:  0          |
...
|   Map is 3D:              Yes                                              |
...
|   Projection: UTM (zone 32)                                                |
|                                                                            |
|               N:        5644999.99    S:           5644000                 |
|               E:         681999.99    W:            681000                 |
|               B:           -43.382    T:           968.621                 |
...

Raster import

Input point cloud can be imported by r.in.lidar into raster map. Flag -e extends current computational region to cover all imported points. Otherwise user needs to set up computational region by g.region. Spatial resolution for output raster map is defined by resolution option. By default, for cells with more points involved, the value is computed by mean value, see method option. Cells covered by no points will get NULL values assigned.

r.in.lidar -oe input=/home/user/geodata/lidar/las_681_5644_1_th_2014-2019.laz output=las_681_5644 resolution=1

Basic metadata about created raster map can be obtained by r.info.

r.info map=las_681_5644
...
|            N:    5645000    S:    5644000   Res:     1                     |
|            E:     682000    W:     681000   Res:     1                     |
|   Range of data:    min = 65.51301  max = 346.671                          |
...
../_images/import-rast-vect.png

Fig. 104 Imported data as vector points. Raster map with 1m resolution in the background.

Filling gaps

There are several GRASS modules for filling gaps in raster maps like r.fillnulls or r.fill.stats. The first module is based on spline interpolation, the second fills gaps with interpolated values using IDW. Let’s use the second module which fill nulls rapidly compared to r.fillnulls. By -k flag it is ensured that original values will be kept untouched. Only cells with no-data value will be modified.

r.fill.stats -k input=las_681_5644 output=las_681_5644_no_gaps
../_images/rast-gaps-fill.png

Fig. 105 NULL values (on left part) filled by r.fill.stats (right part).

Note

Note that only cells in given distance (8 pixels by default, see cells option) are processed, see Fig. 106.

../_images/rast-gaps.png

Fig. 106 Cells out of distance not filled.

High resolution DSM

Digital Surface Model (DSM) can be interpolated by v.surf.rst using regularized spline with tension approximation. Output resolution set to 0.5 meter. The computation can be really slow. The computation time becomes more reasonable by running in parallel, see nprocs option (GRASS 7.4+ only).

g.region vector=las_681_5644 res=0.5 -pa
v.surf.rst input=las_681_5644 elevation=dsm_681_5644 npmin=80 tension=20 smooth=1 nprocs=5

Tip

Try also to set higher npmin to reduce artifacts.

../_images/dsm-3d.png

Fig. 107 DSM in 3D view. Orthophoto downloaded from Geoportal-Th.de (tile 32680_5644).

Tip

GRASS imports/links RGB image as separate bands. Color composition can be displayed using d.rgb. By r.composite it is possible to create color composite as a new raster map.

Note

Similarly can built Digital Terrain Model (DTM) from filtered terrain only points by edge detection.

The filter aims to recognize and extract attached and detached object (such as buildings, bridges, power lines, trees, etc.) in order to create a Digital Terrain Model. (source: v.lidar.edgedetection manual page) Example of simple workflow based on v.lidar.edgedetection, v.lidar.growing and v.lidar.correction below. Note that edge detection is usually a time and RAM consuming task, and the result is not perfect.

v.lidar.edgedetection input=las_681_5644 output=edge_681_5644 ew_step=8 ns_step=8 lambda_g=0.5
v.in.lidar -otb input=las_681_5644_1_th_2014-2019.laz output=las_681_5644_first return_filter=first
v.lidar.growing input=edge_681_5644 output=grow_681_5644 first=las_681_5644_first
v.lidar.correction input=grow_681_5644 output=corr_681_5644 terrain=terr_681_5644
../_images/terrain-only-points.png

Fig. 108 Filtered terrain points only.