Alexander Ilich September 22, 2022




Please cite as

Ilich, Alexander R.; Misiuk, Benjamin; Lecours, Vincent; Murawski, Steven A.; 2021. “MultiscaleDTM”,


This package calculates multi-scale geomorphometric terrain attributes from regularly gridded digital terrain models (DTM; i.e. elevation or bathymetry rasters) via a specified window size.

Figure adapted from Wilson et al. (2007)

Install and Load Package

If you don’t already have remotes installed, use the code install.packages("remotes")

Then to install this package use the code remotes::install_github("ailich/MultiscaleDTM") (If you are using Windows, you may need to install Rtools using the instructions found here).

This package relies on the terra package for handling of spatial raster data.

Main Functions

Slope, Aspect and Curvature

Z = aX^2 + bY^2 +cXY+ dX +eY +f

Figure adapted from Walbridge et al., (2018)


Figure adapted from Sappington et al. (2007)

Figure adapted from Habib (2021)

 = 1-
x = sin()*cos()

Figure adapted from Friedman et al. (2012) and created with

Figure adapted from Jenness (2004)

Figure adapted from Cavalli et al. (2008)

Relative Position

Figure adapted from Lundblad et al., (2006)


In this tutorial we will calculate various terrain attributes using a 5 x 5 cell rectangular window. Any rectangular odd numbered window size however could be used (see figure directly below). Window sizes are specified with a vector of length 2 of c(n_rows, n_cols). If a single number is provided it will be used for both the number of rows and columns. The only metric that does not follow this syntax is BPI which uses an annulus shaped focal window which we will calculate using an inner radius of 4 and an outer radius of 6 cells.

Load packages

library(MultiscaleDTM) #Load MultiscaleDTM package

See package help page


Read in Data

r<- rast(volcano, extent= ext(2667400, 2667400 + ncol(volcano)*10, 6478700, 6478700 + nrow(volcano)*10), crs = "EPSG:27200")

Slope, Aspect, and Curvature

slp_asp<- SlpAsp(r = r, w = c(5,5), unit = "degrees", method = "queen", metrics = c("slope", "aspect", "eastness", "northness"))

qmetrics<- Qfit(r, w = c(5,5), unit = "degrees", metrics = c("elev", "qslope", "qaspect", "qeastness", "qnorthness", "profc", "planc", "twistc", "meanc", "maxc", "minc", "features"), na.rm = TRUE)

To explore these measures in an interactive environment use explore_terrain() or go to this website


vrm<- VRM(r, w=c(5,5), na.rm = TRUE)

Note: multi-scale SAPA is experimental. The established metric by De Preez (2015) would use w=1.

sapa<- SAPA(r, w=c(5,5), slope_correction = TRUE)

adj_SD<- AdjSD(r, w=c(5,5), na.rm = TRUE)

rie<- RIE(r, w=c(5,5), na.rm = TRUE)

Relative Position

tpi<- TPI(r, w=c(5,5), na.rm = TRUE)

rdmv<- RDMV(r, w=c(5,5), na.rm = TRUE, method="range")

bpi<- BPI(r, radius = c(4,6), unit = "cell", na.rm = TRUE)

The annulus window for BPI can be specified in either cell units (number of raster cells) or in map units (e.g. meters) which can be useful if your x and y resolutions are not equal. Additionally, the function annulus_window can be used to verify that you are specifying your window correctly (NA’s are excluded cells and 1’s are included cells) and can be directly supplied to the w argument in the BPI funtion instead of using radius and unit arguments.

annulus_window(radius = c(4,6), unit = "cell")
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]   NA   NA   NA   NA   NA   NA    1   NA   NA    NA    NA    NA    NA
##  [2,]   NA   NA   NA    1    1    1    1    1    1     1    NA    NA    NA
##  [3,]   NA   NA    1    1    1    1    1    1    1     1     1    NA    NA
##  [4,]   NA    1    1    1   NA   NA   NA   NA   NA     1     1     1    NA
##  [5,]   NA    1    1   NA   NA   NA   NA   NA   NA    NA     1     1    NA
##  [6,]   NA    1    1   NA   NA   NA   NA   NA   NA    NA     1     1    NA
##  [7,]    1    1    1   NA   NA   NA   NA   NA   NA    NA     1     1     1
##  [8,]   NA    1    1   NA   NA   NA   NA   NA   NA    NA     1     1    NA
##  [9,]   NA    1    1   NA   NA   NA   NA   NA   NA    NA     1     1    NA
## [10,]   NA    1    1    1   NA   NA   NA   NA   NA     1     1     1    NA
## [11,]   NA   NA    1    1    1    1    1    1    1     1     1    NA    NA
## [12,]   NA   NA   NA    1    1    1    1    1    1     1    NA    NA    NA
## [13,]   NA   NA   NA   NA   NA   NA    1   NA   NA    NA    NA    NA    NA


Cavalli, M., Tarolli, P., Marchi, L., Dalla Fontana, G., 2008. The effectiveness of airborne LiDAR data in the recognition of channel-bed morphology. CATENA 73, 249–260.

Du Preez, C., 2015. A new arc–chord ratio (ACR) rugosity index for quantifying three-dimensional landscape structural complexity. Landscape Ecol 30, 181–192.

Evans, I.S., 1980. An integrated system of terrain analysis and slope mapping. Zeitschrift f¨ur Geomorphologic Suppl-Bd 36, 274–295.

Fleming, M.D., Hoffer, R.M., 1979. Machine processing of landsat MSS data and DMA topographic data for forest cover type mapping (No. LARS Technical Report 062879). Laboratory for Applications of Remote Sensing, Purdue University, West Lafayette, Indiana.

Friedman, A., Pizarro, O., Williams, S.B., Johnson-Roberson, M., 2012. Multi-Scale Measures of Rugosity, Slope and Aspect from Benthic Stereo Image Reconstructions. PLOS ONE 7, e50440.

Habib, M., 2021. Quantifying Topographic Ruggedness Using Principal Component Analysis. Advances in Civil Engineering 2021, e3311912.

Horn, B.K., 1981. Hill Shading and the Reflectance Map. Proceedings of the IEEE 69, 14–47.

Jenness, J.S., 2004. Calculating landscape surface area from digital elevation models. Wildlife Society Bulletin 32, 829–839.;2

Lecours, V., Devillers, R., Simms, A.E., Lucieer, V.L., Brown, C.J., 2017. Towards a Framework for Terrain Attribute Selection in Environmental Studies. Environmental Modelling & Software 89, 19–30.

Lundblad, E.R., Wright, D.J., Miller, J., Larkin, E.M., Rinehart, R., Naar, D.F., Donahue, B.T., Anderson, S.M., Battista, T., 2006. A benthic terrain classification scheme for American Samoa. Marine Geodesy 29, 89–111.

Minár, J., Evans, I.S., Jenčo, M., 2020. A comprehensive system of definitions of land surface (topographic) curvatures, with implications for their application in geoscience modelling and prediction. Earth-Science Reviews 211, 103414.

Misiuk, B., Lecours, V., Dolan, M.F.J., Robert, K., 2021. Evaluating the Suitability of Multi-Scale Terrain Attribute Calculation Approaches for Seabed Mapping Applications. Marine Geodesy 44, 327–385.

Risk, M.J., 1972. Fish Diversity on a Coral Reef in the Virgin Islands. Atoll Research Bulletin 153, 1–4.

Ritter, P., 1987. A vector-based slope and aspect generation algorithm. Photogrammetric Engineering and Remote Sensing 53, 1109–1111.

Sappington, J.M., Longshore, K.M., Thompson, D.B., 2007. Quantifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert. The Journal of Wildlife Management 71, 1419–1426.

Walbridge, S., Slocum, N., Pobuda, M., Wright, D.J., 2018. Unified geomorphological analysis workflows with benthic terrain modeler. Geosciences 8, 94.

Weiss, A., 2001. Topographic Position and Landforms Analysis. Presented at the ESRI user conference, San Diego, CA.

Wilson, M.F., O’Connell, B., Brown, C., Guinan, J.C., Grehan, A.J., 2007. Multiscale Terrain Analysis of Multibeam Bathymetry Data for Habitat Mapping on the Continental Slope. Marine Geodesy 30, 3-35.

Wood, J., 1996. The geomorphological characterisation of digital elevation models (Ph.D.). University of Leicester.