Unit 19 - DSM script parallelization

This unit is focused on parallel computing. Sample script below produces seamless DSM (Digital Terrain Model, see Unit 18 - Lidar, DSM interpolation) from bunch of LAS/LAZ files. Computation will be split into tiles and performed in parallel.

DSM interpolation in parallel

User interface defines two parameters, directory (line 6) for input directory with input LAS/LAZ files, and elevation (line 9) name for output elevation raster map mosaics. The resolution of output DSM is defined by resolution parameter (line 13). And finally number of processes running in parallel will be controlled by nproc (line 18) parameter.

A script consists of three main functions:

1. import_files() to import input LAS/LAZ files (line 33). Import process can be performed in parallel by ParallelModuleQueue from PyGRASS library (see Unit 11 - PyGRASS scripting for PyGRASS introduction), lines 37, 42, 54-55, 57.

2. create_dsm_tiles() to compute DSM per tile (line 64) using v.surf.rst. DSM tiles need to be computed with a reasonable overlap in order to create seamless mosaics, see 70-73. Tiles can be processed in parallel too, see nproc option on line 79.

3. patch_tiles() to patch DSM tiles together by r.series, see 83. From overlapping cell values is computed an average value. This is main reason why r.patch is not used here.

  1#!/usr/bin/env python3
  2
  3# %module
  4# % description: Creates DSM from input LAS tiles.
  5# %end
  6# %option G_OPT_M_DIR
  7# % required: yes
  8# %end
  9# %option G_OPT_R_ELEV
 10# % description: Name for output elevation raster map mosaics
 11# %end
 12# %option
 13# % key: resolution
 14# % description: Output resolution
 15# % type: double
 16# %end
 17# %option
 18# % key: nprocs
 19# % description: Number of processes per tile
 20# % answer: 1
 21# % type: integer
 22# %end
 23
 24import sys
 25import time
 26from copy import deepcopy
 27from pathlib import Path
 28
 29import grass.script as gs
 30
 31from grass.pygrass.modules import Module, ParallelModuleQueue
 32
 33def import_files(directory):
 34    start = time.time()
 35
 36    # queue for parallel jobs
 37    queue = ParallelModuleQueue(int(options['nprocs']))
 38
 39    import_module = Module('v.in.lidar',
 40                           flags='otb',
 41                           overwrite=True,
 42                           run_=False
 43    )
 44
 45    maps = []
 46    for fullname in Path(directory).glob('*.las'):
 47        basename = fullname.name
 48        # '-' is not valid for vector map names
 49        # vector map names cannot start with number
 50        mapname = Path(basename).stem.replace('-', '_')
 51
 52        maps.append(mapname)
 53        gs.message("Importing <{}>...".format(fullname))
 54        import_task = deepcopy(import_module)
 55        queue.put(import_task(input=str(fullname), output=mapname))
 56    
 57    queue.wait()
 58
 59    if not maps:
 60        gs.fatal("No input files found")
 61
 62    return maps
 63
 64def create_dsm_tiles(maps, res, nprocs, offset_multiplier=10):
 65    offset=res * offset_multiplier
 66
 67    for mapname in maps:
 68        Module('g.region',
 69               vector=mapname,
 70               n='n+{}'.format(offset),
 71               s='s-{}'.format(offset),
 72               e='e+{}'.format(offset),
 73               w='w-{}'.format(offset)
 74        )
 75        
 76        Module('v.surf.rst',
 77               input=mapname,
 78               elevation=mapname,
 79               nprocs=nprocs,
 80               overwrite=True
 81        )
 82
 83def patch_tiles(maps, output, resolution):
 84    gs.message("Patching tiles <{}>...".format(','.join(maps)))
 85    Module('g.region', raster=maps, res=resolution)
 86    Module('r.series', input=maps, output=output, method='average', overwrite=True)
 87    Module('r.colors', map=output, color='elevation')
 88
 89def main():
 90    start = time.time()
 91
 92    maps = import_files(options['input'])
 93    create_dsm_tiles(maps,
 94                     float(options['resolution']),
 95                     int(options['nprocs'])
 96    )
 97    patch_tiles(maps,
 98                options['elevation'],
 99                options['resolution']
100    )
101
102    gs.message("Done in {:.0f} min".format((time.time() - start) / 60.))
103    
104    return 0
105
106if __name__ == "__main__":
107    options, flags = gs.parser()
108
109    sys.exit(main())

Sample script to download: create-dsm.py

Note

The script is taking a long time with all the tiles from lidar directory. Choose few tiles for testing.

../_images/dtm_patched.png

Fig. 104 DSM created from all available tiles.

Canopy Height Model

In this session we will compute the Canopy Height Model (CHM), as the difference between interpolated DSM and imported EU-DEM (DTM) from Unit 15 - Data reprojection.

The CHM is computed using r.mapcalc, executing the difference between DSM and DTM.

r.mapcalc expression="chm = dtm_laz - dem"
../_images/chm.png

Fig. 105 The CHM map.

Task

Note that DSM created from Lidar data and EU-DEM have completely different qualitative characteristics. For a reasonable result it would be necessary to substitute DEM from another qualitatively similar source like DSM.