In process_gmted() or calc_gmted() (not clear for me where calculations are made since some of them are actually done inprocess_gmted()), we could add other variables derived from elevation: slope, aspect, flowdir, roughness. It is very simple using terra::terrain() function. See this old code of mine:
#' Add terrain covariates to a terra::SpatRaster (dem included)
#'
#' @param dem_path a path to dem raster
#' @param sp_vect a terra::SpatVector
#' @returns the same terra::SpatVector with terrain covariates
#' (dem, slope, aspect, roughness, flowdir)
#' @export
add_terrain <- function(sp_vect, dem_path) {
dem_rast <- terra::rast(dem_path)
names(dem_rast) <- "dem"
dem_rast$slope <- terra::terrain(dem_rast$dem, "slope")
# aspect is in degrees, clockwise from North
# if no slope: 90
dem_rast$aspect <- terra::terrain(dem_rast$dem, "aspect")
# Roughness is the difference between the maximum and the
# minimum value of a cell and its 8 surrounding cells.
dem_rast$roughness <- terra::terrain(dem_rast$dem, "roughness")
# Values encoded as power of 2 (x is the cell): IT IS A DISCRETE VARIABLE
# 32 64 128
# 16 x 1
# 8 4 2
# Cells are set to 0 if no lower neighboring.
dem_rast$flowdir <- terra::terrain(dem_rast$dem, "flowdir")
# extract at sp_vect locations
sp_vect_cov <- terra::project(sp_vect, terra::crs(dem_rast)) %>%
terra::extract(x = dem_rast, bind = TRUE) %>%
terra::project(terra::crs(sp_vect))
if (any(is.na(sp_vect_cov$dem))) {
warning("NAs found in dem column")
}
if (any(is.na(sp_vect_cov$slope))) {
warning("NAs found in slope column")
}
if (any(is.na(sp_vect_cov$aspect))) {
warning("NAs found in aspect column")
}
if (any(is.na(sp_vect_cov$roughness))) {
warning("NAs found in roughness column")
}
if (any(is.na(sp_vect_cov$flowdir))) {
warning("NAs found in flowdir column")
}
return(sp_vect_cov)
}
In
process_gmted()
orcalc_gmted()
(not clear for me where calculations are made since some of them are actually done inprocess_gmted()
), we could add other variables derived from elevation: slope, aspect, flowdir, roughness. It is very simple usingterra::terrain()
function. See this old code of mine: