From a 1m DSM in the Tiraumea catchment, New Zealand, terrain derivatives will be computed with the following code. The aim is to detect Earthflows, possibly through machine learning, where a model would be fed with the terrain derivatives, and optical information.
The study area is shown below:
# Call mapping library and set options
library(mapview)
mapviewOptions(basemaps = 'OpenStreetMap', console = F, verbose = F, leafletWidth = 800)library(sf, quietly = T, verbose = F, warn.conflicts = F)
tiraumea = st_read('study_area/Tiraumea.shp', quiet = T)
mapview(tiraumea)Earthflows have distinct characteristics that define them, which can somehow be inferred from topographic characteristics. A set of topographic variables have been selected for calculation on the original 1m DSM and on a filter version of it, to 3m. The variables computed are listed below:
list = read.csv('code_list.csv', sep = ';')
library(dplyr, quietly = T)
library(DT)
list %>% arrange(desc(Module), Derivative_Name) %>% datatable()The naming convention will include the stated name above and for the layers computed based on the 3m filtered window, they will have the suffix 3x3.
In order to make the RSAGA workflow work, it is first needed to convert the original DSM into a SAGA compatible format (.sgrd)
dsm = raster('data_rs/2016_HRC_DSM_tiraumea.tif')
writeRaster(dsm, 'terrain/input_dsm/dsm.sgrd', format = 'SAGA', overwrite = T, NAflag = 0, prj = T)Then, we prepare the RSAGA environment. Currently the SAGA files are on the same project, so there is no need to change the paths. Check that SAGA versions are always compatible when transferring to other computers.
library(RSAGA)
env = rsaga.env(
path = "software\\saga-7.6.1_x64",
modules = "software\\saga-7.6.1_x64\\tools"
# cmd = "software\\saga-7.6.1_x64\\saga_cmd.exe" # For some reason this throws an error, but when commented out it still works.
)Verify specified path to SAGA command line program...
Verify specified path to SAGA modules...
DoneInitially, we will do a simple filter of the original DSM with a 3x3 window size, to then be able to compute derivatives on all the values.
rsaga.filter.simple(
in.grid = 'terrain/input_dsm/dsm.sgrd',
out.grid = 'terrain/input_dsm/dsm_3x3.sgrd',
method = 'smooth', radius = 1, mode = 'square',
env = env
)Below, all the used derivatives will be explained with code, and when needed with additional information. Everything will be computed both for the original DSM and for the filtered 3x3 DSM.
For the Topographic Position Index (tpidx) at a filtered stage, the method required a large computational power. Hence, we resampled the DSM to 3m and calculated it based on that.
dsm = raster('data_rs/2016_HRC_DSM_tiraumea.tif')
rgbi = stack('data_rs/HRC_2016_RGBI_mosaic_tiraumea.kea')
crs = crs(rgbi)
crs(dsm) = crs
NAvalue(dsm) = -9999
dsm3 = aggregate(dsm, fact = 3)
writeRaster(dsm3, 'data_rs/2016_HRC_DSM_tiraumea_3m.tif')
writeRaster(
dsm3, 'terrain/input_dsm/dsm_3m.sgrd', format = 'SAGA', overwrite = T, NAflag = 0, prj = T
)With this I will compute overall 9 outputs:
rsaga.slope.asp.curv(
in.dem = "terrain/input_dsm/dsm.sgrd",
out.aspect = "terrain/out_products/aspect.sgrd",
out.slope = "terrain/out_products/slope.sgrd",
out.cgene = "terrain/out_products/cgene.sgrd",
out.cplan = "terrain/out_products/cplan.sgrd",
out.cprof = "terrain/out_products/cprof.sgrd",
out.ccros = "terrain/out_products/ccros.sgrd",
out.clong = "terrain/out_products/clong.sgrd",
out.cmaxi = "terrain/out_products/cmaxi.sgrd",
out.cmini = "terrain/out_products/cmini.sgrd",
unit.slope = "degrees", unit.aspect = "degrees",
env = env,
method = "poly2zevenbergen"
)From these modules on, we need to find out the usage first, because there are no shortcuts for them. We always refer to this table for the morphometry module to get the respective code:
rsaga.get.modules('ta_morphometry', env = env) %>% knitr::kable()
|
rsaga.get.usage('ta_morphometry', 14, env = env)This module will get 5 outputs, and will work with default W, T, E.
rsaga.geoprocessor(
'ta_morphometry', 14,
list(
DEM='terrain/input_dsm/dsm.sgrd',
HO='terrain/out_products/slhgt.sgrd',
HU='terrain/out_products/vldpt.sgrd',
NH='terrain/out_products/nrhgt.sgrd',
SH='terrain/out_products/sthgt.sgrd',
MS='terrain/out_products/mdslp.sgrd'
),
env = env
)It describes wet areas by assuming water accumulation in flat areas is due to upslope, local and downslope topography. It is a quantitative estimation of the hydraulic gradient. Obtained by calculating the downhill distance when water loses a determined quantity of energy from precipitation.
This tool has its own module, and will return one outcome, coded ‘ddgrd’. We get the usage:
rsaga.get.usage('ta_morphometry', 9, env = env)And call it, with default distance (10) and gradient in degrees:
rsaga.geoprocessor(
'ta_morphometry', 9,
list(
DEM='terrain/input_dsm/dsm.sgrd',
GRADIENT='terrain/out_products/ddgrd.sgrd',
OUTPUT = 2
),
env = env
)We assume that negative MBI values represent areas of net deposition such as depressions and floodplains; positive MBI values represent areas of net erosion such as hillslopes, and MBI values close to zero indicate areas where there is a balance between erosion and deposition such as low slopes and plain areas. High positive MBI values occur at convex terrain forms, like upper slopes and crests, while lower MBI values are associated with valley areas and concave zones at lower slopes. Balanced MBI values close to zero can be found in midslope zones and mean a location of no net loss or net accumulation of material.
This tool has its own module, and will return one outcome, coded ‘mbidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 10, env = env)And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 10,
list(
DEM='terrain/input_dsm/dsm.sgrd',
MBI='terrain/out_products/mbidx.sgrd'
),
env = env
)TPI measures the relative topographic position of the central point as the difference between the elevation at this point and the mean elevation within a predetermined neighbourhood. Using TPI, landscapes can be classified in slope position classes. TPI is only one of a vast array of morphometric properties based on neighbouring areas that can be useful in topographic and DEM analysis. Used for roughness determination. The lower the numbers are the lower areas in the landscape. The higher numbers are the higher areas in the landscape. 1
This tool has its own module, and will return one outcome, coded ‘tpidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 18, env = env)And compute with all the defaults, min and max are definitely needed and vary depending on the raster resolution. If not used properly computation will be exhaustive. Hence, for 1 meter resolution we use 20 and for 3 m, 100 as maximum radious
rsaga.geoprocessor(
'ta_morphometry', 18,
list(
DEM='terrain/input_dsm/dsm.sgrd',
TPI='terrain/out_products/tpidx.sgrd',
RADIUS_MIN=0,
RADIUS_MAX=20
),
env = env
)
rsaga.geoprocessor(
'ta_morphometry', 18,
list(
DEM="terrain/input_dsm/dsm_3m.sgrd",
TPI='terrain/out_products/tpidx_3m.sgrd',
RADIUS_MIN=0,
RADIUS_MAX=100
),
env = env
)This module calculates an index of convergence/divergence regarding to overland flow. By its meaning it is similar to plan or horizontal curvature, but gives much smoother results. The calculation uses the aspects of surrounding cells, i.e. it looks to which degree surrounding cells point to the center cell. The result is given as percentages, negative values correspond to convergent, positive to divergent flow conditions. Minus 100 would be like a peak of a cone (a), plus 100 a pit (c), and 0 an even slope (b) 2.
This tool has its own module, and will return one outcome, coded ‘cvidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 1, env = env)And compute with method gradient and neighbours 3x3:
rsaga.geoprocessor(
'ta_morphometry', 1,
list(
ELEVATION='terrain/input_dsm/dsm.sgrd',
RESULT='terrain/out_products/cvidx.sgrd',
METHOD=1,
NEIGHBOURS=1
),
env = env
)Quantitative measure of topographic heterogeneity by calculating the sum change in elevation between a grid cell and its eight neighbor grid cells. This tool works with absolute values by squaring the differences between the target and neighbor cells, then taking the square root. Concave and convex shape areas could have similar values. The value of this metric will vary as a function of the size and complexity of the terrain used in the analysis. The closer you are to 0 the less rugged the terrain likely is. The bigger the number is, eg 105, then the terrain is likely to be more rugged 3
This tool has its own module, and will return one outcome, coded ‘tridx’. We get the usage:
rsaga.get.usage('ta_morphometry', 16, env = env)And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 16,
list(
DEM='terrain/input_dsm/dsm.sgrd',
TRI='terrain/out_products/tridx.sgrd'
),
env = env
)This parameter emphasizes fine versus coarse expression of topographic spacing, or “grain”… Texture is calculated by extracting grid cells (here, informally, “pits” and “peaks”) that outline the distribution of valleys and ridges. It is defined by both relief (feature frequency) and spacing in the horizontal. Each grid cell value represents the relative frequency (in percent) of the number of pits and peaks within a radius of ten cells (Iwahashi and Pike, 2007. pp.412-413). It should be noted that it is not clear that the relative frequency is actually a percentage. According to Iwahashi and Pike (2007, p.30), To ensure statistically robust classes, thresholds for subdividing the images are arbitrarily set at mean values of frequency distributions of the input variables.
This tool has its own module, and will return one outcome, coded ‘textu’. We get the usage:
rsaga.get.usage('ta_morphometry', 20, env = env)And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 20,
list(
DEM='terrain/input_dsm/dsm.sgrd',
TEXTURE='terrain/out_products/textu.sgrd'
),
env = env
)Again I will access the codes for the tools to get the usage:
rsaga.get.modules('ta_hydrology', env = env) %>% knitr::kable()
|
Simple flow accumulation related index, calculated as difference between maximum and minimum elevation in catchment area divided by square root of catchment area size. The calculation is performed for each grid cell, therefore minimum elevation is same as elevation at cell’s position. Due to the discrete character of a single maximum elevation, flow calculation is simply done with Deterministic 8.
This tool has its own module, and will return one outcome, coded ‘mridx’. We get the usage:
rsaga.get.usage('ta_hydrology', 23, env = env)And compute with all the defaults:
rsaga.geoprocessor(
'ta_hydrology', 23,
list(
DEM='terrain/input_dsm/dsm.sgrd',
MRN='terrain/out_products/mridx.sgrd'
),
env = env
)L is the slope length factor, representing the effect of slope length on erosion. It is the ratio of soil loss from the field slope length to that from a 72.6-foot (22.1-meter) length on the same soil type and gradient. Slope length is the distance from the origin of overland flow along its flow path to the location of either concentrated flow or deposition. S is the slope steepness. Represents the effect of slope steepness on erosion. Soil loss increases more rapidly with slope steepness than it does with slope length. L factor and S factor are usually considered together. LS factors = the slope length factor L computes the effect of slope length on erosion and the slope steepness factor S computes the effect of slope steepness on erosion. Values of both L and S equal 1 for the unit plot conditions of 72.6 ft. length and 9 percent steepness. Values of L and S are relative and represent how erodible the particular slope length and steepness is relative to the 72.6 ft long, 9% steep unit plot. Thus some values of L and S are less than 1 and some values are greater than 1.
This particular tool I had to compute using the GUI for SAGA, because RSAGA does not let me use the LS Factor (one step) library that I need to do it right. In the GUI I only need to give the DEM. It is saved as ‘lsfct’.
This tool has its own module, and will return one outcome, coded ‘lsfct’. We get the usage:
rsaga.get.usage('terrain_analysis','LS Factor (One Step)', env = env)And compute:
rsaga.geoprocessor(
'terrain_analysis',
'LS Factor (One Step)',
list(
DEM='terrain/input_dsm/dsm.sgrd',
LS_FACTOR='terrain/out_products/lsfct.sgrd',
LS_METHOD=0,
PREPROCESSING=2,
MINSLOPE=0.0001
),
env = env
)The ‘SAGA Wetness Index’ is, as the name says, similar to the ‘Topographic Wetness Index’ (TWI), but it is based on a modified catchment area calculation (‘Modified Catchment Area’), which does not think of the flow as very thin film. As result it predicts for cells situated in valley floors with a small vertical distance to a channel a more realistic, higher potential soil moisture compared to the standard TWI calculation.
This tool has its own module, and will return one outcome, coded ‘sagaw’. We get the usage:
rsaga.get.usage('ta_hydrology', 15, env = env)And compute with defaults:
rsaga.geoprocessor(
'ta_hydrology', 15,
list(
DEM='terrain/input_dsm/dsm.sgrd',
TWI='terrain/out_products/sagaw.sgrd'
),
env = env
)A measurement of the distance from the origin of overland flow along its flow path to the location of either concentrated flow or deposition. 4
This tool has its own module, and will return one outcome, coded ‘sllgt’. We get the usage:
rsaga.get.usage('ta_hydrology', 7, env = env)And compute with defaults:
rsaga.geoprocessor(
'ta_hydrology', 7,
list(
DEM='terrain/input_dsm/dsm.sgrd',
LENGTH='terrain/out_products/sllgt.sgrd'
),
env = env
)The Stream Power Index (SPI) is a measure of the erosive power of flowing water. SPI is calculated based upon slope and contributing area. SPI approximates locations where gullies might be more likely to form on the landscape.
To compute this index (spidx), some previous steps are required, i.e. computing the Specific Catchment Area (spcar). To do this, we use the Flow Accumulation (One Step) on the ta_hidrology module. This will also generate a flow accumulation layer (flow) that will be used later for the Vertical Distance to Channel Network.
rsaga.get.usage('terrain_analysis','Flow Accumulation (One Step)', env = env)
rsaga.geoprocessor(
'terrain_analysis',
'Flow Accumulation (One Step)',
list(
DEM='terrain/input_dsm/dsm.sgrd',
TCA='terrain/int_products/flow.sgrd',
SCA='terrain/int_products/spcar.sgrd',
PREPROCESSING=0,
FLOW_ROUTING=0
),
env = env
)
rsaga.get.usage('ta_hydrology', 21, env = env)
rsaga.geoprocessor(
'ta_hydrology', 21,
list(
SLOPE='terrain/out_products/slope.sgrd',
AREA='terrain/int_products/spcar.sgrd',
SPI='terrain/out_products/spidx.sgrd'
),
env = env
)The tool has its own module as Topographic Wetness Index (One Step). To get the usage:
rsaga.get.usage('terrain_analysis', 'Topographic Wetness Index (One Step)', env = env)To compute we run the code below. I selected Deterministic 8 [0] for the Flow Distribution. It is coded ‘twidx’.
rsaga.geoprocessor(
'terrain_analysis',
'Topographic Wetness Index (One Step)',
list(
DEM='terrain/input_dsm/dsm.sgrd',
TWI='terrain/out_products/twidx.sgrd',
FLOW_METHOD=0
),
env = env
)Again I will access the codes for the tools to get the usage:
rsaga.get.modules('ta_lighting', env = env) %>% knitr::kable()
|
From this tool we will get the hillshade ‘shade’.
rsaga.geoprocessor(
'ta_lighting', 0,
list(
ELEVATION='terrain/input_dsm/dsm.sgrd',
SHADE='terrain/out_products/shade.sgrd'
),
env = env
)Crucial variable widely used to quantify the characteristics of surface structures and estimate surface radiation budget. The SVF expresses the proportion (ratio) of radiation leaving the sky, assumed isotropic, that is able to reach a ground surface tilted at an arbitrary angle. Its value must vary between the minimum of 0, when the sky is not visible at all and to the maximum of 1, when the ground surface is horizontal and the sky entirely visible.
For this and the next tool, a multiscale parameter is chosen to speed up the computation. From this tool we can get two derivatives:
rsaga.get.usage('ta_lighting', 3, env = env)And compute with defaults:
rsaga.geoprocessor(
'ta_lighting', 3,
list(
DEM='terrain/input_dsm/dsm.sgrd',
VISIBLE='terrain/out_products/visky.sgrd',
SVF='terrain/out_products/svfct.sgrd',
METHOD=1,
DLEVEL=3
),
env = env
)Topographic openness expresses the dominance (positive) or enclosure (negative) of a landscape location. Openness has been related to how wide a landscape can be viewed from any position. It has been proven to be a meaningful input for computer aided geomorphological mapping. Openness is an angular measure of the relation between surface relief and horizontal distance. For angles less than 90”, it is equivalent to the internal angle of a cone, its apex at a DEM location, constrained by neighboring elevations within a specified radial distance. Openness incorporates the terrain line-of-sight, or viewshed, concept and is calculated from multiple zenith and nadir angles-here along eight azimuths. Openness has two viewer perspectives. Positive values, expressing openness above the surface, are high for convex forms, whereas negative values describe this attribute below the surface and are high for concave forms. Openness values are mapped by gray-scale tones. The emphasis of terrain convexity and concavity in openness maps facilitates the interpretation of landforms on the Earth’s surface and its seafloor, and on the planets, as well as features on any irregular surface-such as those generated by industrial procedures.
From this tool we can get two derivatives:
rsaga.get.usage('ta_lighting', 5, env = env)And compute with defaults:
rsaga.geoprocessor(
'ta_lighting', 5,
list(
DEM='terrain/input_dsm/dsm.sgrd',
POS='terrain/out_products/posop.sgrd',
NEG='terrain/out_products/negop.sgrd',
METHOD=0,
DLEVEL=3
),
env = env
)For this some previous steps are needed, see the complete code. The flow accumulation was already computed in the Stream Power Index step.
rsaga.get.usage('ta_channels', 0, env = env)
rsaga.geoprocessor(
'ta_channels', 0,
list(
ELEVATION='terrain/input_dsm/dsm.sgrd',
CHNLNTWRK='terrain/int_products/chnet.sgrd',
INIT_GRID='terrain/int_products/flow.sgrd',
INIT_VALUE=1000000
),
env = env
)
rsaga.get.usage('ta_channels', 3, env = env)
rsaga.geoprocessor(
'ta_channels', 3,
list(
ELEVATION='terrain/input_dsm/dsm.sgrd',
CHANNELS='terrain/int_products/chnet.sgrd',
DISTANCE='terrain/out_products/vdcnw.sgrd'
),
env = env
)Useful function to convert to .tif extension.
library(gdalUtils)
library(stringr)
dir_path_to_translate = "./terrain/out_products/"
files_short = list.files(path = dir_path_to_translate, pattern = ".sdat$", full.names = F)
files_as_tif = str_replace(files_short, ".sdat$", ".tif")
existing_files = c(
list.files(
path = "./terrain/derivatives",
pattern = ".tif$",
full.names = F,
recursive = T
) %>% strex::str_after_first('/'),
list.files(
path = "./terrain/derivatives",
pattern = ".tif$",
full.names = F
)
)
existing_files = existing_files[!is.na(existing_files)]
'%!in%' = function(x,y)!('%in%'(x,y))
files_to_translate = files_as_tif[files_as_tif %!in% existing_files] %>%
str_replace(".tif", '.sdat')
files_for_translation = file.path(paste0(dir_path_to_translate,files_to_translate))
batch_gdal_translate(
infiles = files_for_translation,
outdir = "./terrain/derivatives",
outsuffix = ".tif"
)Finally, there is a need to have every input for eCognition on the same CRS. For some reason, eCognition does not recognize this when done with GDAL. Hence, a Model Builder in ArcGIS was prepared to quickly transform everything to EPSG:2193, which is recognized and avoids conflicts. To run this, refer to the projection_tools.tbx on the ‘toolbox’ folder of the project.
Eisank, C., Hölbling, D., & Friedl, B. (2014). How well do terrain objects derived from pre-event DEM spatially correspond to landslides? Geological Society of America Abstracts with Programs, 46(6), 105. Vancouver, Canada.
Watkins, Russell, L., 2015, Terrain Metrics and Landscape Characterization from Bathymetric Data: SAGA GIS Methods and Command Sequences, Report prepared for the Ecospatial Information Team, Coral Reef Ecosystem Division, Pacific Islands Fisheries Science Center, Honolulu, HI, under NOAA contract number WE-133F-15-SE-0518, 46pp. ftp://ftp.soest.hawaii.edu/pibhmc/website/webdocs/documentation/linkages_project_methods_final.pdf