- See also
- python.packages.forDSM
Aim of script
calculates a raster DSM from point cloud data, optimized for forested areas.
General description
The python script forDSM implements the DSM derivation approach of Hollaus et al, 2010, referred to as "land cover dependent DSM". The surface model is calculated as a combination of a raster containing the highest data point within each cell (DSMmax) and a grid obtained from moving least squares interpolation (DSMmls) with a tilted plane as base model (i.e. moving planes interpolation). The moving planes interpolation has the following advantages compared to the mere rasterization of the DSMmax model: (i) a more flexible point neighbourhood is considered by extending the selection area beyond the domain of the grid cell, and (ii) the sigma0 of least squared interpolation provide a roughness measure. The latter is calculated from the residuals of the plane adjustment (deviation between the heights of the data points and the heights of the best fitting plane at the x/y locations of the data points). High roughness values are hereby used as an indicator of vegetation areas and/or immediate height jumps in the terrain (e.g. at roof boundaries) as the points from both the ground and the canopy (or roof, power line...) are used in the plane interpolation. In these areas the DSMmax is used in the final DSM to to avoid underestimation of the canopy heights. Small roughness values, in turn, indicate smooth surfaces like bare ground or roofs. In this case the smoother DMSmls is preferred to avoid artificial height jumps in the final product. More details about the approach can be found in Hollaus et al, 2010.
Basic algorithm
The algorithm consists the following steps:
- Raster based homogenization of the data distribution using Module Cell , by selecting the highest points within a relatively fine raster (i.e. ca. half of final grid size of the target DSM). Although the processing is raster based, the outcome of this steps is a subset of the original data points which are stored in an OPALS Datamanager (ODM) file. This odm-file is used as input for the following two processing steps.
- Calculation of DSMmax raster model extracting the highest elevations within each raster cell from the ODM using Module Cell . If no data points are availabe in a cell, the respective cell is marked as a void (no data).
- Calculation of DSMmls grid model using Module Grid . For each grid point the height is interpolated by estimating a best fitting tilted plane based on the k-nearest data points in a maximum search radius around the grid point. In addition to the grid heights, a roughness measure (standard deviation σ0 of the plane adjustment) is estimated and stored in a separate raster.
- Calculation if the final land cover dependent DSM by merging DSMmax and DSMmls based on the σ0 roughness layer using Module Algebra . DSMmax is used in rough areas (i.e. σ0 greater than a specified threshold) and DSMmls otherwise. In addition, data voids in the DSMmax model are filled with DSMmls values.
Usage of script
The input (-input) can either be a single ODM file or a list of ODM files (e.g. an entire directory containing multiple ODMs). In case multiple ODMs are specified, the ODMs should not contain redundant data as proper data overlaps are for avoiding border discontinuities are considered by the script internally. This especially applies to large, tiled datasets (e.g. entire flight blocks, project areas, or even country-wide data) which should be provided without overlap. All temporary products (intermediate ODMs and rasters) are stored in a folder which can be specified with -temp. The output folder for storage of the final DSMs is specified by -output. By default, a single land cover dependent DSM is calculated from the entire area of all input data set. However, a user defined limit window (xmin ymin xmax ymax) can be specified via -limit and tile based processing is enabled, when either specifying a tile size (-tileSize) or by providing a shape file containing a set of (rectangular) polygons. In case of tile based DSM creation a comprehensive DMS mosaic can be generated by setting parameter -mosaic to "True".
Please note, that all specified paths are interpreted relative to the project directory (-projectDir) unless absolute paths are provided. Already existing temporary products are either calculated anew or re-calculation is skipped depending on the state of parameter -skipIfExists.
Parameter description
-threshold Threshold for model calculation
Type : Floating-point number
Remark : mandatory
Description: Threshold for model calculation. If the simga0 map is greater than the threshold, the combined model will use the "max" value at that point.
-mosaic Create DSM mosaic from all individual DSMs (default: Boolean(False))
Type : Boolean
Remark : optional, default: Boolean(False)
Description:
Possible values:
| 1 | 1 |
| 0 | 0 |
| true | true |
| false | false |
| yes | yes |
| no | no |
| Boolean(True) | Boolean(True) |
| Boolean(False) | Boolean(False) |
| True | True |
| False | False |
-limit Limit of ouput
Type : String
Remark : optional
Description: Limit output to a specific region ("lower,left,upper,right"). Alternatively, a path to a (polygon) shapefile to read shapes of interest from. (Values will only be calculated for points within the union of all polygons)
-tileSize preTM tiling string
Type : String
Remark : optional
Description: preTM tiling string - either one of: one or two floats (setting a tilesize), a shapefile (with a unique column name given in -a), rows/cols like "3r2c", "scan" to copy tiling from input files, None
-attribute preTM naming string (default: LL_6)
Type : String
Remark : optional, default: 'LL_6'
Description: preTM naming string - depening on the tiling string, this can be "shp.XX" where XX is the name of a unique field in the shapefile, "inp" to copy filenames from the tiling coming from the input files, any of "LL", "LR", "UR", "UL" in combination with a delimiter and a number of digits (e.g. "LR_6" for 6 digits and "_" as delimiter, "INDEX" or "NUMBERED" for adressing via count.
-enduranceMode Enables/Disables the endurance mode. (default: Boolean(False))
Type : Boolean
Remark : optional, default: Boolean(False)
Description: If endurance mode is enabled the program will not stop if it encounters an errorin a tile. Instead it finishes all tiles as best as possible.
Possible values:
| 1 | 1 |
| 0 | 0 |
| true | true |
| false | false |
| yes | yes |
| no | no |
| Boolean(True) | Boolean(True) |
| Boolean(False) | Boolean(False) |
| True | True |
| False | False |
-gridSize Set the grid size of the output raster file. (default: 1)
Type : Floating-point number
Remark : optional, default: 1
Description:
-searchRadius Set the search radius for MSL raster map. (default: 10)
Type : Floating-point number
Remark : optional, default: 10
Description:
-neighbours Set the number of nearest neighbours for grid point interpolation. (default: 8)
Type : Floating-point number
Remark : optional, default: 8
Description:
-overlap Set overlap for tiling (Default: searchRadius [-r])
Type : Floating-point number
Remark : estimable (rule: searchRadius)
Description: Set an overlap for tiling differing from default. Be sure to set this to match the filenames in case of already tiled input, and to ensure enough overlap, if multiple scripts are to be executed (e.g. forDelineate) (Estimation rule: searchRadius)
-nbThreads The number of threads created by preCutting. Use 0 for maximum number of threads. (default: 1)
Type : Integer
Remark : optional, default: 1
Description: The number of threads created by preCutting. Setting it to 0 creates the maximumamount of threads supported by this machine. Setting it to 0 or a high number may impact the stability of the computer.
Common Package/Script Options
Settings concerning the general options for (package) scripts.
-infile input text file, directory or strip file
Type : String
Remark : mandatory
Description: As input either the name of an ASCII file (extension .txt) containing the names of the data files (one file per line) or an expression using wildcards to specify the respective data sets are accepted.
-outfile output directory or settings file (.txt)
Type : String
Remark : optional
Description: Specifies the directory where the final results are stored. If not specified, the results are created in the current working directory.
-temp temporary directory (default: current directory)
Type : String
Remark : optional, default: '.'
Description: Specifies the directory where all intermediate results are stored. If not specified, a TEMP subdirectory is created in the current working directory. Additional subdirectories (one per involved module) are created.
-cfg configuration file
Type : String
Remark : optional
Description: The name of a configuration file containing all relevant calculation parameters is expected. If not specified, the default cfg files as provided by the OPALS distribution are used.
-projectDir project directory
Type : String
Remark : optional
Description: The default project directory is the current working directory, but can be easily changed by this parameter. All path parameters, if not specified as absolute paths, are interpreted relative to the project directory.
-skipIfExists skip processing if result already exists (default: Boolean(True))
Type : Boolean
Remark : optional, default: Boolean(True)
Description: Skip processing if result already exists. In order to re-run current script it is useful to repeat the processing only if the respective output does not already exist. This allows for incremental processing of large projects.
Possible values:
| 1 | 1 |
| 0 | 0 |
| true | true |
| false | false |
| yes | yes |
| no | no |
| Boolean(True) | Boolean(True) |
| Boolean(False) | Boolean(False) |
| True | True |
| False | False |
-gearman address and port of gearman job server
Type :
Remark : optional
Description: currently this option is not yet supported. placed here for future reference
Logging Options
Settings concerning the verbosity level of logging.
-screenLogLevel verbosity level of screen output (default: info)
Type : fromStr
Remark : optional, default: info
Description: The minimum level of log message importance to print on the screen
-fileLogLevel verbosity level of log file output (default: info)
Type : fromStr
Remark : optional, default: info
Description: The minimum level of log message importance to export to the xml-log
Examples
The data used in the following examples can be found in the $OPALS_ROOT/demo/ directory.
As a prerequisite for the subsequent examples, please import the data using the following commands:
Example 1
This example generates a simple, tiled DSM, and additionally a mosaic:
opals forDSM -infile strip11.odm -tileSize 100 -mosaic True -
threshold 0.2
Example 2
When using a shapefile as argument for the -tileSize parameter, a separate DSM is created for each polygon of the shapefile. -attribute sets the output filenames to be determined by the "polygon_id"-Attribute in the shapefile.
forDSM -i strip11.odm -tileSize stripsPolygons.shp -temp forDSM_temp -
threshold 0.2
Example 3
This example shows how to override settings made in the configuration file (.cfg). For the grid calculations, 8 neighbours and a search radius of 10 m are used. The output cell size will be 0.2 meters. Since tiling is not activated, the output will be a single raster file in GeoTiff format.
forDSM -i strip11.odm -searchRadius 10 -gridSize 0.2 -neighbours 8 -temp forDSM_temp -
threshold 0.2
Example 4
The Python script forDSM.py can take and process multiple input files at once. The following example shows two different ways to do so, using the strip*.laz files of the same flight block.
forDSM -i strip1?.odm -
threshold 0.2 -gridSize 0.1
forDSM -i strip11.odm strip21.odm strip31.odm -thr 0.2 -gri 0.1
In the first case, all of the files are chosen using the wildcard character (*), while in the second case only specific files are selected using a whitespace as separator. In the figure that follows, the derived 0.1m DSM grid is shown (using all strips).
Figure 1: DSM derived from the six strip*.laz files
Example 5
This example gives some guidance on how to use the limit window (-limit) in order to select a subset of the entire area. Here, the full-waveform file fullwave.fwf is used.
forDSM -i fullwave.odm -limit "24862 311179 24970 311240" -th 0.2 -gr 0.8
The resulting 0.8m DSM is illustrated in Figure 2:
Figure 2: DSM derived from a subset area of dataset fullwave.fwf
Example 6
A way to facilitate the interpretation of a DSM is to derive the normalized DSM (nDSM). A nDSM is the difference model "DSM minus DTM" and, therefore, represents heights above the terrain. In this example, it is shown how to create the nDSM, when the DTM is provided. The DSM is derived again from the six .laz files used previously in Example 4 and the DTM is already provided in the $OPALS_ROOT/demo/ directory.
forDSM -i
strip*1.odm -th 0.2 -gr 0.1
opalsAlgebra -inFile forDSM/forDSM_all.tif strips_dtm.tif -
formula r[0]-r[1] -outFile nDSM.tif
The resulting nDSM (Fig.3) has gridsize of 0.1 meters and it was calculated using the module opalsAlgebra.
Figure 3: nDSM derived from the six strip*.laz files
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
1.8.17