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.

Download sample data for this unit: jena-sample-data-lidar.7z (541MB)

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

There is also a new generation of Lidar modules based on PDAL library: r.in.pdal and v.in.pdal.

Important

Unfortunately on MS Windows PDAL-based modules are missing due to compilation problems (#1821) and libLAS-based modules lack compression support.

Data can be converted from LAZ to LAS by pdal command-line tool:

pdal translate las_681_5644_1_th_2014-2019.laz las_681_5644_1_th_2014-2019.las

Vector import

Let’s create a new GRASS location jena-lidar defined by EPSG:25832.

../_images/create-location-epsg.png

Fig. 103 Create a new GRASS location by specifying EPSG code (Select CRS from a list by EPSG or description).

Let’s import selected tile into GRASS as vector point map. By -t flag creation of attribute table is skipped. By -b flag also building topology will be skipped. The import process will be significantly faster. Information about data SRS can be missing, in this case use -o flag to skip projection check.

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

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 points 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 would need to set up computational region by g.region before importing data. 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 (no-data) value 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.

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=6

Úkol

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.

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.