---
title: "3 Prepare geological and pedological data"
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 = "#>")
```

## Prepare geological data

The geological [GÜK200](https://www.bgr.bund.de/EN/Themen/Sammlungen-Grundlagen/GG_geol_Info/Karten/Deutschland/GUEK200/guek200_node_en.html) data are accessed as sepearated shapefiles. Those are first identified:

```{r, message=FALSE,warning=FALSE}
library(magrittr)
all_shp <- list.files("../data/raw_data/vector/GUEK200", recursive = TRUE,pattern=".shp$",full.names = TRUE)

arc <- all_shp[grep("_arc.",all_shp)]
no_arc <- all_shp[-grep("_arc.",all_shp)]

```

Then the files are: 

- imported as Simple Feature objects, 
- combined, 
- additional boolean columns for moraines, mires and pleistocene fluvial sands are added and
- the Simple Feature object is stored as Geopackage.

```{r, message=FALSE,warning=FALSE}
library(sf)
library(dplyr)

guek <- Map(function(x){st_read(x) %>% 
    dplyr::select(T1_MatBez,T1_GenTxt,T1_PethTxt,T1_G_Beson,T1_Gruppe,T1_System,T1_Serie)},
    x = no_arc)
guek <- do.call(rbind,guek) %>% 
  group_by(T1_MatBez,T1_GenTxt,T1_PethTxt,T1_G_Beson,T1_Gruppe,T1_System,T1_Serie) %>% 
  summarize() %>% 
  mutate(moraine = grepl("orän",T1_GenTxt)) %>% 
  mutate(mire = grepl("moor",T1_GenTxt)) %>%
  mutate("pleistocene_fluvial_sand" = grepl("Pleisto",T1_Serie) & grepl("and",T1_PethTxt) & grepl("fluv",T1_GenTxt)) %>%
  st_cast(to ="MULTIPOLYGON") %>% 
  st_transform(25832)


st_write(guek,"../data/derived_data/vector/guek.gpkg", "guek", layer_options = c("OVERWRITE=yes"))
```

## Prepare pedological data

The pedological [BÜK200](https://www.bgr.bund.de/EN/Themen/Boden/Projekte/Informationsgrundlagen_laufend/BUEK200/BUEK200_en.html) data are accessed as sepearated shapefiles. Like the geological files, these are:

- listed
- imported as Simple Feature objects, 
- combined, 
- additional boolean columns for mires and open waters are added and
- the Simple Feature object is stored as Geopackage.

```{r, message=FALSE,warning=FALSE}
all_shp_g <- list.files("../data/raw_data/vector/BUEK200", recursive = TRUE,pattern=".shp$",full.names = TRUE)
buek1 <- Map(function(x){st_read(x) %>% 
    dplyr::select(Legende,Hinweis)},
    x = all_shp_g)
buek <- do.call(rbind,buek1) %>% 
  group_by(Legende,Hinweis) %>% 
  summarize() %>% 
  mutate(mireall = grepl("Hn|Hh",Legende),
         Legendhead = unname(unlist(Map(function(x){return(strsplit(x,":")[[1]][1])},x=as.character(Legende))))) %>%
  mutate(mirehead = grepl("HNn|HHn",Legendhead),
         waters = grepl("Gewässerflächen",Legendhead)) %>% 
  st_cast(to ="MULTIPOLYGON") %>% 
  st_transform(25832)

st_write(buek,"../data/derived_data/vector/buek.gpkg", "buek", layer_options = c("OVERWRITE=yes"))
```

## Rasterize the boolean selections

```{r}
library(raster)
library(velox)
library(fasterize)

dem_bathy <- raster("../data/derived_data/raster/dem.tif")

moraine_rast <- fasterize(sf = guek, raster = dem_bathy, field = "moraine", background = 0)
moraine_rast <- moraine_rast*!is.na(dem_bathy)
plot(moraine_rast)
moraine_rast[][moraine_rast[]==0] <- NA

# Does the Bodenübersichtskarte contains traces of mire
mire_all_rast <- fasterize(sf = buek, raster = dem_bathy, field = "mireall", background = 0)
mire_all_rast <- mire_all_rast*!is.na(dem_bathy)
plot(mire_all_rast)
writeRaster(x = mire_all_rast,
              filename = "../data/derived_data/raster/mire_all_rast.tif",
              overwrite = TRUE)


# Is the defining soil in the Bodenübersichtskarte mire
mire_head_rast <- fasterize(sf = buek, raster = dem_bathy, field = "mirehead", background = 0)
mire_head_rast <- mire_head_rast*!is.na(dem_bathy)
plot(mire_head_rast)
writeRaster(x = mire_head_rast,
              filename = "../data/derived_data/raster/mire_head_rast.tif",
              overwrite = TRUE)

# Combination of all and head rast
mire_comb_rast <- (mire_head_rast*.5) + (mire_all_rast*.5)
plot(mire_comb_rast)
writeRaster(x = mire_comb_rast,
              filename = "../data/derived_data/raster/mire_comb_rast.tif",
              overwrite = TRUE)
mire_head_rast[][mire_head_rast[]==0] <- NA


# Is a location on pleistocene_fluvial_sand
pleistocene_fluvial_sand_rast <- fasterize(sf = guek, raster = dem_bathy, field = "pleistocene_fluvial_sand", background = 0)
pleistocene_fluvial_sand_rast <- pleistocene_fluvial_sand_rast*!is.na(dem_bathy)
plot(pleistocene_fluvial_sand_rast)
writeRaster(x = pleistocene_fluvial_sand_rast,
              filename = "../data/derived_data/raster/pleistocene_fluvial_sand_rast.tif",
              overwrite = TRUE)

waters_rast <- fasterize(sf = buek, raster = dem_bathy, field = "waters", background = 0)
waters_rast <- waters_rast*!is.na(dem_bathy)
plot(waters_rast)
waters_rast[][waters_rast[]==0] <- NA

```


## Use GRASS to calculate the distances to the selected features

Since the package `rgrass7` is used to adress [GRASS GIS](https://grass.osgeo.org/) for the distance calculation, like with GDAL it is necessary to define the path of the GRASS installation on windows. The calculation itself is done using the function `r.grow.distance` as shown below.


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

writeRAST(x = as(moraine_rast, "SpatialGridDataFrame"),
          vname="moraine",
          flags = c("overwrite"))

writeRAST(x = as(mire_head_rast, "SpatialGridDataFrame"),
          vname="mire_head",
          flags = c("overwrite"))

writeRAST(x = as(waters_rast, "SpatialGridDataFrame"),
          vname="waters",
          flags = c("overwrite"))

```

### Moraine

```{r}

execGRASS("g.region",
            parameters = list(raster = "moraine",
                              align = "moraine"))

execGRASS("r.grow.distance",
            parameters = list(input = "moraine",
                              distance = "moraine_dist",                            
                              metric = "euclidean"))

moraine_distr <- raster(readRAST("moraine_dist"))

moraine_distr[is.na(dem_bathy[])] <- NA

writeRaster(x = moraine_distr,
              filename = "../data/derived_data/raster/moraine_distr.tif",
              overwrite = TRUE)
plot(moraine_distr)


```

### Mire

```{r}

execGRASS("g.region",
            parameters = list(raster = "mire_head",
                              align = "mire_head"))

execGRASS("r.grow.distance",
            parameters = list(input = "mire_head",
                              distance = "mire_head_dist",                            
                              metric = "euclidean"))

mire_head_distr <- raster(readRAST("mire_head_dist"))

mire_head_distr[is.na(dem_bathy[])] <- NA

writeRaster(x = mire_head_distr,
              filename = "../data/derived_data/raster/mire_head_distr.tif",
              overwrite = TRUE)
plot(mire_head_distr)


```

### Waters

```{r}

execGRASS("g.region",
            parameters = list(raster = "waters",
                              align = "waters"))

execGRASS("r.grow.distance",
            parameters = list(input = "waters",
                              distance = "waters_dist",                            
                              metric = "euclidean"))

waters_dist <- raster(readRAST("waters_dist"))

waters_dist[is.na(dem_bathy[])] <- NA

writeRaster(x = waters_dist,
              filename = "../data/derived_data/raster/waters_dist.tif",
              overwrite = TRUE)
plot(waters_dist)
```

