1. Create the following directory structure in the folder of your choice:
library(mapview)
library(leafpop)
library(plotly)
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(sp)
library(data.table)
library(fasttime)
library(xml2)

source("/Users/eshaahluwalia/Desktop/_/vu/year one/p4/grondwater processes/ProjectICE2/FileReadFunctions2_2026.R")

Program reads and creates a database (grwdata.dat) all hydraulic head files downloaded from www.dinoloket.nl read the elevation of the cross-section you are working on. This notebook belongs to Q1 of ICE2.

Create cross section

ahn <- read.csv("/Users/eshaahluwalia/Desktop/_/vu/year one/p4/grondwater processes/ProjectICE1/cross-section/cross2.csv", header=T)
x3 <- c(ahn[1, 1], ahn[length(ahn[, 1]), 1])
y3 <- c(ahn[1, 2], ahn[length(ahn[, 1]), 2])
Crosssection <- SpatialLines(list(Lines(Line(cbind(x3, y3)), ID = "a")))
proj4string(Crosssection) <- "EPSG:28992"
gw_dir <- "/Users/eshaahluwalia/Desktop/_/vu/year one/p4/grondwater processes/ProjectICE1/Grwdata/DINO_Grondwaterstanden"
# Create a vector of all filenames
allfilesDINO <- list.files(path = "/Users/eshaahluwalia/Desktop/_/vu/year one/p4/grondwater processes/ProjectICE1/Grwdata/DINO_Grondwaterstanden", pattern = ".csv", full.names = TRUE,recursive=T)
allfilesCSVBRO <- list.files(path = "/Users/eshaahluwalia/Desktop/_/vu/year one/p4/grondwater processes/ProjectICE1/Grwdata/BRO_Grondwatermonitoring", pattern = ".csv", full.names = TRUE,recursive=T)
allfilesXMLBRO = sapply(allfilesCSVBRO,FUN=function(x){
  parts=gregexpr("/",x)[[1]]
  dirv=substr(x,1,parts[length(parts)]-1)
list.files(path = dirv, pattern = ".xml", full.names = TRUE,recursive=T)
})

Read all groundwaterdatafiles Calculate average grwhead for each file

# Create new list to store the loop output
sumdat <- vector("list", length = length(allfilesDINO)+length(allfilesCSVBRO))
for (i in 1:length(allfilesDINO)) {
 #i = 1
  suppressWarnings({
   sumdat[[i]]=ReadDINO(allfilesDINO[i])
  })
}

for (i in 1:length(allfilesCSVBRO)) {
# i = 1
  suppressWarnings({
   sumdat[[i+length(allfilesDINO)]]=ReadBRO(allfilesXMLBRO[i],allfilesCSVBRO[i])
  })
}
# Combine them together
sumdat <- do.call(rbind, sumdat)
##clean up
sumdatc=sumdat
sumdatc=sumdatc[sumdatc$nrmet>10,] #atleast 10 days of measurements
sumdatc=sumdatc[!is.na(sumdatc$mv),] #need a surface elevation
sumdatc=sumdatc[sumdatc$ufil!=Inf,] #need a bottom filter
sumdatc=sumdatc[sumdatc$bfil!=-Inf,] #need a top filter

sumdat=sumdatc

##if the center of filters als <2 m apart we asumme them to be the same and 1 i removed
mindist=rep(NA,length(sumdat[,1]))
for(i in 1:length(sumdat[,1])){
x=sumdat[i,]
  mindist[i]=sort(sqrt((x$Xloc-sumdat$Xloc)^2+(x$Yloc-sumdat$Yloc)^2+((x$bfil+x$ufil)/2-(sumdat$bfil+sumdat$ufil)/2)^2))[2]
}
sumdat=sumdat[mindist>2,]

sumdat <- data.frame(ID = rownames(sumdat), sumdat) 
coordinates(sumdat) <- ~ Xloc + Yloc
proj4string(sumdat) <- "EPSG:28992"

Basic selection: > 100 measurements; filter length < 3m; filter length > 0.5 m; surface elevation is known Plot with all selected wells. If you click on a well you get a table with characteristics

selb <- sumdat$nrmet > 100 & 
  (sumdat$bfil - sumdat$ufil) < 3 &
  (sumdat$bfil - sumdat$ufil) > 0.5 & 
  (abs(sumdat$bfil) != Inf) & 
  (abs(sumdat$ufil) != Inf) &
  (abs(sumdat$mv) != Inf) &
  !is.na(sumdat$avnap)

mapview(Crosssection, color = "red") +
  mapview(sumdat[selb, ], zcol = "avusat",layer.name="Grw depth")

Selection criteria

Set selection for shallow, medium and deep groundwater well.

##select wells with a groundwater depth between 
GrwSel=data.frame(Shallow=c(1,2),Medium=c(8,20),Deep=c(25,60))
##and with a filter depth of x m below the avergae groundwater level
FilterSel=data.frame(Shallow=5,Medium=22,Deep=62)

1 shallow groundwater (<2 m below surface)

If you click on the wells you get the well characteristics

## detailed selection 1: Filter within 5m from average hydraulic head; average level < 2 m from surface; more than 365 measurments with daily resolution
sel1 <- (sumdat$avnap - (sumdat$bfil + sumdat$ufil) / 2) < FilterSel$Shallow &
  sumdat$avusat > GrwSel$Shallow[1] & sumdat$avusat < GrwSel$Shallow[2] & sumdat$dtd > 1200

mapview(Crosssection, color = "red") +
  mapview(sumdat[selb & sel1, ], zcol = "mv",layer.name="Elevation")

With basic plot. write down Id of favorite well=170

If you click on the wells you get the well you see the measured timeseries

mapview(Crosssection, color = "red") +
  mapview(sumdat[selb & sel1, ], zcol = "avusat", popup = popupGraph(plotlist), legend = F)

2 medium groundwater depth (>10, <20m below surface)

2: Map with average characteristics

## detailed selection 2: Filter within 10m from average hydraulic head; average
## level > 10m from surface; average level < 20m from surface; more than 365
## measurments with daily resolution
sel2 <- (sumdat$avnap - (sumdat$bfil + sumdat$ufil) / 2) < FilterSel$Medium &
  sumdat$avusat > GrwSel$Medium[1] & sumdat$avusat < GrwSel$Medium[2] & sumdat$dtd > 1825

mapview(Crosssection, color = "red") +
  mapview(sumdat[selb & sel2, ], zcol = "avusat")

2: Map with measured timeseries

mapview(Crosssection, color = "red") +
  mapview(sumdat[selb & sel2, ], zcol = "avusat", popup = popupGraph(plotlist), legend = F)

3 deep groundwater (>30m below surface)

3: Map with average characteristics

## detailed selection 3: Filter within 15m from average hydraulic head; average level > 30m from surface; average level < 20m from surface; more than 365 measurements with daily resolution
sel3 <- (sumdat$avnap - (sumdat$bfil + sumdat$ufil) / 2) < FilterSel$Deep &
  sumdat$avusat > GrwSel$Deep[1] & sumdat$avusat < GrwSel$Deep[2] & sumdat$dtd > 1825

mapview(Crosssection, color = "red") + mapview(sumdat[selb & sel3, ], zcol = "avusat")

3: Map with measured timeseries

mapview(Crosssection, color = "red") + 
  mapview(sumdat[selb & sel3, ], zcol = "avusat", popup = popupGraph(plotlist), legend = F)

Make your own selection fom the selection tables above!!!! 1 of each selection

####!!!Important: add you own IDs!!!
## depth1, depth2, depth3
dsel <- c(170, 192, 178)
finaldat <- lapply(match(dsel,type.convert(sumdat$ID,as.is=T)), FUN = function(x) {
  readgrwhead(x, sumdat)
})

Plot with selected wells

dates <- range(do.call("c", lapply(finaldat, FUN = function(x) x[[3]])))
xlim <- as.POSIXct(range(do.call("c", lapply(finaldat, FUN = function(x) x[[3]]))))
ylim <- range(unlist(lapply(finaldat, FUN = function(x) x[[2]])), na.rm = T)
fd <- unlist(lapply(finaldat, FUN = function(x) {
  round(x[[1]]$mv - (x[[1]]$bfil + x[[1]]$ufil) / 2, 1)
}))
avg <- unlist(lapply(finaldat, FUN = function(x) round(x[[1]]$avusat, 1)))

plot_ly() %>% 
  add_lines(x = as.POSIXct(finaldat[[1]][[3]]),
            y = finaldat[[1]][[2]] / 100,
            color = paste("Filter depth = ", fd[1], "; Average Grw depth = ", avg[1])) %>% 
  add_lines(x = as.POSIXct(finaldat[[2]][[3]]),
            y = finaldat[[2]][[2]] / 100,
            color = paste("Filter depth = ", fd[2], "; Average Grw depth = ", avg[2])) %>% 
  add_lines(x = as.POSIXct(finaldat[[3]][[3]]),
            y = finaldat[[3]][[2]] / 100,
            color = paste("Filter depth = ", fd[3], "; Average Grw depth = ", avg[3])) %>% 
  # Axis limits
  layout(xaxis = list(title = "",
                      range = xlim), 
         yaxis = list(title = "Hydraulic head [m +NAP]",
                      range = ylim / 100),
                  legend = list(orientation = 'h'))
ylim <- range(unlist(lapply(finaldat, FUN = function(x) {
  x[[2]] - mean(x[[2]], na.rm = T)
})), na.rm = T)

plot_ly() %>% 
  add_lines(x = as.POSIXct(finaldat[[1]][[3]]),
            y = (finaldat[[1]][[2]] - mean(finaldat[[1]][[2]], na.rm = T)) / 100,
            color = paste("Filter depth = ", fd[1], "; Average Grw depth = ", avg[1])) %>% 
  add_lines(x = as.POSIXct(finaldat[[2]][[3]]),
            y = (finaldat[[2]][[2]] - mean(finaldat[[2]][[2]], na.rm = T)) / 100,
            color = paste("Filter depth = ", fd[2], "; Average Grw depth = ", avg[2])) %>% 
  add_lines(x = as.POSIXct(finaldat[[3]][[3]]),
            y = (finaldat[[3]][[2]] - mean(finaldat[[3]][[2]], na.rm = T)) / 100,
            color = paste("Filter depth = ", fd[3], "; Average Grw depth = ", avg[3])) %>% 
  # Axis limits
  layout(xaxis = list(title = "",
                      range = xlim), 
         yaxis = list(title = "Hydraulic head [m +NAP]",
                      range = ylim / 100),
         legend = list(orientation = 'h'))

Save data

dst <- 'output'

if(!dir.exists(dst)) dir.create(dst, recursive = TRUE)

for (i in 1:length(finaldat)) {
  write.table(cbind(format(finaldat[[i]][[3]],"%Y-%m-%d"),
                    finaldat[[i]][[2]] / 100,
                    finaldat[[i]][[1]]$mv,
                    (finaldat[[i]][[1]]$bfil + finaldat[[i]][[1]]$ufil) / 2),
              file = file.path(dst, paste("pb", i, ".txt", sep = "")),
              sep = ",",
              row.names = F,
              col.names = c("Date", "Head NAP", "Surface", "Filterdepth"))
}