Geomorphological application: Southwest Flank of Mount Rainier

Study area

Mount Rainier is frequently affected by various types of natural hazards. Those hazards range from lava flows, landslides, rockfalls, earthquakes, pyroclastic flows to lahars. The latter pose the greatest hazard from Mount Rainier. Lahars, also known as volcanic mudflows or debris flows, are large flows of eruption or hydrothermal origin that benefit from abundant ice, loose volcanic rock and surface water at Mount Rainier. Lahars at Mount Rainier vary in their fluid content but can travel as fast as 70-80 km per hour and may flow up to 100 km downstream to Puget Sound. At least 60 lahars have flowed off Mount Rainier into its draining river valleys in the past 10,000 years, with the largest events occurring about every 500 years (Pinsker, 2004). In response to the lahar-hazard, warning systems have been implemented. Before even thinking about such installations, the areas most susceptible to specific hazardous events have to be identified. Physically based modelling approaches are capable of predicting lahar runout. Such models have in common that they need proper and high quality topographic data. High resolution LiDAR produces raw data perfectly suitable for such analyses. In the following chapters some guidance is given for a typical workflow in geomorphological research to produce some ready-to-use topographic raster data with OPALS.

Data

LiDAR techniques have been extensively used for topographic surveys since the mid-1990s. Data acquisition was first performed by means of ALS (airborne laser scanning), whereas TLS (terrestrial laser scanning) emerged as a valuable field in earth sciences since approx. 2000 (Jaboyedoff et al. 2012). Existing geomorphological applications encompass landslide and rockfall dynamics (Teza et al., 2008; Abellán et al., 2009, 2010), coastal cliff erosion (e.g. Rosser et al., 2005), braided rivers evolution (e.g. Milan et al., 2007), river bank erosion (e.g. O’Neal and Pizzuto, 2011) or debris flows impacts (Schürch et al., 2011). Höfle & Rutzinger (2011) published an extensive article about the applications of airborne LiDAR in geomorphological applications.

Data acquisition by LiDAR techniques comes with the benefit of generating full 3D land surfaces. However, when it comes to analyses of multi-temporal datasets, the most common method of point cloud comparison in earth sciences is still the Digital Elevation Model (DEM) of difference (DoD) technique (Lague et al., 2013). That means that two registered point clouds are gridded to generate DEMs which are then differentiated on a pixel-by-pixel basis. Thus, the vertical distance between two corresponding raster cells is taken into account as a measure of surface change. Although a very fast technique that comes with a simple and well known data structure, it has the huge drawback of being only a 2.5D representation of the land surface, thus decreasing information density proportionally to surface steepness (i.e. vertical surfaces such as cliffs or rock faces) cannot be accounted for in a grid representation. Contrary, full 3D analyses are limited to direct point cloud processing techniques (such as C2C, C2M or M3C2) that do not depend on gridding or meshing the original data (Lague et al., 2013). Yet, applications that directly make use of the 3D point cloud are not as abundant, which can be explained by the lack of suitable software tools providing management and analysis functionality as well as direct integration with standard GIS and remote sensing software (Höfle & Rutzinger, 2011). Nonetheless, rasterized point cloud derivatives are appropriate by all means with respect to the research question raised, hence the dominance of raster-based applications.

In this use case we try to exploit the potential of OPALS as a tool for explorative data analysis and the preparation of raster data for subsequent model applications. Pfeifer et al. (2013) give some in-depth insight into the OPALS framework. For any question that arises while working with OPALS, an extensive use of the OPALS module descriptions that are available online is highly recommended. The data used here is obtained from OpenTopography, a high-resolution topography data repository that offers point cloud data under a free license (public domain). The raw data of the southwest flank of Mount Rainier can be directly downloaded from here. Fig. 1 shows the area covered by the point cloud. As for the point cloud download, the LAS file format was chosen (download in a generic ASCII file format available too). The survey that took place from August 28, 2012 to September 1, 2012, covers an area of 123 km². Processing the entire point cloud would be quite time demanding, thus, we only use a subset of the point cloud for fast and easy processing. You can download the subset used for this use case (UseCase_geomorph_MtRainer.7z) from the OPALS download page.

Figure 1: Area covered by point cloud at the southwest flank of Mount Rainier.

Derivation and analysis of topographic raster models

In a first step, we derive different topographic models from our 3D point cloud. The most frequently used model for geomorphological applications is the Digital Terrain Model (DTM) representing the elevation of the ground surface. In contrast to the DTM, a Digital Surface Model (DSM) represents the elevation of the topmost surface (e.g. treetops, roofs). A DSM can be derived by 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). In order to obtain a simple DTM, we will use the lowest point within a raster. Depending on the research question, more sophisticated vegetation removal tools might be required. Lastly, a simple method is applied to fill data gaps in the DTM. This is necessary for every physical based or numerical model applied later on that makes use of the DTM as the topographic base, especially when any kind of flow propagation (overland flow, river discharge, debris flow propagation, etc.) is involved.

Derivation of topographic raster models

The OPALS modules involved here are as follows:

  • opalsImport for generating a new ODM file from the raw point cloud file
  • opalsCell and opalsGrid to derive a DTM, DSMmax and DSMmls from the survey area
  • opalsShade for visualizing the previously generated topographic models as a hillshade

For generating a new ODM:

opalsImport -inFile mtrainier_points.laz -outFile mtrainier_points.odm

Generating a digital elevation model works with both, opalsCell and opalsGrid. The aim of Module Cell is to derive raster models by accumulating specific features (min, max, mean, etc.) of a selected input data attribute (z, amplitude, echo width, etc.), while opalsGrid derives digital elevation models in regular grid structure using interpolation techniques. The aim of Module Grid is to derive digital surface models (DSM/DTM) in regular grid structure using simple interpolation techniques like nearest neighbour or moving planes. This means that opalsCell derives only a single raster pixel value based on all valid values for a selected attribute within a raster cell (in our case the z-coordinate that represents the height), whereas opalsGrid interpolates over all values within a user-defined local neighbourhood. Consequently, the height estimation based on opalsCell may be more accurate in rough areas (e.g. canopy), but less robust with respect to errors resulting from outliers.

For derivation of GIS raster data we can now use:

opalsCell -inFile mtrainier_points.odm -outFile dsm_3m.tif -feature min -feature max -cellSize 3.0
opalsGrid -inFile mtrainier_points.odm -outFile dsm_grid_3m.tif -gridSize 3.0 -interpol movingPlane

Mind that you can specify different features (max, min) in one line. In case no attribute is specified, OPALS uses the default attribute z. The output files are named accordingly (dsm_3m_z_max.tif and dsm_3m_z_min.tif). For the DSM generation with opalsGrid we used a moving (tilted) plane interpolation, yet there are more interpolation methods available in OPALS.

For a better visualization of the digital elevation models, we use opalsShade (cf. Fig. 2). Hill shading is a method to visualize a surface grid model by illumination using an artificial light source shining from a certain direction. However, all subsequent analyses are based on the DSM/DTM only.

Shading a DSM/DTM is straight forward, yet there are some more options available in OPALS if required (with respect to the shading method or the artificial sun position).:

opalsShade -inFile dsm_3m_z_max.tif -outFile dsm_3m_max_shade.tif
opalsShade -inFile dsm_3m_z_min.tif -outFile dsm_3m_min_shade.tif
opalsShade -inFile dsm_grid_3m.tif -outFile dsm_grid_3m_shade.tif
Figure 2: The three hillshades in comparison from opalsCell min (left), opalsCell max (centre), opalsGrid (right).

Analysis of raster data

When working with physical based slope stability models, it is often advisable to separate forested from non-forested areas to account for the stabilizing effects of trees (intercepted precipitation, root cohesion, etc.). Thus, it would be interesting to extract all trees from the study area. Our data only contains x,y and z coordinates and no other attributes that give any information on a potential vegetation cover (e.g. echo with from full-waveform LiDAR data). Therefore we have to make an assumption, namely that we consider forested areas only as such with relative heights greater than 3 m with respect to the lowest point in the corresponding raster cell.

This requires two steps:

  1. Derive a difference raster between minimum and maximum laser point elevation in each cell (already calculated above with opalsCell: DTM and DSMmax) by applying opalsAlgebra
    opalsAlgebra -inFile dsm_3m_z_max.tif -inFile dsm_3m_z_min.tif -outFile ndsm_max_min.tif -formula "r[0]-r[1]" -gridSize 3.0</li
  2. Then extract all areas with relative height greater than 3 m. All other areas will be set to no data (invalid)
    opalsAlgebra -inFile ndsm_max_min.tif -outFile ndsm_max_min_gt3m.tif -formula "r[0]>3.0 ? r[0] : invalid" -gridSize 3.0

In step 1 we produce a normalized DSM (nDSM) by subtracting the DTM from the DSM. The aim of Module Algebra is to derive a new grid dataset by combining multiple input grid datasets. The new grid values are calculated by applying an algebraic formula using the grid values of the respective input grids. The formula is passed as a string in generic filter syntax. Apart from these basic operations, conditional statements are of crucial importance for the construction of complex formulas. For example in step 2 we want to generate a raster that only contains the data from step 1, if the relative height of a raster cell is higher than 3 m, otherwise the raster cell is set to no data. Please note that r[0] represents the first input raster, r[1] the second one. Fig. 3 shows the result in red (superimposed to the DTM in a GIS).

Figure 3: Normalized DSM with underlain DTM.

Smoothing of raster data & filling gaps

When processing point cloud data into a regular raster, gaps may occur due to occlusions in the original point cloud (mainly an issue with terrestrial laser scanners) or due to variations in point density that affect the interpolation (e.g. grid size smaller than point spacing). However, for many geomorphological applications, especially when applying distributed models that rely on a raster input for the base topography, the DTM should be gap-free.

OPALS provides a possibility to remove gaps by smoothing the raster data. This is done by averaging the height in a certain local neighbourhood (kernel). Therefore, we use the opalsStatFilter module. The aim of Module StatFilter is to perform statistical filtering of raster datasets with respect to a user-defined kernel area. However, filtering an entire raster surface would lead to alterations that are undesirable for all non-empty raster cells. Consequently, we want to fill the gaps in the DTM by taking the information from the smoothed DTM while not changing all other cell values at all. This requires three successive steps:

  1. Apply opalsStatFilter on the previously generated DTM (dsm_3m_z_min.tif). A mandatory parameter is the statistical feature. In our case we used the mean of all pixels within our defined kernel size. Please note that the kernel size is to be defined as a radius (in pixels) rather than as number of columns and lines.
    opalsStatFilter -inFile dsm_3m_z_min.tif -outFile dsm_3m_smooth.tif -feature mean -kernelSize 3
  2. Use a conditional statement with the Algebra module to only affect the raster cells that contain gaps. All other raster cells should remain unaltered.
    opalsAlgebra -inFile dsm_3m_z_min.tif -inFile dsm_3m_smooth.tif -outFile dsm_3m_filled.tif -formula "r[0] ? r[0] : r[1]" -gridSize 3.0
  3. Finally, visualize the corrected DTM as a hillshade (Fig. 4).
    opalsShade -inFile dsm_3m_filled.tif -outFile dsm_3m_filled_shade.tif -pixelsize 3.0
Figure 4: Comparison between uncorrected (left) and corrected gap-free DSM subset (centre). Hillshade with corrected data (right).

Derivation of morphometric parameters

Although digital elevation models are the most commonly used LiDAR product, morphometric derivatives (such as slope, aspect, surface roughness) greatly improve the topographic interpretation for natural hazard management. Here some simple steps are shown how to produce such morphometric parameters.

Slope and aspect

Slope and aspect (the compass direction that a slope faces) are to commonly used parameters in geomorphological analyses. Slope is usually the most significant and sensitive parameter in physically based modelling, while aspect can be an evenly important parameter, depending on the research question (e.g. when micro-climatic influences are of importance). The opalsGrid module offers integrated features that easily generate such raster maps (Fig. 5). Those features are slope and expos.

  1. Apply opalsGrid to the basic mtrainier_points.laz file after importing it.
    opalsImport -inFile mtrainier_points.laz -outFile mtrainier_points.odm
    opalsGrid -inFile mtrainier_points.odm -outFile mtrainier_points.tif -gridSize 10 -interpolation movingPlanes -feature sigmaZ slope expos
  2. If only ground points should be considered, a filter can be directly appended to the command line. Filters are specified using their names. For filters depending on parameters, respective filter parameter definitions must be appended, enclosed in square brackets. These must adhere to filter-specific parameter syntax (more information on filters can be found here).
    opalsGrid -inFile mtrainier_points.odm -outFile mtrainier_points_ground.tif -gridSize 10 -interpolation movingPlanes -feature slope expos -filter "Class[Ground]"
  3. A common task in any GIS program (such as QGIS or ArcGIS) is the classification separating the z range of the grid model into distinct classes and assigning a specific color to each class. OPALS offers a built-in module (opalsZColor) to derive color coded visualizations of grid models as geo-coded raster images. The optional parameter palFile contains the definition of a background color.
    opalsZColor -inFile mtrainier_points_expos.tif -palFile aspectPal.xml
    opalsZColor -inFile mtrainier_points_slope.tif -palFile slopePal.xml
Figure 5: Color coded aspect map (left) and slope map (right).

Surface roughness and openness

Surface roughness is an important geomorphological parameter as it heavily influences any flow processes (such as lahars, debris flows or river floods). In OPALS, the roughness (sigmaZ) is expressed as the standard deviation of the interpolated grid height (with movingAverage or movingPlanes interpolator only). Same as with the slope and aspect, it is contained as a feature in the opalsGrid module. Another measure of convexity/concavity is contained in the module opalsOpenness (Fig. 6).

  1. Apply Module Grid to the imported mtrainier_points.odm file and at the same time also apply a region and a class filter, as we want to look at a higher resolution subset of the entire extent of the study area.
    opalsGrid -infile mtrainier_points.odm -outfile mtrainier_point_sub.tif -gridSize 1 -interpolation movingPlanes -feature sigmaZ -filter "Region[589121.69 5185802.3 590380.48 5186766.81] AND Class[Ground]"
  2. The Openness module provides a raster map of local viewsheds (i.e. openness) based on a grid DTM.
    opalsOpenness -inFile mtrainier_points.tif
Figure 6: Openness raster map with the surface roughness subset overlain.
Author
Ekrem Canli
Date
January 2016

References

  • Abellán, A., Calvet, J., Vilaplana, J.M., Blanchard, J., 2010. Detection and spatial prediction of rockfalls by means of terrestrial laser scanner monitoring. Geomorphology 119, 162–171. doi:10.1016/j.geomorph.2010.03.016
  • Abellán, A., Jaboyedoff, M., Oppikofer, T., Vilaplana, J.M., 2009. Detection of millimetric deformation using a terrestrial laser scanner: experiment and application to a rockfall event. Natural Hazards and Earth System Science 9, 365–372.
  • Höfle, B., Rutzinger, M., 2011. Topographic airborne LiDAR in geomorphology: A technological perspective. Zeitschrift für Geomorphologie, Supplementary Issues 55, 1–29. doi:10.1127/0372-8854/2011/0055S2-0043
  • Jaboyedoff, M., Oppikofer, T., Abellán, A., Derron, M.-H., Loye, A., Metzger, R., Pedrazzini, A., 2012. Use of LIDAR in landslide investigations: a review. Natural Hazards 61, 5–28. doi:10.1007/s11069-010-9634-2
  • Lague, D., Brodu, N., Leroux, J., 2013. Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z). ISPRS Journal of Photogrammetry and Remote Sensing 82, 10–26. doi:10.1016/j.isprsjprs.2013.04.009
  • Milan, D.J., Heritage, G.L., Hetherington, D., 2007. Application of a 3D laser scanner in the assessment of erosion and deposition volumes and channel change in a proglacial river. Earth Surface Processes and Landforms 32, 1657–1674. doi:10.1002/esp.1592
  • O’Neal, M.A., Pizzuto, J.E., 2011. The rates and spatial patterns of annual riverbank erosion revealed through terrestrial laser-scanner surveys of the South River, Virginia. Earth Surface Processes and Landforms 36, 695–701. doi:10.1002/esp.2098
  • Pfeifer, N., Mandlburger, G., Otepka, J., Karel, W., 2014. OPALS – A framework for Airborne Laser Scanning data analysis. Computers, Environment and Urban Systems 45, 125–136. doi:10.1016/j.compenvurbsys.2013.11.002
  • Pinsker, L., 2004. Paths of Destruction: The Hidden Threat at Mount Rainier. Geotimes, Apr. 2004, online: http://www.geotimes.org/apr04/feature_MountRainier.html (accessed: Jan 13, 2016)
  • Rosser, N.J., Petley, D.N., Lim, M., Dunning, S.A., Allison, R.J., 2005. Terrestrial laser scanning for monitoring the process of hard rock coastal cliff erosion. Quarterly Journal of Engineering Geology and Hydrogeology 38, 363–375.
  • Schürch, P., Densmore, A.L., Rosser, N.J., Lim, M., McArdell, B.W., 2011. Detection of surface change in complex topography using terrestrial laser scanning: application to the Illgraben debris-flow channel.. Earth Surface Processes and Landforms 36, 1847–1859. doi:10.1002/esp.2206
  • Teza, G., Pesci, A., Genevois, R., Galgaro, A., 2008. Characterization of landslide ground surface kinematics from terrestrial laser scanning and strain field computation. Geomorphology 97, 424–437. doi:10.1016/j.geomorph.2007.09.003
opalsZColor is the executable file of Module ZColor
Definition: ModuleExecutables.hpp:253
opalsOpenness is the executable file of Module Openness
Definition: ModuleExecutables.hpp:153
@ movingPlanes
moving (tilted) plane interpolation
opalsAlgebra is the executable file of Module Algebra
Definition: ModuleExecutables.hpp:13
opalsImport is the executable file of Module Import
Definition: ModuleExecutables.hpp:113
@ slope
steepest slope [%] (deprecated: use slpPerc instead)
Definition: GridFeature.hpp:16
@ mean
Mean.
opalsStatFilter is the executable file of Module StatFilter
Definition: ModuleExecutables.hpp:213
@ kernelSize
output grid size equals entire kernel size
opalsGrid is the executable file of Module Grid
Definition: ModuleExecutables.hpp:93
opalsCell is the executable file of Module Cell
Definition: ModuleExecutables.hpp:23
@ filter
string to be parsed in construction of DM::IFilter (various modules)
@ feature
Use a statistic feature of the boundary gap points for filling.
@ formula
formula string for albegraic grid computations (opalsAlgebra)
@ palFile
palette file (opalsZcolor)
opalsShade is the executable file of Module Shade
Definition: ModuleExecutables.hpp:198