Unit 17 - Lidar

There is open Lidar data for Norway available at hoydedata.no. Go to Download, choose region in map viewer and request data for download.

LiDAR data used in this unit can be found in sample dataset (download 7z-archive), in lidar folder (681-682/5644-5646 tiles). Data is stored in LAZ (compressed LAS format).

Import data

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

Vector import

Let’s try to import selected tile into vector point map first. Flag -t skips creation of attribute table. 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, we will skip projection check by -o flag.

v.in.lidar -otb input=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.

We can also check the point overall point density by v.outlier. Note that v.outlier is working in the current computation region(!). So we also set up the region based on imported data (note that g.region can run for a while. Since we skipped building topology, the map extent must computed by scanning all points in input vector map).

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

Now let’s try import input points 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, see method option. Cells covered by no points will get NULL values assigned.

r.in.lidar -oe input=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                          |
...

#.. figure:: ../images/units/18/import-rast-vect.png

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. We will use the second module which fill nulls rapidly compared to r.fillnulls. By -k flag we ensure that original values will be kept. Only cells with no-data value will be modified.

r.fill.stats -k input=las_681_5644 output=las_681_5644_no_gaps

#.. figure:: ../images/units/18/rast-gaps-fill.png

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

Note that only cells in given distance (8 pixels by default, see cells option) are processed, see lidar-gaps.

#.. figure:: ../images/units/18/rast-gaps.png

Cells out of distance not filled.

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

#.. figure:: ../images/units/18/terrain-only-points.png

Filtered terrain only points.

High resolution DSM

Digital Surface Model (DSM) will interpolated by v.surf.rst using regularized spline with tension approximation. Output resolution will be set to 0.5 meter. The computation can be really slow. You can turn computation time to be more reasonable by running it 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.

#.. figure:: ../images/units/18/dsm-3d.png
class:middle

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

Note

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.

Similarly we can build Digital Terrain Model (DTM) from filtered terrain only points, see Edge detection for details.