Note: this is a bookdown project, supposed to be run from within the src/update_shallowgroundwater subfolder. You can use the update_shallowgroundwater.Rproj RStudio project file in this subfolder to run it.

1 Introduction

1.1 Current, intermediate state of the shallowgroundwater data source

A shallowgroundwater data source was prepared by Dries Adriaens and colleagues in order to delineate the areas in Flanders where the Mean Lowest Watertable (‘GLG’) is located higher (less deep) than approx. 1.5 - 2 m below soil surface. This data source is needed to restrict the target population of several MNE monitoring schemes in the groundwater compartiment. It is more specifically applied to restrict the ‘locally groundwater dependent’ types, i.e. to discard the sites without shallow groundwater.

Based on an evaluation that uses the ‘(almost) everywhere groundwater dependent’ types, the current state still mismatches a proportion of those types. Hence, the area it covers needs to be made broader.

1.2 Aim and method

In this document, we append several layers to shallowgroundwater. The aim is to select areas which shallowgroundwater doesn’t cover yet and where ‘(almost) everywhere groundwater dependent’ types (here also called ‘obliggwdep’ types) are present. Such ‘areas of interest’ are e.g. anthropogenic soil type polygons and polygons with ‘dunes’ as soil texture type.

  • In the case of anthropogenic and dune soil polygons, we will subseqently limit those areas to buffers (of 100 m) around habitatmap polygons with ‘obliggwdep’ types – still subsequently clipping along the borders of the aforementioned areas. So in practice, we first select the habitatmap polygons with ‘obliggwdep’ types that intersect the anthropogenic and dune soil polygons, we buffer them, and we clip the result along the areas of interest.
    • The buffers around habitatmap polygons with ‘obliggwdep’ types are applied since it is expected that shallow groundwater (with potentially ‘locally groundwater dependent’ types) will occur in the vicinity of those polygons.

The areas of interest are:

  • all anthropogenic soil type polygons: i.e. where the bsm_mo_soilunitype column of the soilmap_simple data source starts with "O";
  • all dune texture polygons: i.e. where the bsm_mo_tex column equals "X".

Regarding the anthropogenic soil type polygons, from earlier evaluation it appeared that the narrow ones that contain ‘obliggwdep’ types are interesting as a whole. Hence we select them as a whole instead of selecting buffers around the habitatmap polygons with ‘obliggwdep’ types.

  • An appropriate algorithm had to be made that selects meaningful polygons - at first it resulted in complete cities meeting the narrowness requirement (since they are part of anthropogenic soil type polygons with long, narrow parts).
  • Hence an additional condition was required: a minimal amount of ‘obliggwdep’ type must be present in such polygons.

Finally, we stitch the parts. The resulting new version of shallowgroundwater is written here as shallowgroundwater.gpkg.

Hence we will have three separate sources, which we will document as three extra TRUE/FALSE column inside shallowgroundwater.gpkg:

  • anthrop_gwdep
  • narrowanthrop_gwdep
  • dunes_gwdep

Where an area overlapped between several of these sources, more than one column is set as TRUE.

We do most geoprocessing tasks using sf (mostly using GEOS as backend). However in some cases we use qgisprocess (QGIS as backend) when the QGIS algorithm is either faster than, or absent from, sf.

2 Preparation

2.1 Data setup and checks

local_root <- find_root(has_file("update_shallowgroundwater.Rproj"))
datapath <- file.path(local_root, "data")
if (!dir.exists(datapath)) dir.create(datapath)
sgpath <- find_root_file("n2khab_data/10_raw/shallowgroundwater", 
                            criterion = has_dir("n2khab_data"))

if (!dir.exists(sgpath)) dir.create(sgpath)
inputpath <- file.path(tempdir(), "shallowgroundwater")
dir.create(inputpath, showWarnings = FALSE)
drive_auth(email = TRUE)
## → Using an auto-discovered, cached token
##   To suppress this message, modify your code or options to clearly consent to the use of a cached token
##   See gargle's "Non-interactive auth" vignette for more details:
##   <https://gargle.r-lib.org/articles/non-interactive-auth.html>
## → The googledrive package is using a cached token for 'floris.vanderhaeghe@inbo.be'
# download "shallowgroundwater_20211129.gpkg.zip", which only contains the 
# "ZOG_20211129_diss" layer from ZonesOndiepGrondwater_20211129.gpkg":
as_id("1bCQUMm1s0iDwr6HWe12bWQpBRiY-YS3H") %>% 
  drive_download(file.path(inputpath, "shallowgroundwater.gpkg.zip"), 
                 overwrite = TRUE)
## Auto-refreshing stale OAuth token.
## File downloaded:
##   * shallowgroundwater_20211129.gpkg.zip
## Saved locally as:
##   * /tmp/RtmpbAYnwN/shallowgroundwater/shallowgroundwater.gpkg.zip
unzip(file.path(inputpath, "shallowgroundwater.gpkg.zip"),
      exdir = inputpath)
unlink(file.path(inputpath, "shallowgroundwater.gpkg.zip"))

Verification of several input data sets.

Checksums of shallowgroundwater (using default algorithm of n2khab::checksum()):

checksums <- 
  list.files(inputpath, full.names = TRUE) %>% 
  checksum
checksums %>% 
  as.matrix
##                         [,1]              
## shallowgroundwater.gpkg "29660a8672412c95"
if (!file.exists(file.path(inputpath, "shallowgroundwater.gpkg"))) {
  stop("The input GeoPackage file is missing from ", inputpath)
}
if (any(checksums != 
        c(shallowgroundwater.gpkg = "29660a8672412c95"))
) stop("Beware, your version of shallowgroundwater is not the one for which this code was intended!")
sspath <- find_root_file("n2khab_data/20_processed/soilmap_simple/soilmap_simple.gpkg", 
                            criterion = has_dir("n2khab_data"))
if (!file.exists(sspath)) stop("Please organize soilmap_simple location. See e.g. https://inbo.github.io/n2khab/articles/v022_example.html")
if (checksum(sspath) != "9f4204b476506031") stop("Incorrect soilmap_simple version.")

hmtpath <- find_root_file("n2khab_data/20_processed/habitatmap_terr/habitatmap_terr.gpkg", 
                            criterion = has_dir("n2khab_data"))
if (!file.exists(hmtpath)) stop("Please organize habitatmap_terr location. See https://inbo.github.io/n2khab")
if (checksum(hmtpath) != "72d3144016e816ed") stop("Incorrect habitatmap_terr version.")

Reading data sets.

soilmap <- read_soilmap()
sg <- read_sf(file.path(inputpath, "shallowgroundwater.gpkg"), 
              layer = params$input_layername,
              crs = 31370)

2.2 Preparing derived input data

Constructing the geospatial and attribute data of habitatmap_terr polygons that contain ‘obliggwdep’ types, relevant in the MNE monitoring programme.

scheme_types_mne_expanded <-
  read_scheme_types(extended = TRUE) %>%
  filter(programme == "MNE") %>%
  distinct(scheme, type) %>%
  group_by(scheme) %>%
  expand_types %>%
  ungroup
types_obliggwdep <-
  read_types() %>% 
  filter(groundw_dep == "GD2") %>%
  select(type) %>% 
  # following step does 2 things: 
  # - drop 6 (main) type codes that don't appear in habitatmap_terr
  # - avoid 7110, which we don't recognize within MNE (one location in habitatmap_terr)
  semi_join(scheme_types_mne_expanded, by = "type")
hmt <- read_habitatmap_terr(keep_aq_types = FALSE)
hmt_occ_obliggwdep <- 
  hmt$habitatmap_terr_types %>% 
  semi_join(types_obliggwdep,  
            by = "type")
hmt_pol_obliggwdep <- 
  hmt$habitatmap_terr_polygons %>% 
  select(polygon_id, description_orig) %>% 
  semi_join(hmt_occ_obliggwdep,
            by = "polygon_id")
st_agr(hmt_pol_obliggwdep) <- "constant"

Creating a geospatial object of anthropogenic soil type polygons, and calculate a property ‘thinness’ (from here) – lower means more narrow – as well as the areal fraction occupied by habitatmap_terr polygons having ‘obliggwdep’ types.

calculate_obliggwdep_occupation <- function(x) {
  st_intersection(x, hmt_pol_obliggwdep) %>% 
  mutate(subarea = st_area(.)) %>% 
  st_drop_geometry %>% 
  group_by(bsm_poly_id) %>% 
  summarise(obliggwdep_pol_frac = drop_units(sum(subarea) / first(area)))
}
soil_anthrop <- 
  soilmap %>% 
  filter(str_sub(bsm_mo_soilunitype, 1, 1) == "O") %>% 
  select(bsm_poly_id, bsm_mo_soilunitype) %>% 
  mutate(perimeter = lwgeom::st_perimeter(.),
         area = st_area(.),
         thinness =  4 * pi * (area / perimeter ^2) %>% drop_units()) %>% 
  st_make_valid()
st_agr(soil_anthrop) <- "constant"
soil_anthrop_occupied <- 
  soil_anthrop %>% 
  calculate_obliggwdep_occupation
soil_anthrop <- 
  soil_anthrop %>% 
  left_join(soil_anthrop_occupied,
            by = "bsm_poly_id")

The distribution of both calculated variables:

soil_anthrop %>% 
  st_drop_geometry %>% 
  ggplot(aes(x = thinness)) + 
  geom_histogram(binwidth = 0.01, fill = "white", colour = "grey70")

soil_anthrop %>% 
  st_drop_geometry %>% 
  ggplot(aes(x = obliggwdep_pol_frac)) + 
  geom_histogram(binwidth = 0.01, fill = "white", colour = "grey70")

Creating a geospatial object of dune texture polygons, with the areal fraction occupied by habitatmap_terr polygons having ‘obliggwdep’ types.

soil_dunes <- 
  soilmap %>% 
  filter(bsm_mo_tex == "X") %>% 
  select(bsm_poly_id, bsm_mo_soilunitype) %>% 
  mutate(area = st_area(.)) %>% 
  st_make_valid()
st_agr(soil_dunes) <- "constant"
soil_dunes_occupied <- 
  soil_dunes %>% 
  calculate_obliggwdep_occupation
soil_dunes <- 
  soil_dunes %>% 
  left_join(soil_dunes_occupied,
            by = "bsm_poly_id")

The distribution of that areal fraction:

soil_dunes %>% 
  st_drop_geometry %>% 
  ggplot(aes(x = obliggwdep_pol_frac)) + 
  geom_histogram(binwidth = 0.01, fill = "white", colour = "grey70")

Cleaning shallowgroundwater geometry.

sg_repairedgeom <- 
  st_make_valid(sg) %>% 
  st_cast("MULTIPOLYGON")
write_sf(sg_repairedgeom,
         file.path(datapath, 
                   "sg_repairedgeom.gpkg"),
         delete_dsn = TRUE)
sg_repairedgeom_union <- 
  sg_repairedgeom %>% 
  st_union %>% 
  st_as_sf
write_sf(sg_repairedgeom_union,
         file.path(datapath, 
                   "sg_repairedgeom_union.gpkg"),
         delete_dsn = TRUE)

3 Execution

3.1 Generating the areas of interest (where shallow groundwater is more probable)

3.1.1 Buffer 100 m around obliggwdep types, within soil_anthrop

We create this geospatial object and also write it as a GPKG file, since the calculation takes several minutes.

# following step takes several minutes!
buffered_obliggwdep_within_soilanthrop <-
  hmt_pol_obliggwdep[soil_anthrop, ] %>% 
  st_make_valid(.) %>% 
  st_buffer(params$buffer_obliggwdep) %>% 
  st_union() %>% 
  st_intersection(soil_anthrop)
write_sf(buffered_obliggwdep_within_soilanthrop, 
         file.path(datapath, 
                   "buffered_obliggwdep_within_soilanthrop.gpkg"),
         delete_dsn = TRUE)

Resulting object: buffered_obliggwdep_within_soilanthrop.

3.1.2 Buffer 100 m around obliggwdep types, within soil_dunes

We create this geospatial object and also write it as a GPKG file, since the calculation takes several minutes.

# following step takes several minutes!
buffered_obliggwdep_within_soildunes <-
  hmt_pol_obliggwdep[soil_dunes, ] %>% 
  st_make_valid(.) %>% 
  st_buffer(params$buffer_obliggwdep) %>% 
  st_union() %>% 
  st_intersection(soil_dunes)
write_sf(buffered_obliggwdep_within_soildunes, 
         file.path(datapath, 
                   "buffered_obliggwdep_within_soildunes.gpkg"),
         delete_dsn = TRUE)

Resulting object: buffered_obliggwdep_within_soildunes.

3.1.3 Narrow anthropogenic soil polygons with ‘obliggwdep’ types

After several iterations, the below used conditions on thinness and obliggwdep_pol_frac appeared satisfying.

narrow_soilanthrop_with_obliggwdep <- 
  soil_anthrop %>% 
  filter(thinness < 0.15,
         obliggwdep_pol_frac > 0.04) %>% 
  st_cast("MULTIPOLYGON")

Interactive map:

mapview(narrow_soilanthrop_with_obliggwdep,
        zcol = "thinness",
        color = "pink3",
        lwd = 2)

Resulting object: narrow_soilanthrop_with_obliggwdep.

3.2 Combine the resulting layers in new_areas.gpkg

plan(multisession) # future kept here as demonstration
if (gpkg_newareas_exists) unlink(file.path(datapath, "new_areas.gpkg"))

list(buffered_obliggwdep_within_soilanthrop =
       buffered_obliggwdep_within_soilanthrop %>% 
       st_union %>% 
       st_as_sf,
     buffered_obliggwdep_within_soildunes =
       buffered_obliggwdep_within_soildunes %>% 
       st_union %>% 
       st_as_sf,
     narrow_soilanthrop_with_obliggwdep =
       narrow_soilanthrop_with_obliggwdep %>% 
       st_union %>% 
       st_as_sf) %>% 
  walk2(names(.),
        ~write_sf(.x,
                  dsn = file.path(datapath, "new_areas.gpkg"),
                  layer = .y))

We write the resulting layers into new_areas.gpkg:

st_layers(file.path(datapath, "new_areas.gpkg"))
## Driver: GPKG 
## Available layers:
##                               layer_name geometry_type features fields
## 1 buffered_obliggwdep_within_soilanthrop Multi Polygon        1      0
## 2   buffered_obliggwdep_within_soildunes Multi Polygon        1      0
## 3     narrow_soilanthrop_with_obliggwdep Multi Polygon        1      0

Each layer consists of just one Multipolygon. That is the result of unioning steps (obtained with st_union(); corresponding QGIS terminology is native:aggregate). The unioning is done in order to get the simplest possible intersection result in the next step.

3.3 Union the new zones (keeping intersections)

We extend shallowgroundwater as described in section 1.2.

# qgis_arguments("native:union")
# use of future unneeded in this case (runs quickly)
f <- future({
  elapsed <- 
    system.time(
      result <- qgis_run_algorithm("native:union",
                                   INPUT = 
                                     buffered_obliggwdep_within_soilanthrop %>% 
                                     mutate(anthrop_gwdep = TRUE),
                                   OVERLAY = 
                                     narrow_soilanthrop_with_obliggwdep %>% 
                                     mutate(narrowanthrop_gwdep = TRUE))
    )
  list(result = result, elapsed = elapsed)
}, seed = NULL)
# resolved(f)
value(f, stdout = FALSE, signal = FALSE)$elapsed
##    user  system elapsed 
##   9.213   0.512  10.490
res_filepath <- value(f, stdout = FALSE)$result$OUTPUT %>% as.character
## Argument `OVERLAY_FIELDS_PREFIX` is unspecified (using QGIS default value).
## Using `OUTPUT = qgis_tmp_vector()`
## qt5ct: using qt5ct plugin
## Problem with SAGA installation: SAGA was not found or is not correctly installed
f <- future({
  elapsed <- 
    system.time(
      result <- qgis_run_algorithm("native:union",
                                   INPUT = res_filepath,
                                   OVERLAY = 
                                     buffered_obliggwdep_within_soildunes %>% 
                                     mutate(dunes_gwdep = TRUE))
    )
  list(result = result, elapsed = elapsed)
}, seed = NULL)
# resolved(f)
value(f, stdout = FALSE, signal = FALSE)$elapsed
##    user  system elapsed 
##   6.126   0.330   7.207
res_filepath <- value(f, stdout = FALSE)$result$OUTPUT %>% as.character
## Argument `OVERLAY_FIELDS_PREFIX` is unspecified (using QGIS default value).
## Using `OUTPUT = qgis_tmp_vector()`
## qt5ct: using qt5ct plugin
## Problem with SAGA installation: SAGA was not found or is not correctly installed
f <- future({
  elapsed <- 
    system.time(
      result <- 
        if(params$union_sg) {
          qgis_run_algorithm("native:union",
                             INPUT = 
                               res_filepath,
                             OVERLAY = 
                               sg_repairedgeom_union %>% 
                               mutate(sg = TRUE))
        } else {
          qgis_run_algorithm("native:union",
                             INPUT = 
                               sg_repairedgeom %>% 
                               select(-OBJECTID, 
                                      -Shape_Length,
                                      -Shape_Area),
                             OVERLAY = 
                               res_filepath)
        }
    )
  list(result = result, elapsed = elapsed)
}, seed = NULL)
# resolved(f)
value(f, stdout = FALSE, signal = FALSE)$elapsed
##     user   system  elapsed 
## 1686.300    3.900 1687.585
res_filepath2 <- value(f, stdout = FALSE)$result$OUTPUT %>% as.character
## Argument `OVERLAY_FIELDS_PREFIX` is unspecified (using QGIS default value).
## Using `OUTPUT = qgis_tmp_vector()`
## qt5ct: using qt5ct plugin
## Problem with SAGA installation: SAGA was not found or is not correctly installed

The result is written as shallowgroundwater.gpkg.

sg_extended <-
  read_sf(res_filepath2) %>% 
  select(-starts_with("fid")) %>% 
  {if (params$union_sg) . else {
    select(.,
           geomorph_wcoast = K_GMK_, 
           anthrop_gwdep,
           narrowanthrop_gwdep,
           drainage = BOD_DRA, 
           dunes_gwdep,
           peat_profile = BOD_PROV, 
           peat_substr = BOD_SUB, 
           peat_parentmat = BOD_MOEV, 
           peat_texture = BOD_TEX, 
           phys_system = FYS_, 
           zwin = K_ZWIN_, 
           habitat_1130 = HAB_1130, 
           gwdepth_coast = K_GWD_, 
           gwdepth_local = GWM_, 
           seepage = BDS_, 
           peat_survey = VEK_, 
           duneslack = K_DUINV_)
  }} %>% 
  mutate(across(where(is.integer), ~ifelse(is.na(.) | . == 0L, FALSE, TRUE))) %>%
  mutate(across(where(is.logical), ~ifelse(is.na(.), FALSE, .)))
system.time(
st_write(sg_extended,
         file.path(sgpath, output_filename),
         layer = "shallowgroundwater",
         delete_dsn = TRUE,
         quiet = TRUE)
)
##    user  system elapsed 
##   1.129   0.465   2.190
plan(sequential)

4 Check results

filepath <- file.path(sgpath, output_filename)
  • Checksums:
c("xxh64", "md5", "sha256") %>% 
  map_dfr(~tibble(algorithm = ., 
                  checksum = checksum(filepath, .))) %>% 
  kable
algorithm checksum
xxh64 702a7fb1981d2c18
md5 5f60ca125dd535bb64b60d60c1701f3d
sha256 9280ff0059fb77f1115fcbdaa03efd17925ab65d0e36b3b0b9178a48a3022cc3
  • Contents:
(pol <- read_sf(filepath))
## Simple feature collection with 307 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 22273.46 ymin: 153000 xmax: 258633.9 ymax: 244024.8
## Projected CRS: Belge 1972 / Belgian Lambert 72
## # A tibble: 307 x 18
##    geomorph_wcoast anthrop_gwdep narrowanthrop_gwdep drainage dunes_gwdep
##    <lgl>           <lgl>         <lgl>               <lgl>    <lgl>      
##  1 FALSE           TRUE          TRUE                FALSE    FALSE      
##  2 FALSE           TRUE          FALSE               FALSE    FALSE      
##  3 FALSE           FALSE         TRUE                FALSE    FALSE      
##  4 FALSE           TRUE          FALSE               FALSE    FALSE      
##  5 TRUE            TRUE          TRUE                FALSE    FALSE      
##  6 TRUE            FALSE         TRUE                FALSE    FALSE      
##  7 TRUE            FALSE         TRUE                FALSE    FALSE      
##  8 TRUE            TRUE          TRUE                FALSE    FALSE      
##  9 TRUE            TRUE          FALSE               FALSE    FALSE      
## 10 TRUE            FALSE         TRUE                FALSE    FALSE      
## # … with 297 more rows, and 13 more variables: peat_profile <lgl>,
## #   peat_substr <lgl>, peat_parentmat <lgl>, peat_texture <lgl>,
## #   phys_system <lgl>, zwin <lgl>, habitat_1130 <lgl>, gwdepth_coast <lgl>,
## #   gwdepth_local <lgl>, seepage <lgl>, peat_survey <lgl>, duneslack <lgl>,
## #   geom <MULTIPOLYGON [m]>
glimpse(pol)
## Rows: 307
## Columns: 18
## $ geomorph_wcoast     <lgl> FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE…
## $ anthrop_gwdep       <lgl> TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRUE,…
## $ narrowanthrop_gwdep <lgl> TRUE, FALSE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, …
## $ drainage            <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ dunes_gwdep         <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ peat_profile        <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ peat_substr         <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ peat_parentmat      <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ peat_texture        <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ phys_system         <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ zwin                <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ habitat_1130        <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ gwdepth_coast       <lgl> TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE,…
## $ gwdepth_local       <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ seepage             <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ peat_survey         <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, F…
## $ duneslack           <lgl> FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, TRUE, FAL…
## $ geom                <MULTIPOLYGON [m]> MULTIPOLYGON (((35286.71 20..., MULTI…
pol %>% 
  st_drop_geometry %>% 
  summary
##  geomorph_wcoast anthrop_gwdep   narrowanthrop_gwdep  drainage      
##  Mode :logical   Mode :logical   Mode :logical       Mode :logical  
##  FALSE:281       FALSE:178       FALSE:200           FALSE:130      
##  TRUE :26        TRUE :129       TRUE :107           TRUE :177      
##  dunes_gwdep     peat_profile    peat_substr     peat_parentmat 
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:283       FALSE:248       FALSE:239       FALSE:296      
##  TRUE :24        TRUE :59        TRUE :68        TRUE :11       
##  peat_texture    phys_system        zwin         habitat_1130   
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:266       FALSE:178       FALSE:300       FALSE:235      
##  TRUE :41        TRUE :129       TRUE :7         TRUE :72       
##  gwdepth_coast   gwdepth_local    seepage        peat_survey    
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:272       FALSE:302       FALSE:211       FALSE:255      
##  TRUE :35        TRUE :5         TRUE :96        TRUE :52       
##  duneslack      
##  Mode :logical  
##  FALSE:291      
##  TRUE :16
  • Are all attribute combinations unique? This is a desirable property: each multipolygon stands for another combination of data origins.
pol %>% 
  st_drop_geometry %>% 
  {nrow(.) == nrow(distinct(.))}
## [1] TRUE
  • Is the CRS conform EPSG:31370?
 st_crs(pol) == st_crs(31370) 
## [1] TRUE
  • Sensible results?
include_graphics(image)
_Appended new areas, in a QGIS screenshot. Zoomed into a specific area. Light colour = original data source, dark colour = appended area._

Figure 4.1: Appended new areas, in a QGIS screenshot. Zoomed into a specific area. Light colour = original data source, dark colour = appended area.

5 Used environment

  • version: R version 4.1.2 (2021-11-01)
  • os: Linux Mint 20
  • system: x86_64, linux-gnu
  • ui: X11
  • language: nl_BE:nl
  • collate: nl_BE.UTF-8
  • ctype: nl_BE.UTF-8
  • tz: Europe/Brussels
  • date: 2021-12-14
  • GEOS: 3.9.0
  • GDAL: 3.2.1
  • PROJ: 7.2.1
  • GDAL_with_GEOS: true
  • USE_PROJ_H: true
  • QGIS: 3.22.1-Białowieża
Table 5.1: Loaded R packages
package loadedversion date source
askpass 1.1 2019-01-13 CRAN (R 4.1.0)
assertthat 0.2.1 2019-03-21 CRAN (R 4.1.0)
base64enc 0.1-3 2015-07-28 CRAN (R 4.1.0)
bookdown 0.22 2021-04-22 CRAN (R 4.1.0)
brew 1.0-6 2011-04-13 CRAN (R 4.1.0)
class 7.3-19 2021-05-03 CRAN (R 4.1.0)
classInt 0.4-3 2020-04-07 CRAN (R 4.1.0)
cli 2.5.0 2021-04-26 CRAN (R 4.1.0)
codetools 0.2-18 2020-11-04 CRAN (R 4.1.0)
colorspace 2.0-1 2021-05-04 CRAN (R 4.1.0)
crayon 1.4.1 2021-02-08 CRAN (R 4.1.0)
crosstalk 1.1.1 2021-01-12 CRAN (R 4.1.0)
curl 4.3.1 2021-04-30 CRAN (R 4.1.0)
DBI 1.1.1 2021-01-15 CRAN (R 4.1.0)
digest 0.6.27 2020-10-24 CRAN (R 4.1.0)
dplyr 1.0.6 2021-05-05 CRAN (R 4.1.0)
e1071 1.7-7 2021-05-23 CRAN (R 4.1.0)
ellipsis 0.3.2 2021-04-29 CRAN (R 4.1.0)
evaluate 0.14 2019-05-28 CRAN (R 4.1.0)
fansi 0.5.0 2021-05-25 CRAN (R 4.1.0)
farver 2.1.0 2021-02-28 CRAN (R 4.1.0)
forcats 0.5.1 2021-01-27 CRAN (R 4.1.0)
fs 1.5.0 2020-07-31 CRAN (R 4.1.0)
future 1.21.0 2020-12-10 CRAN (R 4.1.0)
gargle 1.1.0 2021-04-02 CRAN (R 4.1.0)
generics 0.1.0 2020-10-31 CRAN (R 4.1.0)
ggplot2 3.3.3 2020-12-30 CRAN (R 4.1.0)
git2r 0.28.0 2021-01-10 CRAN (R 4.1.0)
git2rdata 0.3.1 2021-01-21 CRAN (R 4.1.0)
globals 0.14.0 2020-11-22 CRAN (R 4.1.0)
glue 1.4.2 2020-08-27 CRAN (R 4.1.0)
googledrive 1.0.1 2020-05-05 CRAN (R 4.1.0)
gtable 0.3.0 2019-03-25 CRAN (R 4.1.0)
highr 0.9 2021-04-16 CRAN (R 4.1.0)
htmltools 0.5.1.1 2021-01-22 CRAN (R 4.1.0)
htmlwidgets 1.5.3 2020-12-10 CRAN (R 4.1.0)
httr 1.4.2 2020-07-20 CRAN (R 4.1.0)
jsonlite 1.7.2 2020-12-09 CRAN (R 4.1.0)
KernSmooth 2.23-20 2021-05-03 CRAN (R 4.1.0)
knitr 1.33 2021-04-24 CRAN (R 4.1.0)
labeling 0.4.2 2020-10-20 CRAN (R 4.1.0)
lattice 0.20-44 2021-05-02 CRAN (R 4.1.0)
leafem 0.1.6 2021-05-24 CRAN (R 4.1.0)
leaflet 2.0.4.1 2021-01-07 CRAN (R 4.1.0)
leaflet.providers 1.9.0 2019-11-09 CRAN (R 4.1.0)
leafpop 0.1.0 2021-05-22 CRAN (R 4.1.0)
lifecycle 1.0.0 2021-02-15 CRAN (R 4.1.0)
listenv 0.8.0 2019-12-05 CRAN (R 4.1.0)
lubridate 1.7.10 2021-02-26 CRAN (R 4.1.0)
lwgeom 0.2-6 2021-04-02 CRAN (R 4.1.0)
magrittr 2.0.1 2020-11-17 CRAN (R 4.1.0)
mapview 2.10.0 2021-06-05 CRAN (R 4.1.0)
munsell 0.5.0 2018-06-12 CRAN (R 4.1.0)
n2khab 0.5.0 2021-07-27 Github ()
openssl 1.4.4 2021-04-30 CRAN (R 4.1.0)
parallelly 1.27.0 2021-07-19 CRAN (R 4.1.0)
pillar 1.6.1 2021-05-16 CRAN (R 4.1.0)
pkgconfig 2.0.3 2019-09-22 CRAN (R 4.1.0)
plyr 1.8.6 2020-03-03 CRAN (R 4.1.0)
png 0.1-7 2013-12-03 CRAN (R 4.1.0)
processx 3.5.2 2021-04-30 CRAN (R 4.1.0)
proxy 0.4-26 2021-06-07 CRAN (R 4.1.0)
ps 1.6.0 2021-02-28 CRAN (R 4.1.0)
purrr 0.3.4 2020-04-17 CRAN (R 4.1.0)
qgisprocess 0.0.0.9000 2021-08-06 Github ()
R6 2.5.0 2020-10-28 CRAN (R 4.1.0)
rappdirs 0.3.3 2021-01-31 CRAN (R 4.1.0)
raster 3.4-10 2021-05-03 CRAN (R 4.1.0)
Rcpp 1.0.6 2021-01-15 CRAN (R 4.1.0)
renv 0.13.2 2021-03-30 CRAN (R 4.1.0)
rlang 0.4.11 2021-04-30 CRAN (R 4.1.0)
rmarkdown 2.8 2021-05-07 CRAN (R 4.1.0)
rprojroot 2.0.2 2020-11-15 CRAN (R 4.1.0)
rstudioapi 0.13 2020-11-12 CRAN (R 4.1.0)
satellite 1.0.2 2019-12-09 CRAN (R 4.1.0)
scales 1.1.1 2020-05-11 CRAN (R 4.1.0)
sessioninfo 1.1.1 2018-11-05 CRAN (R 4.1.0)
sf 0.9-8 2021-03-17 CRAN (R 4.1.0)
sp 1.4-5 2021-01-10 CRAN (R 4.1.0)
stringi 1.6.2 2021-05-17 CRAN (R 4.1.0)
stringr 1.4.0 2019-02-10 CRAN (R 4.1.0)
svglite 2.0.0 2021-02-20 CRAN (R 4.1.0)
systemfonts 1.0.2 2021-05-11 CRAN (R 4.1.0)
tibble 3.1.2 2021-05-16 CRAN (R 4.1.0)
tidyr 1.1.3 2021-03-03 CRAN (R 4.1.0)
tidyselect 1.1.1 2021-04-30 CRAN (R 4.1.0)
units 0.7-1 2021-03-16 CRAN (R 4.1.0)
utf8 1.2.1 2021-03-12 CRAN (R 4.1.0)
uuid 0.1-4 2020-02-26 CRAN (R 4.1.0)
vctrs 0.3.8 2021-04-29 CRAN (R 4.1.0)
webshot 0.5.2 2019-11-22 CRAN (R 4.1.0)
withr 2.4.2 2021-04-18 CRAN (R 4.1.0)
xfun 0.23 2021-05-15 CRAN (R 4.1.0)
yaml 2.2.1 2020-02-01 CRAN (R 4.1.0)