---
title: "4. Calculate Geomorphological variables"
author: "Wolfgang Hamer, Daniel Knitter"
output: html_document
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  eval = TRUE,
  results = "show",
  #echo = FALSE,
  comment = "#>")
```

# {.tabset .tabset-fade .tabset-pills}

## Function definition

```{r calc_geomporhon}
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)
}
```

## Function application

### Calculation of geommorphometric variables using function `calc_geomorphon` 

```{r}
library(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)))


writeRaster(x = geomorph,
            filename = paste0("../data/derived_data/raster/geomorph/",names(geomorph),".tif"),
            bylayer = TRUE,
            format = "GTiff",
            overwrite = TRUE)
```

### Additional computation of Geomorphons for the fuzzification of valleys

Since the package `rgrass7` is used to adress [GRASS GIS](https://grass.osgeo.org/) for the calculation of `r.geomorphon` it is necessary to define the path of the GRASS installation on windows. 

```{r}
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)
```

