Unit 19 - Lidar scripting

In this unit we will recover our knowledge of PyGRASS scripting (see Unit 14 - PyGRASS scripting). Let’s create a fancy script to produce seamless DEM (Digital Elevation Model) from bunch of LAS/LAZ files. Why so fancy? We will do our work in parallel!

Let’s start with defining user interface. There will be one input:

  • directory (line 6) - directory with input LAS/LAZ files

and one output:

  • elevation (line 9) - name for output elevation raster map mosaics

The resolution of output DEM will be defined by resolution parameter (line 13). And finally user will be able to control number of processes running in parallel by nproc (line 18).

Our script consists of three main functions:

1. import_files() to import input LAS/LAZ files (line 33). Import process can be done in parallel by ParallelModuleQueue, see lines 37, 42, 56, 57, 59.

2. create_dem_tiles() to compute DEM per tile (line 66) using v.surf.rst. We need to compute DEMs with reasonable overlap in order to create seamless mosaics, see 72-75. Tiles can be processed in parallel too, see line 81.

3. patch_tiles() to patch DEM tiles together by r.series, see 88. In overlaps cell values will be computed as average, it is a main reason why we do not use r.patch.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
#!/usr/bin/env python

#%module
#% description: Creates DEM from input LAS tiles.
#%end
#%option G_OPT_M_DIR
#% required: yes
#%end
#%option G_OPT_R_ELEV
#% description: Name for output elevation raster map mosaics
#%end
#%option
#% key: resolution
#% description: Output resolution
#% type: double
#%end
#%option
#% key: nprocs
#% description: Number of processes per tile
#% answer: 1
#% type: integer
#%end

import os
import sys
import time
from copy import deepcopy

import grass.script as gs

from grass.pygrass.modules import Module, ParallelModuleQueue

def import_files(directory):
    start = time.time()

    # queue for parallel jobs
    queue = ParallelModuleQueue(int(options['nprocs']))

    import_module = Module('v.in.lidar',
                           flags='otb',
                           overwrite=True,
                           run_=False
    )

    maps = []
    for f in os.listdir(directory):
        if os.path.splitext(f)[1] != '.laz':
            continue
        fullname = os.path.join(directory, f)
        basename = os.path.basename(f)
        # '-' is not valid for vector map names
        mapname = os.path.splitext(basename)[0].replace('-', '_')
        
        maps.append(mapname)
        gs.message("Importing <{}>...".format(fullname))
        import_task = deepcopy(import_module)
        queue.put(import_task(input=fullname, output=mapname))
    
    queue.wait()

    if not maps:
        gs.fatal("No input files found")

    return maps

def create_dem_tiles(maps, res, nprocs, offset_multiplier=10):
    offset=res * offset_multiplier

    for mapname in maps:
        Module('g.region',
               vector=mapname,
               n='n+{}'.format(offset),
               s='s-{}'.format(offset),
               e='e+{}'.format(offset),
               w='w-{}'.format(offset)
        )
        
        Module('v.surf.rst',
               input=mapname,
               elevation=mapname,
               nprocs=nprocs,
               overwrite=True
        )

def patch_tiles(maps, output, resolution):
    gs.message("Patching tiles <{}>...".format(','.join(maps)))
    Module('g.region', raster=maps, res=resolution)
    Module('r.series', input=maps, output=output, method='average', overwrite=True)
    Module('r.colors', map=output, color='elevation')

def main():
    start = time.time()

    maps = import_files(options['input'])
    create_dem_tiles(maps,
                     float(options['resolution']),
                     int(options['nprocs'])
    )
    patch_tiles(maps,
                options['elevation'],
                options['resolution']
    )

    gs.message("Done in {:.0f} min".format((time.time() - start) / 60.))
    
    return 0

if __name__ == "__main__":
    options, flags = gs.parser()

    sys.exit(main())