calc_geomorphon <- function(
elevation,
gdalpath = "C:/OSGeo4W64/bin",
grasspath = "C:/Program Files/GRASS GIS 7.4.2",
writepath = FALSE,
variables = c("Smooth_DEM",
"Flow_accumulation",
"Terrain_wetness",
"Terrain_position",
"pcurvature",
"tcurvature",
"Slope",
"Positive_openness",
"Negative_openness"),
variable_parameters = FALSE){
if(length(strsplit(as.character(crs(elevation)),
"epsg:4326")[[1]]) > 1){
warning("Projected coordinate systems with meters as a unit are recommended when using this function!")
}
## Define for windows where GDAL is located
if(.Platform$OS.type == "windows"){
gdal.dir <- shortPathName(gdalpath)
gdal_translate <- paste0(gdal.dir, "/gdal_translate.exe")
gdalwarp <- paste0(gdal.dir, "/gdalwarp.exe")
gdalinfo <- paste0(gdal.dir, "/gdalinfo.exe")
}else{
gdal_translate = "gdal_translate"
gdalwarp = "gdalwarp"
gdalinfo = "gdalinfo"
}
if(.Platform$OS.type == "windows"){
grasspath <- grasspath
}else{
grasspath <- "/opt/grass"
}
loc <- rgrass7::initGRASS(grasspath,
home = tempdir(),
mapset = "PERMANENT",
override = TRUE)
execGRASS("g.proj",
flags = c("c"),
parameters = list(proj4 = as.character(crs(elevation))))
writeRAST(x = as(elevation, "SpatialGridDataFrame"),
vname="dem",
flags = c("overwrite"))
execGRASS("g.region",
parameters = list(raster = "dem",
align = "dem"))
# Smoothing elevation model
if(any(is.element("Smooth_DEM",variables))){
smlist <- list(size = 11,
input = "dem",
output = "dem_smooth")
if(!is.logical(variable_parameters)){
if(is.element("Smooth_DEM",names(variable_parameters))){
smlist <- modifyList(smlist,variable_parameters$Smooth_DEM)
smlist <- modifyList(smlist,list(input = "dem",
output = "dem_smooth"))
}
}
execGRASS("r.param.scale",
parameters = smlist,
flags = c("overwrite")
)
Smooth_DEM <- raster(readRAST("dem_smooth"))
} else {
writeRAST(x = as(elevation, "SpatialGridDataFrame"),
vname="dem_smooth",
flags = c("overwrite")
)
Smooth_DEM <- elevation
}
if(any(is.element(c("pcurvature","tcurvature"),variables))){
curvlist <- list(elevation = "dem_smooth",
pcurvature = "pcurvature",
tcurvature = "tcurvature")
if(!is.logical(variable_parameters)){
if(is.element("pcurvature",names(variable_parameters))){
curvlist <- modifyList(curvlist,variable_parameters$pcurvature)
curvlist <- modifyList(curvlist,list(elevation = "dem_smooth",
pcurvature = "pcurvature",
tcurvature = "tcurvature"))
}
if(is.element("tcurvature",names(variable_parameters))){
curvlist <- modifyList(curvlist,variable_parameters$tcurvature)
curvlist <- modifyList(curvlist,list(elevation = "dem_smooth",
pcurvature = "pcurvature",
tcurvature = "tcurvature"))
}
}
execGRASS("r.slope.aspect",
flags = c("overwrite"),
parameters = curvlist
)
pcurvature <- raster(readRAST("pcurvature"))
tcurvature <- raster(readRAST("tcurvature"))
}
if(any(is.element("Slope",variables))){
slopelist <- list(elevation = "dem_smooth",
slope = "slope")
if(!is.logical(variable_parameters)){
if(is.element("Slope",names(variable_parameters))){
slopelist <- modifyList(slopelist,variable_parameters$Slope)
slopelist <- modifyList(slopelist,list(elevation = "dem_smooth",
slope = "slope"))
}
}
execGRASS("r.slope.aspect",
flags = c("overwrite"),
parameters = slopelist
)
Slope <- readRAST("slope") %>%
raster::raster()
}
if(any(is.element("Terrain_position",variables))){
tpilist <- list(w = 5)
if(!is.logical(variable_parameters)){
if(is.element("Terrain_position",names(variable_parameters))){
tpilist <- modifyList(tpilist,variable_parameters$Terrain_position)
}
}
tpiw <- function(x, w = 5) {
m <- matrix(1 / (w^2 - 1), nc = w, nr = w)
m[ceiling(0.5 * length(m))] <- 0
f <- focal(x, m)
x - f
}
Terrain_position <- tpiw(Smooth_DEM, w = tpilist$w)
names(Terrain_position) <- "tpi"
}
#### Hydrology related parameter
if(any(is.element(c("Flow_accumulation","Terrain_wetness"),variables))){
wat <- list(elevation = "dem_smooth",
accumulation = "fac",
tci = "twi")
if(!is.logical(variable_parameters)){
if(is.element("Flow_accumulation",names(variable_parameters))){
wat <- modifyList(wat,variable_parameters$Flow_accumulation)
wat <- modifyList(wat,list(elevation = "dem_smooth",
accumulation = "fac",
tci = "twi"))
}
if(is.element("Terrain_wetness",names(variable_parameters))){
wat <- modifyList(wat,variable_parameters$Terrain_wetness)
wat <- modifyList(wat,list(elevation = "dem_smooth",
accumulation = "fac",
tci = "twi"))
}
}
execGRASS("r.watershed",
flags = c("overwrite", "b"),
parameters = wat
)
Flow_accumulation <- raster(readRAST(vname = "fac"))
Terrain_wetness <- raster(readRAST(vname = "twi"))
}
#### Topographic Openess
if(any(is.element(c("Positive_openness","Negative_openness"),variables))){
tmp <- tempdir()
writeRaster(Smooth_DEM,
paste0(tmp,"/dem_saga.sgrd"),
overwrite = TRUE)
env <- rsaga.env()
topolist <- list(DEM = paste0(tmp,"/dem_saga.sgrd"),
POS = paste0(tmp,"/pos.sgrd"),
NEG = paste0(tmp,"/neg.sgrd"))
if(!is.logical(variable_parameters)){
if(is.element("Positive_openness",names(variable_parameters))){
topolist <- modifyList(topolist,variable_parameters$Positive_openness)
}
if(is.element("Negative_openness",names(variable_parameters))){
topolist <- modifyList(topolist,variable_parameters$Negative_openness)
}
topolist <- modifyList(topolist,list(DEM = paste0(tmp,"/dem_saga.sgrd"),
POS = paste0(tmp,"/pos.sgrd"),
NEG = paste0(tmp,"/neg.sgrd")))
}
rsaga.geoprocessor(lib = "ta_lighting", module = "Topographic Openness",
param = topolist,
env = env)
Positive_openness <- raster(paste0(tmp,"/pos.sdat"))
Negative_openness <- raster(paste0(tmp,"/neg.sdat"))
unlink(tmp)
Positive_openness <- mask(resample(Positive_openness,
Smooth_DEM,
method="ngb"),
Smooth_DEM)
Negative_openness <- mask(resample(Negative_openness,
Smooth_DEM,
method="ngb"),
Smooth_DEM)
}
parameters <- lapply(X = as.list(variables),
FUN = function(x) {
eval(as.name(x))}) %>%
raster::stack()
if(!is.logical(writepath)){
dir.create(writepath, showWarnings = FALSE)
writeRaster(parameters,
paste0(writepath,"/",names(parameters),".tif"),
bylayer = T,
overwrite=TRUE)
}
return(parameters)
}calc_geomorphonlibrary(magrittr)
library(raster)
library(rgrass7)
library(RSAGA)
dem_bathy <- raster("../data/derived_data/raster/dem.tif")
geomorph <- calc_geomorphon(dem_bathy,
variable_parameters =
list(Negative_openness=
list(RADIUS = 1000)))
#> Search for SAGA command line program and modules...
#> Done
#>
#>
#> SAGA Version: 2.3.2 (64 bit)
#>
#> library path: C:\OSGEO4~1\apps\saga-ltr\modules\
#> library name: ta_lighting
#> library : Lighting, Visibility
#> tool : Topographic Openness
#> author : O.Conrad (c) 2012
#> processors : 8 [8]
#>
#> Load grid: C:\Users\hamer\AppData\Local\Temp\RtmpQVh1Bp/dem_saga.sgrd...
#>
#>
#> Parameters
#>
#>
#> Grid system: 25.03983; 6430x 6828y; 489135.073215x 5913474.474488y
#> Elevation: dem_saga
#> Positive Openness: Positive Openness
#> Negative Openness: Negative Openness
#> Radial Limit: 1000.000000
#> Method: sectors
#> Multi Scale Factor: 3.000000
#> Number of Sectors: 8
#>
#> Save grid: C:\Users\hamer\AppData\Local\Temp\RtmpQVh1Bp/pos.sgrd...
#> Save grid: C:\Users\hamer\AppData\Local\Temp\RtmpQVh1Bp/neg.sgrd...
writeRaster(x = geomorph,
filename = paste0("../data/derived_data/raster/geomorph/",names(geomorph),".tif"),
bylayer = TRUE,
format = "GTiff",
overwrite = TRUE)Since the package rgrass7 is used to adress GRASS GIS for the calculation of r.geomorphon it is necessary to define the path of the GRASS installation on windows.
if(.Platform$OS.type == "windows"){
grasspath <- "C:/Program Files/GRASS GIS 7.4.2"
}else{
grasspath <- "/opt/grass"
}
loc <- rgrass7::initGRASS(grasspath,
home = tempdir(),
mapset = "PERMANENT",
override = TRUE)
execGRASS("g.proj",
flags = c("c"),
parameters = list(proj4 = as.character(crs(dem_bathy))))
writeRAST(x = as(dem_bathy, "SpatialGridDataFrame"),
vname="dem",
flags = c("overwrite"))
execGRASS("g.region",
parameters = list(raster = "dem",
align = "dem"))
geomlist <- list(elevation = "dem",
forms = "geomorphons",
search = 5000)
execGRASS("r.geomorphon",
flags = c("overwrite", "verbose", "m"),
parameters = geomlist
)
geomorphons <- raster(readRAST(vname = "geomorphons"))
writeRaster(x = geomorphons,
filename = paste0("../data/derived_data/raster/geomorph/geomorphons.tif"),
bylayer = TRUE,
format = "GTiff",
overwrite = TRUE)