Skip to content
Snippets Groups Projects
Commit b1043b60 authored by Helfenstein, Anatol's avatar Helfenstein, Anatol :speech_balloon:
Browse files

derive dynamic land use covariates

parent 85c1f651
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 130 additions and 24 deletions
00_TUTORIAL_script.R
+ 130
− 24
View file @ b1043b60
......@@ -20,8 +20,8 @@
# Load required packages --------------------------------------------------
pkgs <- c("tidyverse", "foreach", "raster", "sf", "mapview", "ranger", "caret",
"viridis")
pkgs <- c("tidyverse", "foreach", "raster", "terra", "sf", "mapview", "ranger",
"caret", "viridis")
# make sure 'mapview' pkg installed from github to avoid pandoc error:
# remotes::install_github("r-spatial/mapview")
lapply(pkgs, library, character.only = TRUE)
......@@ -70,6 +70,18 @@ ls_r_cov <- foreach(cov = 1:length(v_dir_cov)) %do%
r_stack_cov <- stack(ls_r_cov)
# locate, read in and stack covariates from which we will derive dynamic (2D+T
# and 3D+T) covariates
v_cov_dyn_names <- dir("out/data/covariates/final_stack_dyn",
pattern = "\\.grd$", recursive = TRUE) %>%
.[c(1:7, 12:16, 8:11)] # change order of HGNs and LGNs by year of creation
ls_r_cov_dyn <- foreach(cov = 1:length(v_cov_dyn_names)) %do%
raster(paste0("out/data/covariates/final_stack_dyn/",
v_cov_dyn_names[[cov]]))
r_stack_cov_dyn <- stack(ls_r_cov_dyn)
# 9) Set TARGET variable min, max, range and breaks for plotting purposes
MIN = min(tbl_regmat_target_cal[TARGET])
MAX = max(tbl_regmat_target_cal[TARGET])
......@@ -78,7 +90,7 @@ BREAKS = unique(round(seq(MIN, MAX, RANGE/10)))
# Soil point data: SOM ----------------------------------------------------
# Exploratory analysis: soil point data -----------------------------------
# Histogram of all SOM measurements
(p_hist <- tbl_regmat_target %>%
......@@ -206,7 +218,7 @@ tbl_regmat_target_val %>%
ggplot() +
theme_bw() +
geom_hex(aes(x = year, y = X), bins = 50) +
scale_fill_viridis(option = "viridis") +
scale_fill_viridis(option = "inferno", direction = -1) +
labs(x = "Year", y = "x-coordinate [km]") +
labs(fill = "Number of lab measurements"))
......@@ -216,15 +228,25 @@ tbl_regmat_target_val %>%
ggplot() +
theme_bw() +
geom_hex(aes(x = year, y = Y), bins = 50) +
scale_fill_viridis(option = "viridis") +
scale_fill_viridis(option = "inferno", direction = -1) +
labs(x = "Year", y = "y-coordinate [km]") +
labs(fill = "Number of lab measurements"))
# Static covariates -------------------------------------------------------
# Exploratory analysis: covariates ----------------------------------------
# Names of static covariates included for this practical (NOTE: this does not
# include all covariates used in BIS-4D to save some time...)
names(r_stack_cov)
# View some covariates
# Names of categorical covariates (the others are continuous)
names(r_stack_cov)[is.factor(r_stack_cov)]
# View classes of one categorical covariate (land use 2021)
levels(r_stack_cov$lgn2021_1km)
# Map some static covariates
# Climate: Long-term mean precipitation
mapView(r_stack_cov$precip_yearly_1981_2010_1km,
col.regions = mako(n = 100, direction = -1),
......@@ -235,7 +257,7 @@ mapView(r_stack_cov$precip_yearly_1981_2010_1km,
mapView(r_stack_cov$s2_GSI_2018_11nov_1km,
col.regions = rocket(n = 100),
layer.name = "Grain Size Index")
# some mosaic artifacts are clearly visible
# Some mosaic artifacts are clearly visible
# Relief: AHN (national digital elevation model of the Netherlands)
mapView(r_stack_cov$ahn3_1km,
......@@ -267,6 +289,7 @@ mapView(r_stack_cov$lgn2021_1km,
"#386cb1", "#e6e1bc", "#e6e1bd", "#e6e1be", "#e6e1bf",
"#e6e1bb", "#e6e1ba"),
layer.name = "Land use in 2021")
# what differences do you see between the land use map from 1960 and 2021?
# Peat classes based on peat thickness and starting depth from original national
# soil map of the Netherlands (Bodemkaart 1:50'000)
......@@ -282,20 +305,103 @@ mapView(r_stack_cov$bodem50_2021_peatcode_1km,
col.regions = c("#77128e", "#da6dd1", "#3c46a5", "#44b3f4", "#1b829e",
"#98d4f6", "#256c57", "#b9f8a1", "white"),
layer.name = "Peat classes (2013-2021)")
# view all point locations colored by measurement year using mapview
tbl_regmat_target %>%
# convert to simple feature (SF) object with geometry type POINT, a spatial object
st_as_sf(., coords = c("X", "Y"), crs = "EPSG:28992") %>%
st_jitter(., factor = 0.00001) %>% # so we can see samples from same location
mapview(.,
zcol = TARGET,
layer.name = TARGET,
col.regions = magma(n = 1000),
legend = TRUE,
viewer.suppress = FALSE)
# what differences do you see between the original and updated maps of peat classes?
# Dynamic covariates: land use (LU) ---------------------------------------
# Date of data used to make each national LU map (LGN),
# see https://www.wur.nl/nl/Onderzoek-Resultaten/Onderzoeksinstituten/Environmental-Research/Faciliteiten-tools/Kaarten-en-GIS-bestanden/Landelijk-Grondgebruik-Nederland/Versies-bestanden.htm
# LGN1: 1984-1987
# LGN2: 1990-1994
# LGN3: 1995-1997
# LGN4: 1999-2000
# LGN5: 2003-2004
# LGN6: 2007-2008
# LGN7: 2012
# LGN2018: 2018
# LGN2019: 2019
# LGN2020: 2020
# LGN2021: 2021
# Historic LU maps (HGN) reconstructed for specific year (+/- 5 years) using old maps
# HGN1900: 1900
# HGN1960: 1960
# HGN1970: 1970
# HGN1980: 1980
# HGN1990: 1990
# To create dynamic LU covariates, first create temporary cols of
# a) LU for each year at sampling location (all the way back to LU 40 yrs
# before first year we predict for (see delta40 below), i.e. 1950 - 40 = 1910
# b) using a, we can then make cols for LU at each sampling location t-1 all
# the way up to t-40 yrs ago
# step a (see above)
r_stack_LU_1910_present <- stack(
stack(replicate(length(1910:1930), r_stack_cov_dyn$hgn_1900_filled_1km)),
stack(replicate(length(1931:1965), r_stack_cov_dyn$hgn_1960_1km)),
stack(replicate(length(1966:1975), r_stack_cov_dyn$hgn_1970_1km)),
stack(replicate(length(1976:1982), r_stack_cov_dyn$hgn_1980_1km)),
stack(replicate(length(1983:1988), r_stack_cov_dyn$lgn1_1km)),
stack(replicate(length(1989:1990), r_stack_cov_dyn$hgn_1990_1km)),
stack(replicate(length(1991:1994), r_stack_cov_dyn$lgn2_1km)),
stack(replicate(length(1995:1998), r_stack_cov_dyn$lgn3_1km)),
stack(replicate(length(1999:2001), r_stack_cov_dyn$lgn4_1km)),
stack(replicate(length(2002:2005), r_stack_cov_dyn$lgn5_1km)),
stack(replicate(length(2006:2010), r_stack_cov_dyn$lgn6_1km)),
stack(replicate(length(2011:2015), r_stack_cov_dyn$lgn7_1km)),
stack(replicate(length(2016:2018), r_stack_cov_dyn$lgn2018_1km)),
stack(replicate(1, r_stack_cov_dyn$lgn2019_1km)),
stack(replicate(1, r_stack_cov_dyn$lgn2020_1km)),
stack(replicate(length(2021:format(Sys.Date(), "%Y")), r_stack_cov_dyn$lgn2021_1km))
)
# rename newly created rasters
names(r_stack_LU_1910_present) <- paste0("LU_", seq(1910, format(Sys.Date(), "%Y"), 1))
# calculating mode using raster pkg takes 100 times longer than using terra,
# so convert to spatRaster using terra and then perform terra::modal()
sr_LU_1910_present <- rast(r_stack_LU_1910_present)
# Obtain land use for every location, with coordinates x and y for every year t
# between 1953 and 2022, by assigning the same class as in the temporally nearest
# year for which a map is available.
# LU_xyt_1km for year 1953 & 2022
sr_LU_xyt_1km <- imap(c(1953, 2022), ~ {
i <- .x
sr_LU_1910_present[[paste0("LU_", i)]]
}) %>% rast(.)
# Sidenote: since output of foreach is list, rast() combines into multilayer SpatRaster
# In the same manner, assign the land use class that occurred most frequently in
# the 5, 10, 20 and 40 years prior to and including year t...
# LU_xyt_delta5_1km for year 1953 & 2022
sr_LU_xyt_delta5_1km <- imap(c(1953, 2022), ~ {
i <- .x
modal(sr_LU_1910_present[[paste0("LU_", seq(i - 4, i))]], na.rm = TRUE)
}) %>% rast(.)
# LU_xyt_delta10_1km for year 1953 & 2022
sr_LU_xyt_delta10_1km <- imap(c(1953, 2022), ~ {
i <- .x
modal(sr_LU_1910_present[[paste0("LU_", seq(i - 9, i))]], na.rm = TRUE)
}) %>% rast(.)
# LU_xyt_delta20_1km for year 1953 & 2022
sr_LU_xyt_delta20_1km <- imap(c(1953, 2022), ~ {
i <- .x
modal(sr_LU_1910_present[[paste0("LU_", seq(i - 19, i))]], na.rm = TRUE)
}) %>% rast(.)
# LU_xyt_delta40_1km for year 1953 & 2022
sr_LU_xyt_delta40_1km <- imap(c(1953, 2022), ~ {
i <- .x
modal(sr_LU_1910_present[[paste0("LU_", seq(i - 39, i))]], na.rm = TRUE)
}) %>% rast(.)
# Dynamic covariates: peat occurrence -------------------------------------
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment