Module DSM
See also
opals::IDSM

Aim of module

Calculates a land cover dependent DSM raster model based on the approach of Hollaus et. al. (2010).

General description

In contrast to the algorithms for determining digital terrain models (DTMs) from Airborne Laser Scanning (ALS) data, those for the generation of Digital Surface Models (DSMs) are rather straightforward. A common way to derive DSMs is the determination of the highest point within a defined raster cell (DSMmax) or the determination of the DSM heights based on moving least squares interpolation e.g. moving planes (DSMmls). For low ALS point densities void pixels can occur in the DSMmax and for inclined smooth surfaces the DSMmax shows artificial roughness mainly due to the irregular point distribution. These disadvantages of the DSMmax can be reduced by using an interpolated DSMmls grid. On the other hand, the DSMmls introduces smoothing effects along surface discontinuities e.g. building borders, tree tops, forest gaps, etc. Module DSM therefore implements a combined approach for DSM generation making use of the strengths of, both, the DSMmax and the DSMmls (Hollaus et al., 2010).

Starting with the 3D point cloud stored in a single or multiple OPALS Datamanger (ODM) files as basic input (inFile) the final DSM (outFile) is derived as a regular grid (gridSize) and stored in a GDAL (GeoData Abstraction Library) supported format (oFormat). In general, the entire point cloud is used for processing but, however, the input data set can be restricted by specifying an appropriate filter string (filter). As stated above, the main idea is the calculation of (i) a raster containing the highest elevation within a raster cell (DSMmax) and (ii) a (smooth) surface grid derived via moving planes interpolation. For the latter the definition of a proper neighbourhood is required. Module DSM hereby uses a generation based neigbourhood defined by the desired number of data points (neighbours) used for plane fitting constrained by a maximum search distance (searchRadius). The point selection strategy is illustrated in Figure 1.

Fig.1: Generation based point selection

The generation based point selection works as follows. First, a sub-sampling grid is constructed (thin blue lines) such that each grid cell of the final DSM (thick blue lines and red circles, respectively) is composed of four sub-sampling grid cells. Within each cell of the sub-sampling grid the highest point is determined (black point surrounded by blue circle) and stored as the representative point of this cell. As a by-product also the lowest points are detected and and stored (black points surrounded by green circles). Depending on the point density and point distribution, some of the sub-sampling cells may not contain points and, thus, are marked as void. Please note, that although the (highest) points are now organized in a raster, allowing fast nearest neighbour point queries, still the original points with their arbitrary 3D positions and attributes are maintained. To interpolate the height of a specific grid point (green circle), the k-nearest neigbouring points are queried from the sub-sampling raster (highest points per cell) in a generation based manner. The first generation contains all cells whose cell centers are located in a circle with a radius equal to the DSM grid size. If the number of valid data points is less than the requested neighbour count (neighbours) the search radius is extended by the size of the sub-sampling grid (i.e. half of the DSM grid size). All (cell) points of the resulting circular ring (i.e. 2nd generation) are added. The procedure is repeated unit either the requested neighbour count or the maximum search radius is reached. The DSMmls grid height is subsequently interpolated by fitting a tilted plane (cf. Module Grid for further details). The DSMmax height, in turn, is strictly derived from the four sub-sampling grid cells sorrounding the grid point by selecting the highest cell. The DSMmax grid cell remains void if all four adjactent sub-sampling grid cells are void. The DSM height for the current position is finally obtained by choosing either the height of DSMmax or DSMmls, respecively, depending on the local surface roughness (maxSigma). The latter is obtained as a by-product of the moving planes interpolation. The DSMmax is chosen if a valid maximum height exists at the current location and if the local sourface roughness is greater than the specified threshold, otherwise the DSMmls is used. For more details please refer to Hollaus et al., 2010.

The specification of user defined limit window (limit), and a user defined void marker identifying invalid cells in the DSM output grid (noData) is optional. By default, the intermediate products (min, max, mls, sigma models) are stored as separate files on the disk. If, however, parameter multiBand is specified, the output file is constructed as a multi-band raster file containing the DSM (i.e. first band) and all intermediate rasters.

Parameter description

-inFileinput data manager file name(s)
Type: opals::Vector<opals::Path>
Remarks: mandatory
The path(s) to the opals datamanager(s) whose point data are used to derive the surface grid model. In case multiple ODMs are specified, the data are queried from all (potentially overlapping) ODMs.
-outFileoutput DSM grid file name
Type: opals::Vector<opals::Path>
Remarks: estimable
Path of output DSM grid file in GDAL supported format.
Estimation rule: The current directory and the name (body) of the input file are used as file name basis. Additionally, the postfix '_dsm' and the default extension of the output format are appended.
-oFormatDSM grid file format [GTiff,AAIGrid...]
Type: opals::String
Remarks: estimable
Use GDAL driver names like GTiff, AAIGrid, USGSDEM, SCOP... .
Estimation rule: The output format is estimated based on the extension of the output file (*.tif->GTiff, *.dem->USGSDEM, *.dtm->SCOP...).
-tileSizetile (block) size
Type: uint32
Remarks: default=64
The internal processing is organized in rectangular tiles. The size of these computing units(number of grid lines per tile) can be specified allowing control of the memory consumption of the module. Recommended values are in the range of 64 (default) if the average point distance approximately matches the grid size. In case coarser of finer grids are to be derived, the tile size may need to be adapted. Finer grids allow for larger tile sizes (and vice versa).
-gridSizeDSM grid width x/y
Type: double
Remarks: default=1
Size of the rectangular grid cell in x- and y-direction. If only one value is passed, quadratic cells are used.
-neighboursminimum nr of neighbours for grid interpolation
Type: int32
Remarks: default=8
The neighbour search is performed raster based in 'generations'. First, a regular raster featuring half of the final cell size is derived (storing the maximum elevation in each cell). For the grid interpolation (cf. general description for deatils) all (raster) points within a search radius of the final grid size are taken into account first. If there are less valid cells available than the required number of neighours, the search radius is extendended by half of the grid size. This procedure is repeated until either the required neighbour count is exceeded or the maximum search radius is reached.
-searchRadiusmaximum search radius for point selection
Type: double
Remarks: estimable
Only points within the given search radius are considered for the interpolation of a single grid post. If the search area contains too few points for successful interpolation, the respective grid point is omitted (void).Estimation rule: searchRadius = 3 * grid size.
-maxSigmaterrain roughness threshold
Type: float
Remarks: default=0.25
This parameter describes the maximum allowed roughness derived as the std.dev. of the plane fit carried out for each grid point. If the std.dev. is below the threshold or if no max-value is available for this grid point, the interpolated point is used in the final DSM. Otherwise the the maximum height within the grid cell is preferred.
-filtermodify the input using a (tree of) filter(s)
Type: opals::String
Remarks: optional
A filter string in EBNF syntax as described in section 'Filter syntax' can be passed to restrict the set of input points (e.g. to consider last echoes only).
-limit2D clipping window
Type: opals::GridLimit
Remarks: optional
If no user defined limits are specified or -limit is even skipped, the OPALS Data Manager (ODM) extents are used. Rounding is enabled by default in case the limits are derived from the ODM and can be enforced for user defined limits (keyword: round).
GridLimit: Defines xy-limits for output grid/raster datasets. The limit coordinates specified by left/lower/right/upper either refer to pixel centers (default) or corners, and may optionally get rounded to multiples of the grid size.
Syntax: ['center'|'corner']['round']['(' left lower right upper ')']
-noDatanodata value indicating void pixels
Type: opals::NoDataType
Remarks: default=max
Value representing an undefined value in the output DSM grid model.
NoDataType: Defines or disables (use 'none') the nodata value for the output raster. The nodata values can be a specific value or a type independent label, such as, 'min', 'max' and 'nan'. Those labels will be translate to values based on the raster band type. See the OPALS Docu (opals::NoDataType) for more Details.
Syntax: none | min | max | nan | <value>
-multiBandCreate multi-band output raster
Type: bool
Remarks: default=0
When activated the resulting DSM grid and all intermediate products (min, max, mls, sigma0 models) will be stored as multiple bands of the same raster file instead of multiple files. Features will be in the above order.
-extrapolationCheckavoid extrapolation flag
Type: bool
Remarks: default=1
Unfavorable point distributions can lead to extrapolations, where the interpolated value clearly exceeds the local data values. Using the local data range extended by a buffer factory, such extrapolation are detected and rejected if this option is enabled.

Examples

The data used in the following examples can be found in the $OPALS_ROOT/demo/ directory. As a prerequisite import the data (strip??.laz) into an OPALS data manager using the following command:

opalsImport -inFile strip??.laz -outFile dsmDemo.odm -tileSize 75 -coord_ref_sys EPSG:25833

Using the default values, a DSM raster in 1m-resolution can be derived with:

opalsDSM -inFile dsmDemo.odm -outFile dsm.tif

Please note, that in addition to the DSM raster (dsm.tif) additional intermediate raster products (max/min/movingPlanes/sigma0/pointCount-models) are created.

More examples coming soon ....

References

Hollaus, M., Mandlburger, G., Pfeifer, N., Mücke, W., 2010. Land cover dependent derivation of digital surface models from airborne laser scanning data. PCV 2010, Paris; in: "IAPRS Volume XXXVIII Part 3A", ISSN: 1682-1750; 221

Author
gm
Date
14.11.2016
@ coord_ref_sys
default coordinate reference system (EPSG Code, WKT string or PRJ-File)
opalsImport is the executable file of Module Import
Definition: ModuleExecutables.hpp:113
@ strip
strip ID of an image (opalsStripAdjust)
opalsDSM is the executable file of Module DSM
Definition: ModuleExecutables.hpp:43