DOI: 10.5281/zenodo.17243697
title: Test DTM for "Caldera Blanca" (Lanzarote Island), used in the submitted paper on surface roughness analysis in geomorphometry

description: The uploaded DTM “Lanzarote.tif” is used as testing surface for the following paper:

Trevisani, S., & Guth, P. L. (2025). Surface Roughness in Geomorphometry: From Basic Metrics Toward a Coherent Framework. Remote Sensing, 17(23), 3864. 

https://doi.org/10.3390/rs17233864

 The DTM has been derived by means of interpolation of lidar point cloud available at: https://centrodedescargas.cnig.es/CentroDescargas/lidar-tercera-cobertura. At the cited portal you will find the tiles of the LiDAR 3D point clouds with national coverage and colored with 4 bands (RGBI), corresponding to the third topographic lidar survey for Spain (2022-2025). From the same lidar point cloud RGB imagery has been derived at 1 m resolution ("lanzaroteRGB.tif").

The shapefiles included represent areas considered in the submitted paper to perform the analysis and design figures. The statistical correlations between different surface roughness indices have been computed for the spatial domain represented by "AOIbasics.shp". 



Some characteristics of the lidar point cloud tiles derived from the web portal:

GRS: ETRS89 in the mainland Spain, Balearic Islands, Ceuta and Melilla, and REGCAN95 in the Canary Islands (both systems compatible with WGS84). UTM projection in the corresponding zone. Orthometric heights.




Format: LAZ (LAS compression file format) file.

Average point density: 5 points/m² or greater


The 2 m DTM has been derived interpolating the ground points using the function "LAS data set to raster" of ArcGIS Pro Esri (version 3.4.0); parameters used: binning, average and void filling via natural neighbor interpolation.  

 

Given that the point clouds have been downloaded from the website of the General Directorate of the National Geographic Institute of Spain and are covered by CC-BY 4.0 license, the derived digital terrain model  and RGB imagery are a “Obra derivada de LiDAR-PNOA-cob3 2022-2025 CC-BY 4.0 scne.es”

 

 

Below for your convenience some of the R key commands (commented) used for the paper with this dataset

 
description: #February 2026

It has been published the new version of SurfRough in R Cran (v1.2). These are the updates:

# SurfRough 0.0.1.2 (January 2026)
## New Features
* Added new functions and kernels with Rcpp implementations:
  - `TRIbi`: Like TRI but using bilinear interpolation along diagonals
  - `RRIcore`: Using only the inner four second order directional differences (i.e., inner part of RRI kernel)
  - `RRIk4`: RRI with fourth order directional differences (filtering smooth curvature)
  - `k1ck4`: Kernels lag 1 pixel and fourth order differences (for filtering curvature) to be used with geostatistical estimators (`Madscan` and `Meanscan`)
* `MadScanL`: A version of Madscan designed to work with large files, with better management of disk space (C++ version planned for future)
* Added functionality for computing roughness indices based on measures of dispersion of local surface parameters (e.g., slope, profile curvature, residual surface, etc.)
* Two fast functions using terra's `focalcpp` have been implemented:
  - `iqrST`: Fast C++ version of IQR for computing focal statistics (method type 7)
  - `stdST`: Fast C++ version of standard deviation for computing focal statistics (population variance, divided by n)

###############

 

#November 2025
#Coded provided without any warranty!
#Basic R commands to replicate the main analyses performed in the paper
#For the computation of basic surface roughness (SR) indices (and some more complex one)
#These are the required libraries. Others will be considered later#I noticed that with the newer version of terra the functions using focal() #function (e.g. MadScan) there could be issues with large files, #likely related to problems in writing files in the temporary folder. In my pc, apparently, #these settings resolved the issues (subtitute the tempdir with your own directory):#terraOptions(memfrac=0.9,tempdir="G:/RTemp",todisk=T)
library(SurfRough)
library(viridis)
library(MultiscaleDTM)

#Plotting utilities for terra
#Plotting related
hillshade=function(alt,azimuth=45,direction=315,...){
  slope=terrain(alt, "slope", unit="radians")
  aspect=terrain(alt, "aspect", unit="radians")
  shade(slope, aspect, azimuth, direction)
}
plotShade=function(dem,...){plot(hillshade(dem),col=grey(0:100/100),legend=F,...)
}
#End plotting utilities


#Set you folder with the DEM dowloaded from Zenodo including also the shapefiles
#of the AOIs to cover the same extent of the paper and figures

#load the DEM
dem=rast("lanzarote.tif")
#this covers a larger extent with respect to one corresponding to the presented computations
plotShade(dem)
plot(dem,col=viridis(100),alpha=0.65,add=T)

#Crop to get the same extent
aoi1=vect("aoiBasics.shp")
dem=crop(dem, aoi1)

#plot
plotShade(dem)
plot(dem,col=viridis(100),alpha=0.65,add=T)

#Compute SR indices
#ris will contain SR indices
ris=dem

#Radial Roughness index (using DDs of order 2)
ris$RRI=RRI(dem)

#plot
plotShade(dem)
plot(ris$RRI,breaks=10,breakby="cases",alpha=0.65,add=T)

#Dispersion of normal vectors to surface (aka VRM)
#you can use both Surfrough as well as MultiscaleDTM
#packages. MultiscaleDTM uses only square windows for VRM.
#surfRough can use custom windows (e.g., both circular as well as
#rectangular). Using the same window, MultiscaleDTM and SurfRough functions
#provide the same result.#SurfRough function
#NVD with circular windows
w=KernelCircular(3)
ris$NVD=circularDispersionNV(dem,w)
gc()

#plot
plotShade(dem)
plot(ris$NVD,breaks=10,breakby="cases",alpha=0.65,add=T)

#If you like you can check that this
#function provides the same results of VRM
#in MultiscaleDTM package#MultiscaleDTM VRM function
NVD1=VRM(dem,w=c(5,5))
w=KernelRectangular(5,5)
NVD2=circularDispersionNV(dem,w)
plot(NVD1-NVD2)
rm(NVD1,NVD2,w)
gc()

#use MultiscaleDTM for surface area method
ris$SA=SAPA(dem, w=c(5,5), slope_correction = TRUE)
gc()
#plot
plotShade(dem)
plot(ris$SA,breaks=10,breakby="cases",alpha=0.65,add=T)


#Another approach is to compute the dispersion/variability of local surface
#parameters (LSPs) such as slope and residual surface

#The standard deviation of residual surface
#There are different methods (many) to derive a residual surface
#this from MUltiscaleDTM fits a planar surface for smoothing
ris$STDResAdjsd=AdjSD(dem, w=c(5,5), na.rm = TRUE)
gc()

#Computing the residual via moving means
#define the kernel
w=KernelCircular(3)
#compute the mean
meanDEM=focal(dem, w,fun=mean,expand=F,na.rm=F)
#compute the residual
res=dem-meanDEM
#NB the following, as is, it is not particularly fast;
#it can improved a lot with some coding using terra focalCpp function
#You can compute in any GIS permitting raster algebra and focal statistics
#the same is true for IQRslope and IQRres.
#We will implement the faster versions on the next updates of SurfRough package#If you need now to compute things fast use the new functions (much faster).
#ris$stdRes=sqrt(focal(res, w,fun=var,expand=F,na.rm=F))
#ris$iqrRes=(focal(res, w,fun=IQR,expand=F,na.rm=T))#use the new versions#You can consider other LSPs such as slope, profile curvature, etc.ris$stdRes=stdST(res,w=w)ris$iqrRes=iqrST(res,w=w)
rm(res,meanDEM)

#plot
plotShade(dem)
plot(ris$STDResAdjsd,breaks=10,breakby="cases",alpha=0.65,add=T)
plotShade(dem)
plot(ris$stdRes,breaks=10,breakby="cases",alpha=0.65,add=T)
plotShade(dem)
plot(ris$iqrRes,breaks=10,breakby="cases",alpha=0.65,add=T)
par(mfrow=c(1,1))

#standard deviation of Slope
#slope in %
ris$slope=100*tan(terrain(dem, "slope", unit="radians"))
#std
w=KernelCircular(3)
#ris$STDSlope=sqrt(focal(ris$slope, w,fun=var,expand=F,na.rm=F))ris$STDSlope=stdST(ris$slope,w=w)
gc()

#plot
plotShade(dem)
plot(ris$STDSlope,breaks=10,breakby="cases",alpha=0.65,add=T)
#for robust use IQR

#As simplified geostatistical approach you can use MAD with second order DDs
#you can compute lags 0.5, 1 and 2 pixels respectively with the kernels
#k05ck2, k1ck2, k2ck2
#Here an example with a lag of 1 pixel and a circular window with radius of 3 pixels
#Now we consider the simplified ones
w=KernelCircular(3)
#lag1
MAD1ck2=Madscan(dem,k1ck2, w)
gc()
#it provides 3 values: omnidirectional SR, anisotropy strength and direction (direction of max #continuity)
#Here we do not consider direction but only strength
ris$OmniMAD1ck2=MAD1ck2$IsoRough
ris$AnisoMAD1ck2=MAD1ck2$AnisoR
rm(MAD1ck2)
gc()

#plot
plotShade(dem)
plot(ris$OmniMAD1ck2,breaks=10,breakby="cases",alpha=0.65,add=T)
plotShade(dem)
plot(ris$AnisoMAD1ck2,breaks=10,breakby="cases",alpha=0.65,add=T)

#Having introduced SR anisotropy, other aspects of ST can be describe by TT approach
library(Rsagacmd)
#search_saga()
# initiate a saga object. You need to set the location of Saga binaries
saga <- saga_gis(saga_bin="C:/Program Files/SAGA/saga_cmd.exe",raster_backend = "terra")
#saga_version("C:/Program Files/SAGA/saga_cmd.exe")
ris$stIwaT1cm=saga$ta_morphometry$terrain_surface_texture(dem,scale=5,epsilon =0.01,method =0)
ris$stIwaT5cm=saga$ta_morphometry$terrain_surface_texture(dem,scale=5,epsilon =0.05,method =0)
ris$stIwaT10cm=saga$ta_morphometry$terrain_surface_texture(dem,scale=5,epsilon =0.1,method =0)
ris$stIwaT15cm=saga$ta_morphometry$terrain_surface_texture(dem,scale=5,epsilon =0.15,method =0)

plot(ris$stIwaT1cm)
plot(ris$stIwaT5cm)
plot(ris$stIwaT10cm)
plot(ris$stIwaT15cm)


#some stats
library(ggplot2)
#install.packages("ggcorrplot")
library(ggcorrplot)
summary(ris)
#correlation
forCor=na.omit(as.data.frame(ris))
corr=round(cor(forCor),2)
ggcorrplot(corr[nrow(corr):1,], hc.order = FALSE, method = "square", type = "full", lab = T,
           digits = 4, tl.srt = 45, lab_size = 4, tl.cex = 10) +
  scale_fill_gradientn(colors = c("#2166ac", "#92c5de"  , "white", "#f4a582", "#b2182b"),
                       values = scales::rescale(c(-1, -0.5, 0, 0.5, 1)),  # Control positions
                       limits = c(-1, 1))

#for Haralick GLCM
#install.packages("GLCMTextures")
library(GLCMTextures) #Load GLCMTextures package
#Compute on residual surface
#It is not fast.
w=KernelCircular(3)
meanDEM=focal(dem, w,fun=mean,expand=F,na.rm=F)
res=dem-meanDEM
rm(meanDEM)
ht=glcm_textures(res, w = c(5,5), n_levels = 16, quant_method = "prob", shift = list(c(1, 0), c(1, 1), c(0, 1), c(-1, 1))) 

#You like you can compute MAD on residuals and consider longer lags than MADk2
#So you can see analogies and differences between geostatistical approaches and GLCM
#MAD is fast
#compute MAD indices for different lag distances
w=KernelCircular(3)
MAD1c=Madscan(res,k1c, w)
MAD2c=Madscan(res,k2c, w)
MAD4c=Madscan(res,k4c, w)
ht$MAD1cIso=MAD1c$IsoRough
ht$MAD1cAniso=MAD1c$AnisoR
#Relative roughness
ht$MADRelRough=MAD1c$IsoRough/((MAD1c$IsoRough+MAD2c$IsoRough+MAD4c$IsoRough+0.0001)/3)
#ht$anisoMean=(MAD1c$AnisoR+MAD2c$AnisoR+MAD4c$AnisoR)/3

#
summary(ht)
#correlation
forCor=na.omit(as.data.frame(ht))
corr=round(cor(forCor),2)
ggcorrplot(corr[nrow(corr):1,], hc.order = FALSE, method = "square", type = "full", lab = T,
           digits = 4, tl.srt = 45, lab_size = 4, tl.cex = 10) +
  scale_fill_gradientn(colors = c("#2166ac", "#92c5de"  , "white", "#f4a582", "#b2182b"),
                       values = scales::rescale(c(-1, -0.5, 0, 0.5, 1)),  # Control positions
                       limits = c(-1, 1))