Merge pull request #6 from ac-willeke/main

Updates requested by Magni K. in march 2023

.gitignore

@@ -1,7 +1,12 @@
+# R
 .Rproj.user
 .Rhistory
 .RData
 .Ruserdata
 data_exple/
 docker_run_rstudio.sh
 rstudio-prefs.json
+
+# ignore base session data files
+data/*
+!data/gran_dataset.csv

Dockerfile

@@ -2,7 +2,8 @@
 FROM rocker/geospatial
 
 RUN apt-get update -qq && \
-    apt-get install -qq libxt-dev r-cran-cairo
+    update-ca-certificates && \
+    apt-get install -qq libxt-dev r-base-dev r-cran-cairo 
 
 RUN install2.r --error \
     --deps TRUE \
@@ -26,7 +27,4 @@ COPY . ./home/rstudio/app/
 RUN mkdir -p ./home/rstudio/.config/rstudio/
 RUN cp ./home/rstudio/app/rstudio-prefs.json ./home/rstudio/.config/rstudio/
 
-CMD Rscript ./home/rstudio/app/app.R
-
-
-
+CMD Rscript ./home/rstudio/app/app.R

Rscripts/exploration.R

@@ -0,0 +1,85 @@
+library(tidyverse)
+library(data.table)
+library(gstat)
+
+# Define file paths
+shapefile_path <- "/home/rstudio/app/data_exple/test_dataset/test_shapefile.shp"
+csv_path <- "/home/rstudio/app/data_exple/test_dataset/test_depths_samples.csv"
+
+# Read shapefile and convert to simple features object
+shp <- rgdal::readOGR(shapefile_path) %>% 
+  sf::st_as_sf()
+
+input_crs = st_crs(shp)
+target_crs = 25833
+
+# Read CSV data
+df <- open_csv(csv_path)
+
+# Transform CSV data to spatial format and set its CRS to match the shapefile
+dfs <- transform_to_sf(df) %>% 
+  sf::st_set_crs(sf::st_crs(input_crs))
+
+# Transform both datasets to the target CRS
+dfs <- dfs %>% sf::st_transform(target_crs)
+shp <- shp %>% sf::st_transform(target_crs)
+
+myGrid <- starsExtra::make_grid(shp, 1)
+myGrid <- sf::st_crop(myGrid, shp)
+
+powerRange <- 1:6
+nmax <- 20
+temp <- data.frame(power = powerRange, MAE = rep(NA, length(powerRange)))
+vol <- numeric()
+
+# Calc. the Volume for kriging and IDW with different power values
+# Cross-validation used to calc. MAE and optimal power value
+# Optimal power value is used to calc the volume
+
+
+# Final interpolation and volume calc using optimal power value
+for(i in 2){
+
+  # Get the MAE
+  temp2 <- gstat::krige.cv(torvdybde_cm ~ 1, dfs, set = list(idp=i), nmax = nmax)
+  temp$MAE[temp$power==i] <- mean(abs(temp2$residual))
+
+  #  Get the volume in cm3 
+  vol_temp <- gstat::idw(torvdybde_cm ~ 1, dfs, 
+                         newdata=myGrid, 
+                         nmax=nmax, 
+                         idp=i)
+  # vol in m3
+  vol <- c(vol, sum(vol_temp$var1.pred, na.rm=T))
+  vol_m3 <- vol/100
+
+}
+
+idwe <- gstat::idw(formula = torvdybde_cm ~ 1, 
+                        locations = dfs, 
+                        newdata = myGrid, 
+                        idp=4,
+                        nmax = nmax)
+
+
+idwe_r <- as(idwe, "Raster")
+
+library(leafem)
+
+pal = colorNumeric(palette = "magma", 
+                   values(idwe_r),  na.color = "transparent",
+                   reverse = TRUE)
+
+leaflet() %>% addTiles() %>%
+  addStarsImage(idwe, colors = pal) %>% 
+  addLegend(pal = pal, values = values(idwe_r),
+            title = "Torvdybde (cm)") %>% 
+  #addControl("title", position = "topleft", className="map-title")
+
+leaflet() %>% addTiles() %>%
+  addRasterImage(idwe_r, colors = pal) #%>% 
+  #addLegend(pal = pal, values = values(idwe_r),
+  #          title = "torvdybde_cm") %>% 
+  #addControl(title, position = "topleft", className="map-title")
+
+

app.R

@@ -8,6 +8,8 @@ source("/home/rstudio/app/appScripts/server.R")
 options(shiny.maxRequestSize=100*1024^2)
 options(warn = -1)
 
+options(shiny.autoreload=TRUE)
+
 app <- shinyApp(ui = ui, server = server)
 sessionInfo()
 runApp(app, host ="0.0.0.0", port = 8999, launch.browser = TRUE)

appScripts/global.R

@@ -9,14 +9,33 @@ compute_volume_slider <- function(peatDepths,
   myGrid <- starsExtra::make_grid(peatlandDelimination, 1)
   myGrid <- sf::st_crop(myGrid, peatlandDelimination)
 
-  vol <- gstat::idw(dybde ~ 1, peatDepths, 
+  vol <- gstat::idw(torvdybde_cm ~ 1, peatDepths, 
                              newdata=myGrid, 
                              nmax=30, 
                              idp=power)
 
   return(sum(vol$var1.pred, na.rm=T))
 }
 
+#' Interpolation of peat depths
+#'
+#' Interpolation of peat depths using Inverse Distance Weighting (IDW). 
+#' The optimal power value is determined by minimizing the Mean Absolute Error (MAE) 
+#' using cross-validation (IDW vs. kriging).
+#' 
+
+#' @param peatDepths peat depth points.
+#' @param peatlandDelimination study area polygon.
+#' @param powerRange range of power values: Default = 1:6.
+#' @param nmax max number of points to use in IDW and kriging. Default = 30.
+#'
+#' @return A list containing the following elements:
+#' - v: volume of peat in the study area.
+#' - idweights: A Spatial object containing the interpolated values.
+#' - best: The best power value.
+#' - gg_out: ggplot showing MAE against power.
+#' - gg_out_vol: A ggplot showing relative volume against power.
+
 interpolation <- function(peatDepths,
                           peatlandDelimination,
                           powerRange = 1:6,
@@ -36,15 +55,16 @@ interpolation <- function(peatDepths,
     for(i in powerRange){
 
       # Get the MAE
-      temp2 <- gstat::krige.cv(dybde ~ 1, peatDepths, set = list(idp=i), nmax = nmax)
+      temp2 <- gstat::krige.cv(torvdybde_cm ~ 1, peatDepths, set = list(idp=i), nmax = nmax)
       temp$MAE[temp$power==i] <- mean(abs(temp2$residual))
 
       #  Get the volume
-      vol_temp <- gstat::idw(dybde ~ 1, peatDepths, 
+      vol_temp <- gstat::idw(torvdybde_cm ~ 1, peatDepths, 
                              newdata=myGrid, 
                              nmax=nmax, 
                              idp=i)
 
+      # vol in m3
       vol <- c(vol, sum(vol_temp$var1.pred, na.rm=T))
 
       Sys.sleep(0.5)
@@ -60,7 +80,7 @@ interpolation <- function(peatDepths,
   best <- temp %>% filter(best == "best")
   best <- best$power
 
-  
+
   ###############
   # Plot volume #
   ###############
@@ -69,8 +89,9 @@ interpolation <- function(peatDepths,
   vol_df$relative_volume <- vol_df$volume/mean(vol_df$volume)*100
   v <- vol_df %>% filter(power == best)
   v <- v$volume
-
-  idweights <- gstat::idw(formula = dybde ~ 1, 
+  v_m3 <- v/100  
+
+  idweights <- gstat::idw(formula = torvdybde_cm ~ 1, 
              locations = peatDepths, 
              newdata = myGrid, 
              idp=best,
@@ -81,25 +102,30 @@ interpolation <- function(peatDepths,
   ######################
   # Plot MAE vs power
   gg_out <- ggplot(temp, aes(x = power, y = MAE,
-                             colour = best,
-                             shape = best))+
-    geom_point(size=10)+
-    theme_bw(base_size = 12)+
+                            colour = best,
+                            shape = best))+
+    geom_point(size=6)+
+    theme_bw(base_size = 14)+
     scale_x_continuous(breaks = powerRange)+
     guides(colour="none",
-           shape = "none")+
+          shape = "none")+
     scale_color_manual(values = c("darkgreen","grey"))+
-    scale_shape_manual(values = c(18, 19))
-
+    scale_shape_manual(values = c(18, 19))+
+    labs(title = "Comparison of Mean Absolute Error (MAE) for Different Power Values in IDW",
+        subtitle = "The optimal power value is visualized by the green diamond.",
+        x = "Power",
+        y = "Mean Absolute Error (MAE)")
+
   # Plot relative volume vs power
   gg_out_vol <- ggplot(vol_df, aes(x = factor(power), y = relative_volume))+
-    geom_point(size=8)+
-    xlab("power")+
-    ylab("Peat volume as a percentage of\nmean predicted peat volume")+
-    theme_bw(base_size = 12)
+    geom_point(size=6)+
+    theme_bw(base_size = 14)+
+    labs(title = "Impact of Power Parameter on Estimated Peat Volume",
+       x = "Power",
+       y = "Peat volume as a percentage of\nmean predicted peat volume")
 
   # Return the interpolation and the result volume
-  l_results <- list(v, idweights, best, gg_out, gg_out_vol)
+  l_results <- list(v_m3, idweights, best, gg_out, gg_out_vol)
   return(l_results)
 }
 
@@ -159,7 +185,7 @@ print_error_incompatible_file <- function(){
     "Please upload a dataset before clicking on 'load dataset'. 
           The dataset should be a zip file containing both a shapefile with the 
           extent of the area of interest and a csv file containing peat depth measures 
-          (in m) taken at the site with coordinates for each measure (given in UTM 32 N, EPSG:25832). 
+          (in m) taken at the site with coordinates for each measure (given in UTM 33 N, EPSG:25833). 
           For more information refer to the README.",
     easyClose = TRUE
   ))
@@ -169,11 +195,19 @@ print_error_csv_columns <- function(){
   showModal(modalDialog(
     title = "CSV input error",
     "The CSV file uploaded is not formatted correctly. Please make sure that it contains at least 
-    the coordinates (x and y) and the depth of the sampled sites (dybde). In any doubts, refer to the README.",
+    the coordinates (x and y) and the depth of the sampled sites (torvdybde_cm). In any doubts, refer to the README.",
     easyClose = TRUE
   ))
 }
 
+print_corrupt_shp <- function(){
+  showModal(modalDialog(
+    title = "SHP input error",
+    "The SHP file uploaded is corrupt. Please make sure that it contains the polygon(s) of the area of interest. 
+          In any doubts, refer to the README. The session is reloaded.",
+    easyClose = TRUE
+  ))
+}
 
 # Error message if files are not opened properly