# Unit 19 - DTM script parallelization¶

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

## DTM 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 DTM 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, 57-58, 60.

2. create_dtm_tiles() to compute DTM per tile (line 67) using v.surf.rst. DTM tiles need to be computed with a reasonable overlap in order to create seamless mosaics, see 73-76. Tiles can be processed in parallel too, see nproc option on line 82.

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

  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 112 #!/usr/bin/env python3 #%module #% description: Creates DTM 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 # vector map names cannot start with number 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_dtm_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_dtm_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()) 

Note

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

Fig. 109 DTM created from all available tiles.

## DTM comparision¶

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"


Fig. 110 The CHM map.