3 Prepare geological and pedological data

Prepare geological data

The geological GÜK200 data are accessed as sepearated shapefiles. Those are first identified:

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.

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)
#> Reading layer `gk1518' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek1518\gk1518.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 3904 features and 55 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3445886 ymin: 6029819 xmax: 3569318 ymax: 6119318
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
#> Reading layer `gk2310' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek2310\gk2310.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 676 features and 39 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3389680 ymin: 5940834 xmax: 3478354 ymax: 6031096
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
#> Reading layer `gk2318' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek2318\gk2318.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 3149 features and 39 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3477935 ymin: 5940783 xmax: 3566194 ymax: 6030276
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
#> Reading layer `gk2326' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek2326\gk2326.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 4455 features and 39 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3564620 ymin: 5941247 xmax: 3654442 ymax: 6054579
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
#> Reading layer `gk3118' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek3118\gk3118.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 3658 features and 55 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3477520 ymin: 5851762 xmax: 3567438 ymax: 5941247
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
#> Reading layer `gk3126' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\GUEK200\guek3126\gk3126.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 3121 features and 55 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 3566194 ymin: 5852229 xmax: 3657346 ymax: 5943314
#> epsg (SRID):    NA
#> proj4string:    +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs
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"))
#> Updating layer `guek' to data source `../data/derived_data/vector/guek.gpkg' using driver `GPKG'
#> options:        OVERWRITE=yes 
#> features:       209
#> fields:         10
#> geometry type:  Multi Polygon

Prepare pedological data

The pedological BÜK200 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.

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)
#> Reading layer `buek200_1518mg_utm_v10_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_1518\buek200_1518mg_utm_v10_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 106 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 449704.6 ymin: 6027845 xmax: 566869.2 ymax: 6101358
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
#> Reading layer `buek200_2310mg_utm_v11_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_2310\buek200_2310mg_utm_v11_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 42 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 389653.5 ymin: 5938898 xmax: 478294.1 ymax: 6027922
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
#> Reading layer `buek200_2318mg_utm_v12_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_2318\buek200_2318mg_utm_v12_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 129 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 477873.3 ymin: 5938828 xmax: 566097.3 ymax: 6028297
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
#> Reading layer `buek200_2326mg_utm_v11_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_2326\buek200_2326mg_utm_v11_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 89 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 564526.5 ymin: 5939310 xmax: 654309.5 ymax: 6050582
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
#> Reading layer `buek200_3118mg_utm_v11_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_3118\buek200_3118mg_utm_v11_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 97 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 477457.6 ymin: 5849860 xmax: 567339.3 ymax: 5939310
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
#> Reading layer `buek200_3126mg_utm_v13_poly' from data source `E:\Scales of Transformation\GitLab\b1\data\raw_data\vector\BUEK200\buek200_3126\buek200_3126mg_utm_v13_poly.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 119 features and 8 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: 566097.1 ymin: 5850327 xmax: 657210.7 ymax: 5941374
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
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"))
#> Updating layer `buek' to data source `../data/derived_data/vector/buek.gpkg' using driver `GPKG'
#> options:        OVERWRITE=yes 
#> features:       577
#> fields:         6
#> geometry type:  Multi Polygon

Rasterize the boolean selections

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 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.

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


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


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


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)