renamed R script and added additional soil point data exploratory analysis plots

00_TUTORIAL_script_template.R → 00_TUTORIAL_script.R

@@ -80,7 +80,7 @@ BREAKS = unique(round(seq(MIN, MAX, RANGE/10)))
 
 # Soil point data: SOM ----------------------------------------------------
 
-# histogram of all SOM measurements
+# Histogram of all SOM measurements
 (p_hist <- tbl_regmat_target %>%
    ggplot(aes(SOM_per)) +
    geom_histogram(binwidth = 1) +
@@ -95,6 +95,46 @@ BREAKS = unique(round(seq(MIN, MAX, RANGE/10)))
          axis.title.x = element_blank(),
          legend.position = "none"))
 
+# Separate tibble into GlobalSoilMap (GSM) depth layers for boxplots
+tbl_regmat_target <- tbl_regmat_target %>%
+  filter(d_mid < 200) %>% # remove observations below 2m
+  # separate into GSM depth layers
+  mutate(d_gsm = cut(d_mid,
+                     breaks = c(0, 5, 15, 30, 60, 100, 200),
+                     labels = c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200"),
+                     right = FALSE))
+
+# Assign levels to depth layers so that they are in reverse order for boxplots
+depth_order <- c("100-200", "60-100", "30-60", "15-30", "5-15", "0-5")
+tbl_regmat_target$d_gsm <- factor(x = tbl_regmat_target$d_gsm,
+                                  levels = depth_order)
+
+# Boxplots of all SOM measurements separated into GSM depth layers
+(p_boxplot <- tbl_regmat_target %>%
+    ggplot(aes(x = SOM_per, y = d_gsm)) +
+    geom_boxplot(outlier.shape = 21) +
+    xlab(as.expression(paste(TARGET_EXP))) +
+    ylab(expression("Depth [cm]")) +
+    scale_x_continuous(breaks = BREAKS,
+                       limits = c(MIN - 0.01 * RANGE,
+                                  MAX + 0.01 * RANGE)) +
+    theme_bw() +
+    theme(strip.background = element_blank(),
+          strip.text = element_blank(),
+          legend.position = "none"))
+
+# Plot density over depth of all observations
+(p_target_density_d <- ggplot() +
+    theme_bw() +
+    geom_hex(data = tbl_regmat_target %>% 
+               filter(d_mid <= 200), # remove observations below 2m
+             aes(x = SOM_per, y = d_mid), bins = 50) +
+    scale_fill_viridis(option = "inferno", direction = -1) +
+    scale_y_continuous(trans = "reverse") +
+    xlab(TARGET_EXP) +
+    ylab("Depth [cm]") +
+    labs(fill = "Number of lab measurements"))
+
 # view all SOM measurements at point locations using mapview
 tbl_regmat_target %>% 
   # convert to simple feature (SF) object with geometry type POINT, a spatial object
@@ -103,7 +143,7 @@ tbl_regmat_target %>%
   mapview(.,
           zcol = TARGET,
           layer.name = TARGET,
-          col.regions = magma(n = 1000),
+          col.regions = hcl.colors(n = 200, palette = "YlOrBr", rev = TRUE),
           legend = TRUE,
           viewer.suppress = FALSE)
 
@@ -115,9 +155,70 @@ tbl_regmat_target_val %>%
   mapview(.,
           zcol = TARGET,
           layer.name = TARGET,
-          col.regions = magma(n = 100),
+          col.regions = hcl.colors(n = 10, palette = "YlOrBr", rev = TRUE),
+          legend = TRUE,
+          viewer.suppress = FALSE)
+
+# Since we are modelling and space AND time, let's next look at the distribution
+# of soil point data in time
+# Histogram of years during which TARGET variable was measured
+(p_hist_year <- tbl_regmat_target %>%
+    ggplot(aes(year)) +
+    geom_histogram(binwidth = 1) +
+    scale_y_continuous() +
+    scale_x_continuous(breaks = seq(1960, 2020, 10)) +
+    labs(x = "Year", y = "Count") + 
+    theme_bw() +
+    theme(strip.background = element_blank(),
+          strip.text = element_blank(),
+          axis.title.x = element_blank(),
+          legend.position = "none"))
+  
+# view all years during which TARGET variable was measured at point locations
+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 = "year",
+          layer.name = "year",
+          col.regions = viridis(n = 50),
           legend = TRUE,
           viewer.suppress = FALSE)
+# mapview legend with commas in the years is a bit annoying!
+
+# view years during which TARGET variable was measured for the first time
+# at point locations that were re-visited in 2022
+tbl_regmat_target_val %>% 
+  # 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 = "year",
+          layer.name = "year",
+          col.regions = viridis(n = 10),
+          legend = TRUE,
+          viewer.suppress = FALSE)
+
+# Density plot space (x-coordinate) vs. time, i.e. spatio-temporal distribution
+(p_target_density_t <- tbl_regmat_target %>% 
+    mutate(X = X/1000) %>% # convert coordinate units to km
+    ggplot() +
+    theme_bw() +
+    geom_hex(aes(x = year, y = X), bins = 50) +
+    scale_fill_viridis(option = "viridis") +
+    labs(x = "Year", y = "x-coordinate [km]") + 
+    labs(fill = "Number of lab measurements"))
+
+# Density plot space (y-coordinate) vs. time, i.e. spatio-temporal distribution
+(p_target_density_t <- tbl_regmat_target %>% 
+    mutate(Y = Y/1000) %>% # convert coordinate units to km
+    ggplot() +
+    theme_bw() +
+    geom_hex(aes(x = year, y = Y), bins = 50) +
+    scale_fill_viridis(option = "viridis") +
+    labs(x = "Year", y = "y-coordinate [km]") + 
+    labs(fill = "Number of lab measurements"))