2 Preparing Raster Data

3 Load data on study area

3.1 Soil samples

Code

# Load soil data from sampling locations
bern_data <- readr::read_csv(
  here::here("data-raw/soildata/berne_soil_sampling_locations.csv")
  )

# Display data
head(bern_data) |> knitr::kable()

site_id_unique	timeset	x	y	dataset	dclass	waterlog.50	waterlog.100	ph.0.10	ph.10.30	ph.30.50	ph.50.100
4_26-In-005	d1968_1974_ptf	2571994	1203001	validation	poor	0	1	6.071733	6.227780	7.109235	7.214589
4_26-In-006	d1974_1978	2572149	1202965	calibration	poor	1	1	6.900000	6.947128	7.203502	7.700000
4_26-In-012	d1974_1978	2572937	1203693	calibration	moderate	1	1	6.200000	6.147128	5.603502	5.904355
4_26-In-014	d1974_1978	2573374	1203710	validation	well	0	0	6.600000	6.754607	7.200000	7.151129
4_26-In-015	d1968_1974_ptf	2573553	1203935	validation	moderate	0	1	6.272715	6.272715	6.718392	7.269008
4_26-In-016	d1968_1974_ptf	2573310	1204328	calibration	poor	0	1	6.272715	6.160700	5.559031	5.161655

The dataset on soil samples from Bern holds 13 variables for 1052 entries (more information here):

site_id_unique: The location’s unique site id.
timeset: The sampling year and information on sampling type for soil pH (no label: CaCl\(_2\) laboratory measurement, field: indicator solution used in field, ptf: H\(_2\)O laboratory measurement transferred by pedotransfer function).
x: The x (easting) coordinates in meters following the (CH1903/LV03) system.
y: The y (northing) coordinates in meters following the (CH1903/LV03) system.
dataset: Specification whether a sample is used for model calibration or validation (this is based on randomization to ensure even spatial coverage).
dclass: Soil drainage class
waterlog.30, waterlog.50, waterlog.100: Specification whether soil was water logged at 30, 50, or 100 cm depth (0 = No, 1 = Yes).
ph.0.10, ph.10.30, ph.30.50, ph.50.100: Average soil pH between 0-10, 10-30, 30-50, and 50-100 cm depth.

3.1.1 Covariate data

Now, let’s load all the covariates that we want to produce our soil maps with.

Code

# Get a list with the path to all raster files
list_raster <- 
  base::list.files(
    here::here("data-raw/geodata/covariates/"),
    full.names = T
  )

# Take a random subset
set.seed(3)
list_raster_subset <- list_raster |> sample(15)

# Display data (lapply to clean names)
lapply(
  list_raster_subset, 
  function(x) sub(".*/(.*)", "\\1", x)
  ) |> 
  unlist() |> 
  head(5) |> 
  print()

[1] "geo500h1id.tif"           "Se_n_aspect2m.tif"       
[3] "mt_rr_y.tif"              "Se_curvplan2m_std_5c.tif"
[5] "vdcn25.tif"

The output above shows the first five raster files with rather cryptic names. The meaning of all 91 raster files are given in Chapter 6. So, make sure to have a look at the list there as it will help you to interpret your model results later on. Let’s look at one of these raster files to get a better feeling for our data. Specifically, let’s look at the slope profile at 2m resolution:

Code

# Load a raster file as example: Picking the slope profile at 2m resolution
raster_example <- terra::rast(list_raster[74])
raster_example

class       : SpatRaster 
dimensions  : 986, 2428, 1  (nrow, ncol, nlyr)
resolution  : 20, 20  (x, y)
extent      : 2568140, 2616700, 1200740, 1220460  (xmin, xmax, ymin, ymax)
coord. ref. : CH1903+ / LV95 
source      : Se_slope2m.tif 
name        : Se_slope2m 
min value   :    0.00000 
max value   :   85.11286

As shown in the output, a raster object has the following properties (among others, see ?terra::rast):

class: The class of the file, here a SpatRaster.
dimensions: The number of rows, columns, years (if temporal encoding).
resolution: The resolution of the coordinate system, here it is 20 in both axes.
extent: The extent of the coordinate system defined by min and max values on the x and y axes.
coord. ref.: Reference coordinate system. Here, the raster is encoded using the LV95 geodetic reference system from which the projected coordinate system CH1903+ is derived.
source: The name of the source file.
names: The name of the raster file (mostly the file name without file-specific ending)
min value: The lowest value of all cells.
max value: The highest value of all cells.

Tip

The code chunks filtered for a random sub-sample of 15 variables. As described in Chapter 5, your task will be to investigate all covariates and find the ones that can best be used for your modelling task.

4 Mapping the study area

Now, let’s look at a visualisation of this raster file. Since we selected the slope at 2m resolution, we expect a relief-like map with a color gradient that indicates the steepness of the terrain.

Code

# Plot raster example
terra::plot(raster_example)

Code

# To have some more flexibility, we can plot this in the ggplot-style as such:
ggplot2::ggplot() +
  tidyterra::geom_spatraster(data = raster_example) +
  ggplot2::scale_fill_viridis_c(
    na.value = NA,
    option = "magma",
    name = "Slope (%) \n"
    ) +
  ggplot2::theme_bw() +
  ggplot2::scale_x_continuous(expand = c(0, 0)) +  # avoid gap between plotting area and axis
  ggplot2::scale_y_continuous(expand = c(0, 0)) +
  ggplot2::labs(title = "Slope of the Study Area")

Tip

Note that the second plot has different coordinates than the upper one. That is because the data was automatically projected to the World Geodetic System (WGS84, ESPG: 4326).

This looks already interesting but we can put our data into a bit more context. For example, a larger map background would be useful to get a better orientation of our location. Also, it would be nice to see where our sampling locations are and to differentiate these locations by whether they are part of the calibration or validation dataset. Bringing this all together requires some more understanding of plotting maps in R. So, don’t worry if you do not understand everything in the code chunk below and enjoy the visualizations:

Code

# To get our map working correctly, we have to ensure that all the input data
# is in the same coordinate system. Since our Bern data is in the Swiss 
# coordinate system, we have to transform the sampling locations to the 
# World Geodetic System first.
# To look up EPSG Codes: https://epsg.io/
# World Geodetic System 1984:  4326
# Swiss CH1903+ / LV95: 2056

# For the raster:
r <- terra::project(raster_example, "+init=EPSG:4326")

# Let's make a function for transforming the sampling locations:
change_coords <- function(data, from_CRS, to_CRS) {
  
  # Check if data input is correct
  if (!all(names(data) %in% c("id", "lat", "lon"))) {
    stop("Input data needs variables: id, lat, lon")
  }
  
  # Create simple feature for old CRS
  sf_old_crs <- sf::st_as_sf(data, coords = c("lon", "lat"), crs = from_CRS)
  
  # Transform to new CRS
  sf_new_crs     <- sf::st_transform(sf_old_crs, crs = to_CRS)
  sf_new_crs$lat <- sf::st_coordinates(sf_new_crs)[, "Y"]
  sf_new_crs$lon <- sf::st_coordinates(sf_new_crs)[, "X"]
  
  sf_new_crs <- sf_new_crs |> dplyr::as_tibble() |> dplyr::select(id, lat, lon)
  
  # Return new CRS
  return(sf_new_crs)
}

# Transform dataframes
coord_cal <- 
  bern_data |> 
  dplyr::filter(dataset == "calibration") |> 
  dplyr::select(site_id_unique, x, y) |> 
  dplyr::rename(id = site_id_unique, lon = x, lat = y) |> 
  change_coords(
    from_CRS = 2056, 
    to_CRS = 4326
    )

coord_val <- 
  bern_data |> 
  dplyr::filter(dataset == "validation") |> 
  dplyr::select(site_id_unique, x, y) |> 
  dplyr::rename(id = site_id_unique, lon = x, lat = y) |> 
  change_coords(
    from_CRS = 2056, 
    to_CRS = 4326
    )

Code

# Notes: 
# - This code may only work when installing the development branch of {leaflet}:
# remotes::install_github('rstudio/leaflet')
# - You might have to do library(terra) for R to find functions needed in the backend
library(terra)

# Let's get a nice color palette now for easy reference
pal <- leaflet::colorNumeric(
  "magma",
  terra::values(r),
  na.color = "transparent"
  )

# Next, we build a leaflet map
leaflet::leaflet() |> 
  # As base maps, use two provided by ESRI
  leaflet::addProviderTiles(leaflet::providers$Esri.WorldImagery, group = "World Imagery") |>
  leaflet::addProviderTiles(leaflet::providers$Esri.WorldTopoMap, group = "World Topo") |>
  # Add our raster file
  leaflet::addRasterImage(
    r,
    colors = pal,
    opacity = 0.6,
    group = "raster"
    ) |>
  # Add markers for sampling locations
  leaflet::addCircleMarkers(
    data = coord_cal,
    lng = ~lon,  # Column name for x coordinates
    lat = ~lat,  # Column name for y coordinates
    group = "training",
    color = "black"
  ) |>
    leaflet::addCircleMarkers(
    data = coord_val,
    lng = ~lon,  # Column name for x coordinates
    lat = ~lat,  # Column name for y coordinates
    group = "validation",
    color = "red"
  ) |>
  # Add some layout and legend
  leaflet::addLayersControl(
    baseGroups = c("World Imagery","World Topo"),
    position = "topleft",
    options = leaflet::layersControlOptions(collapsed = FALSE),
    overlayGroups = c("raster", "training", "validation")
    ) |>
  leaflet::addLegend(
    pal = pal,
    values = terra::values(r),
    title = "Slope (%)")