Update functions replace_existing_points and draw_sample

source/R/draw_sample.R

@@ -22,13 +22,18 @@
 #' @importFrom sf st_coordinates st_drop_geometry
 draw_sample <- function(
   sampling_frame,
-  sample_size = round(nrow(sampling_frame) / 20),
+  sample_size,
   sample_size_multiplication = 1,
   balance = c("X", "Y"),
   ips,
   seed = 1234,
   ...
 ) {
+  if (missing(sample_size)) {
+    sample_size <- ifelse(nrow(sampling_frame) < 2000,
+                          nrow(sampling_frame) / 5,
+                          nrow(sampling_frame) / 20)
+  }
 
   sample_size_extra <- sample_size * sample_size_multiplication
 

source/R/steekproeftrekking_nabehandeling.R

@@ -1,5 +1,5 @@
 path_to_existing <- function(file) {
-  file.path(mbag_dir, "data", "processed", file)
+  file.path(mbag_dir, "data", "SOVON", file)
 }
 
 steekproef_uitdunnen <- function(
@@ -64,16 +64,34 @@ thin_sample <- function(sample, thin_dist) {
 
 replace_by_existing <- function(sample,
                                 existing_points,
+                                id_existing_points,
+                                gebied,
                                 overlap_prop = 0.5,
                                 sbp_file) {
+  stopifnot(sf::st_crs(existing_points)$input != "EPSG:28992")
+  # Determine where HOL and OL is in gebied
+  ol <- selectie_openheid(
+    gebied = gebied,
+    ol_strata = c("OL"))
+
+  hol <- selectie_openheid(
+    gebied = gebied,
+    ol_strata = c("HOL"))
   # Recalculate sbp stratum existing points
   old_points <- existing_points %>%
+    st_transform(crs = 31370) %>%
+    st_filter(gebied) %>%
     mutate(is_sbp = st_intersects(.,
                                   st_union(sbp_file),
                                   sparse = FALSE) %>%
              as.logical()
     ) %>%
-    mutate(openheid_klasse = ifelse(grepl("HOL", stratum), "HOL", "OL")) %>%
+    mutate(openheid_klasse = ifelse(st_within(., hol) %>% as.logical(),
+                                    "HOL",
+                                    ifelse(st_within(., ol) %>% as.logical(),
+                                           "OL",
+                                           NA))) %>%
+    rename(definitief_punt = id_existing_points) %>%
     select(definitief_punt, openheid_klasse, is_sbp)
 
   # Add buffers
@@ -104,7 +122,8 @@ replace_by_existing <- function(sample,
     inner_join(intersect, by = "definitief_punt") %>%
     select(definitief_punt, pointid) %>%
     inner_join(sample %>% st_drop_geometry(), by = "pointid") %>%
-    select(pointid = definitief_punt, all_of(columns))
+    select(pointid = definitief_punt, all_of(columns)) %>%
+    mutate(pointid = as.character(pointid))
 
   # Add to sample
   sample_out <- sample %>%

source/markdown/steekproefontwerp_zandstreek.Rmd

@@ -54,7 +54,7 @@ Dat bestand kan je openen met de volgende code:
 
 ## Perimeter bepalen
 
-Om de perimeter te bepalen voer je onderstaande code uit in mas_steekproef_zn_lm.Rproj, na laden van globals in targets_steekproef_zn_lm.Rmd.
+Om de perimeter te bepalen voer je onderstaande code uit in mas_steekproef_zn_lm.Rproj, na laden van globals in targets_steekproef_zL_lm_za_kp_po.Rmd.
 
 De code maakt een intersectie tussen landbouwstreek Zandstreek (heel België) en contour Vlaanderen en voegt deze toe  aan file "zandleemstreek_perimeters.gpkg".
 
@@ -889,7 +889,9 @@ allocatie <- tar_read(allocatie_df, store = targets_store) %>%
   ))
 
 allocatie %>%
-  mutate(steekproeftrekking = round(popsize / 20)) %>%
+  mutate(steekproeftrekking = ifelse(popsize < 2000,
+                                     round(popsize / 5),
+                                     round(popsize / 20))) %>%
   select(-c(openheid_sbp, targetsize, excess)) %>%
   kable(digits = 3)
 ```
@@ -1098,52 +1100,43 @@ Bovendien zijn de data en de metadata in een standaard en open formaat, waardoor
 
 ## Overlap met bestaande telpunten
 
-We willen een lijst met alle bestaande punten en hun historie (vanaf wanneer beginnen tellen).
-Zo kunnen we een regel verzinnen om bestaande telpunten over te nemen bij een nieuwe trekking.
-Bijvoorbeeld: Als een nieuw punt voor meer dan 50 % overlapt met een bestaand punt en deze behoort tot hetzelfde stratum, dan wordt het nieuwe punt vervangen door het bestaande punt.
-
-> Opvragen bij SOVON?
-
-We bekijken overlap met de pilootstudie.
-Er is niet veel overlap.
 We vervangen alle telpunten uit de steekproef van de Zandstreek door reeds bestaande punten indien ze voor minstens 50 % overlappen en ze tot hetzelfde stratum behoren.
 Indien er meerdere bestaande punten overlappen, nemen we het bestaande punt dat het meest overlapt.
 Omdat we nu een nieuwe laag van SBP gebruiken, moeten we dit stratum opnieuw berekenen voor de bestaande punten.
 
 ```{r}
+ol <- selectie_openheid(
+  gebied = perimeters_data,
+  ol_strata = c("OL"))
+
+hol <- selectie_openheid(
+  gebied = perimeters_data,
+  ol_strata = c("HOL"))
+
 # inlezen bestaande punten
-oude_cirkels <- read_csv(here(
+oude_cirkels <- read_sf(here(
   "data",
-  "mas_piloot",
-  "steekproef_piloot_avimap.csv"
-)) %>%
-  select(
-    plotnaam = definitief_punt,
-    regio,
-    stratum,
-    X, Y
-  ) %>%
-  mutate(
-    openheid_klasse = ifelse(grepl("HOL", stratum), "HOL", "OL"),
-    sbp = ifelse(grepl("binnen", stratum), "binnen", "buiten"),
-    stratum = gsub("L\\s{1}", "L\\\n", stratum),
-    jaar = 2022L
-  ) %>%
-  select(-c(stratum, sbp)) %>%
-  st_as_sf(coords = c("X", "Y"), crs = 31370)
-
-# herbereken sbp
-oude_punten_sbp <- add_stratum_sbp(
-  punten_sf = oude_cirkels,
-  sbp = sbp_akkervogels
-) %>%
-  mutate(
-    sbp = ifelse(is_sbp == TRUE, "binnen", "buiten"),
-    openheid_sbp = paste(openheid_klasse, sbp, "plan", sep = " ")
+  "SOVON",
+  "avimap_601_0_MAS_Vlaanderen_telpunten_xy.shp"
+)) %>% 
+  st_transform(crs = 31370) %>%
+  st_filter(perimeters_data) %>%
+  mutate(is_sbp = st_intersects(.,
+                                st_union(sbp_akkervogels),
+                                sparse = FALSE) %>%
+           as.logical()
   ) %>%
-  select(-sbp)
-
-oude_cirkels_sbp <- oude_punten_sbp %>%
+  mutate(openheid_klasse = ifelse(st_within(geometry, hol, sparse = FALSE),
+                                  "HOL",
+                                  ifelse(st_within(geometry, ol, sparse = FALSE),
+                                         "OL",
+                                         NA))) %>%
+  rename(definitief_punt = id) %>%
+  select(definitief_punt, openheid_klasse, is_sbp) %>% 
+  mutate(sbp = ifelse(is_sbp == TRUE, "binnen", "buiten"),
+         openheid_sbp = paste(openheid_klasse, sbp, "plan", sep = " "))
+
+oude_cirkels_buffer <- oude_cirkels %>%
   st_buffer(300)
 
 # buffer nieuwe punten
@@ -1158,32 +1151,32 @@ steekproef_raw <- tar_read(steekproef_thinned, store = targets_store) %>%
     )
   ))
 
-cirkels_Zandstreek <- steekproef_raw %>%
+cirkels_zandstreek <- steekproef_raw %>%
   st_buffer(300)
 ```
 
 ```{r}
 # visualisatie
-mapview(oude_cirkels_sbp,
+mapview(oude_cirkels_buffer,
   col.regions = "red", layer = "bestaande punten"
 ) +
-  mapview(cirkels_Zandstreek,
+  mapview(cirkels_zandstreek,
     layer = "nieuwe punten"
   )
 ```
 
-Er zijn geen punten die overlappen voor minstens 50 % en tot hetzelfde stratum behoren. Daarom voeren we onderstaande code niet uit.
+Er zijn een aantal punten die overlappen voor minstens 50 % en tot hetzelfde stratum behoren.
 
-```{r, eval=FALSE}
+```{r}
 # welke punten overlappen minstens 50 %
-intersect <- st_intersection(oude_cirkels_sbp, cirkels_Zandstreek) %>%
+intersect <- st_intersection(oude_cirkels_buffer, cirkels_zandstreek) %>%
   mutate(intersect_area = st_area(.) %>%
            units::drop_units()) %>%
   filter(
     intersect_area >= 0.5 * 300 * 300 * pi, # 50 %
     openheid_sbp == openheid_sbp.1
   ) %>% # zelfde stratum
-  select(plotnaam, pointid, intersect_area) %>%
+  select(definitief_punt, pointid, intersect_area) %>%
   st_drop_geometry() %>%
   group_by(pointid) %>%
   filter(intersect_area == max(intersect_area)) %>%
@@ -1197,15 +1190,6 @@ intersect %>%
   )
 ```
 
-```{r, eval=FALSE}
-oude_punten_sbp %>%
-  st_drop_geometry() %>%
-  mutate(overlap = ifelse(plotnaam %in% intersect$plotnaam, "yes", "no")) %>%
-  count(regio, overlap) %>%
-  kable()
-```
-
-
 ## Output als KML
 
 We nemen een steekproef van 50 cirkels per stratum (eerste set) waarbij we ook de volgende 50 punten per stratum includeren voor als er een punt niet zou kunnen geteld worden in het veld (reserve set).