mhq_terr_popunits seems to miss many terrestrial (GRTS) sampling units from MHQ

[florisvdh]

It seems that quite a lot of terrestrial sampling units in MHQ are not represented in the mhq_terr_popunits data source. Maybe I'm doing something wrong, maybe something is out of date.

> mhq_terr_datapath <- file.path(dirname(gitroot), "n2khab-sample-admin/data/mhq_terr/rapportage2025")
> mhq_terr_popunits_grts <-
    read_vc("mhq_terr_popunits", root = mhq_terr_datapath) %>%
    as_tibble()
> 
> # mhq_samples_textdatapath <- "/substitute/with/correct/filepath"
> mhq_samples_heath_haymeadow <-
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_terr_cyclus2_anb.csv")
    )
> mhq_samples_mire_grassland <-                                                                                    
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_terr_cyclus2_inbo.csv")
    )
> mhq_samples_dunes <-                                                                                             
    read_csv_mhq_samples(file.path(
      mhq_samples_textdatapath,
      "mhq_terr_cyclus2_duinen_inbo_2023-05-11.csv"
    ))
> mhq_samples_forest <-                                                                                            
    read_csv_mhq_samples(
      file.path(mhq_samples_textdatapath, "mhq_forests_cycle1.csv"),
      grts_var = "grts_ranking"
    )
> mhq_terr_samples_grts <-                                                                                         
    bind_rows(
      mhq_samples_heath_haymeadow,
      mhq_samples_mire_grassland,
      mhq_samples_dunes,
      mhq_samples_forest,
    ) %>% 
    # reset some peculiarities in read_csv_mhq_samples()
    rename(type = stratum, grts_ranking_draw = grts_address) %>% 
    arrange(type, grts_ranking_draw)
> 
> # number of terrestrial MHQ sampling units
> nrow(mhq_terr_samples_grts)
[1] 1657
> 
> # data in terrestrial MHQ samples but missing from mhq_terr_popunits_grts
> 
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>% 
    nrow()
[1] 380
> 
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>% 
    count(type, sort = TRUE) %>% 
    print(n = Inf)
# A tibble: 32 × 2
   type            n
   <fct>       <int>
 1 91E0_vo        40
 2 6510_hu        39
 3 2330_bu        33
 4 2310           31
 5 2330_dw        22
 6 91E0_va        21
 7 91E0_vc        20
 8 2170           19
 9 91E0_vn        19
10 91E0_vm        16
11 6410_mo        15
12 4030           14
13 6230_hn        10
14 6510_hus       10
15 2190_overig     8
16 4010            8
17 1330_hpr        7
18 7140_meso       7
19 2130_had        5
20 2330            4
21 6230_ha         4
22 6510_hua        4
23 6510_huk        4
24 9160            4
25 2130_hd         3
26 2190_mp         3
27 9130_end        3
28 6230_hmo        2
29 7140_oli        2
30 2120            1
31 6410_ve         1
32 9130_fm         1

Definition of function read_csv_mhq_samples() (which is also set up to cater for aquatic data sources)

#' Read and tidy csv file with MHQ samples
#'
#' @param path File path.
#' @param grts_var Column name to be used as GRTS address.
#' @param single_type Optional string to set a single type that represents all
#'   rows.
read_csv_mhq_samples <- function(path,
                                 grts_var = "grts_ranking_draw",
                                 single_type = NULL) {
  (
    if (is.null(single_type)) {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer(),
          habitattype = col_character()
        )
      ) %>%
        select(
          stratum = habitattype,
          grts_address = {{ grts_var }}
        )
    } else {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer()
        )
      ) %>%
        mutate(stratum = single_type) %>%
        select(
          stratum,
          grts_address = {{ grts_var }}
        )
    }
  ) %>%
    mutate(
      # shortcut a complication for forests, where > 1 type is sometimes noted
      stratum = str_extract(stratum, "^\\w+(\\+$)?"),
      stratum = parse_factor(stratum, levels = levels(n2khab_strata$stratum))
    ) %>%
    arrange(stratum, grts_address)

[ToonHub]

I discovered 3 causes of mismatch between the selected of terrestrial sampling units in MHQ and the the mhq_terr_popunits data source.

It helps to run following code.

mhq_terr_assessments <-
    read_vc("mhq_terr_assessments", root = mhq_terr_datapath) %>%
    as_tibble()

mhq_terr_assessments_latest <- mhq_terr_assessments %>%
    group_by(point_code) %>%
    filter(assessment_date == max(assessment_date)) %>%
    ungroup() %>%
    select(point_code, type_assessment = type, is_present, assessment_date, not_measurable, inaccessible)

check_missing <- mhq_terr_samples_grts %>%
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>%
    left_join(mhq_terr_assessments_latest, join_by(point_code))

• With mhq_terr_popunits I aim to give the latest status of each sampling unit. It includes information from the field assessments (after the terrestrial sampling units were selected). So one cause of mismatch is that the assessed type is different from the type in the drawn sampling units or the assesed type is not present. The information on the type is included in mhq_terr_popunits and if the type is not present the popunit is removed. Therefore the combination grts_ranking_draw/type will not be present in the orignially selected sampling units.
• A second cause is that I forget to include the selected sampling units that are based on habitatmap 2020 and were not selected based on a previos version of habitatmap.
• Finally I found only 4 mismatches for the forest types. Maybe something went wrong with your code when reading the forest sampling units. I did not use the function read_csv_mhq_samples.

Let's discuss further how we can better organise the data to avoid this confusion.

[ToonHub]

I discovered an additional cause of missing sampling units in mhq_terr_popunits.

In mhq_terr_popunits I removed sampling units that were not yet assessed and for which te type does not correspond with habitatmap 2023. The selection of sampling units for cycle 2 was based on habitatmap 2020.

I added the sampling units that were wissing due to the second cause explained above. See this branch.

I end up with 6 sampling units that are still missing in mhq_terr_popunits. I will investigate this further.

[florisvdh]

Thanks for pursuing this, much appreciated. I see that the read_csv_mhq_samples() cannot be run without the n2khab_strata object, sorry for that. I've uploaded it as an rds file for completeness: n2khab_strata.rds.zip; after unzipping you can run below code to load it. I've also included an update of the read function in order to capture the point_code column so that your join can be made.

Updated code for read_csv_mhq_samples()

n2khab_strata <- readRDS("n2khab_strata.rds")

#' Read and tidy csv file with MHQ samples
#'
#' @param path File path.
#' @param grts_var Column name to be used as GRTS address.
#' @param single_type Optional string to set a single type that represents all
#'   rows.
read_csv_mhq_samples <- function(path,
                                 grts_var = "grts_ranking_draw",
                                 single_type = NULL) {
  (
    if (is.null(single_type)) {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer(),
          habitattype = col_character(),
          point_code = col_character()
        )
      ) %>%
        select(
          stratum = habitattype,
          grts_address = {{ grts_var }},
          point_code
        )
    } else {
      read_delim(
        file = path,
        delim = ";",
        col_types = cols_only(
          {{ grts_var }} := col_integer(),
          point_code = col_character()
        )
      ) %>%
        mutate(stratum = single_type) %>%
        select(
          stratum,
          grts_address = {{ grts_var }},
          point_code
        )
    }
  ) %>%
    mutate(
      # shortcut a complication for forests, where > 1 type is sometimes noted
      stratum = str_extract(stratum, "^\\w+(\\+$)?"),
      stratum = parse_factor(stratum, levels = levels(n2khab_strata$stratum))
    ) %>%
    arrange(stratum, grts_address)
}

Then, the code of my initial chunk can be run, followed by the chunk you provided.

I add some findings below.

> # several point codes have multiple types for the same assessment date, but this
> # may make sense if these are negative assessments
> mhq_terr_assessments_latest %>% 
    mutate(n_assessed_latest = n(), .by = point_code) %>% 
    filter(n_assessed_latest > 1) %>% 
    count(is_present)
# A tibble: 2 × 2
  is_present     n
  <lgl>      <int>
1 FALSE        241
2 NA            56
> 
> # joining the unique list of assessed point codes since we are still neglecting
> # the 'type' in the assessments data source
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>% 
    left_join(
      mhq_terr_assessments_latest %>% 
        distinct(point_code, is_present) %>% 
        mutate(assessed = TRUE),
      join_by(point_code),
      unmatched = "drop"
    ) %>% 
    count(assessed, is_present)
# A tibble: 4 × 3
  assessed is_present     n
  <lgl>    <lgl>      <int>
1 TRUE     FALSE         53
2 TRUE     TRUE          93
3 TRUE     NA             6
4 NA       NA           228
> # above, assessed 'NA' means that these MHQ sample units, missed by the popunits
> # data frame, are not present in the assessment data source (but effectively
> # drawn)
> 
> # some missed sampling units are present in the assessment data source with a
> # type match (this is a small subset of the TRUE x TRUE combination above):
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>% 
    inner_join(
      mhq_terr_assessments_latest %>% 
        distinct(point_code, type_assessment, is_present) %>% 
        filter(is_present),
      join_by(type == type_assessment, point_code),
      relationship = "one-to-one",
      unmatched = "drop"
    )
# A tibble: 7 × 4
  type        grts_ranking_draw point_code is_present
  <fct>                   <int> <chr>      <lgl>     
1 2190_overig           2200577 2200577_1  TRUE      
2 2190_overig           3700913 3700913_1  TRUE      
3 2190_overig           4761777 4761777_1  TRUE      
4 6230_ha                 34374 34374_1    TRUE      
5 6230_hn                304710 304710_1   TRUE      
6 6510_hu                857281 857281_1   TRUE      
7 6510_hu                901945 901945_1   TRUE      

[florisvdh]

For clarity I should add that I'm still using the objects referred in the beginning of this issue (not yet trying the ones from the correction branch)

Finally I found only 4 mismatches for the forest types. Maybe something went wrong with your code when reading the forest sampling units.

Regarding the forest types, I've used the file mhq_forests_cycle1.csv; in R I get below so this seems OK to me when comparing with the file:

> mhq_samples_forest %>% rename(type = stratum, grts_ranking_draw = grts_address)
# A tibble: 176 × 3
   type     grts_ranking_draw point_code
   <fct>                <int> <chr>     
 1 9130_end             10129 10129_1   
 2 9130_end             12665 12665_1   
 3 9130_end             15833 15833_1   
 4 9130_end             18321 18321_1   
 5 9130_end             20281 20281_1   
 6 9130_end             21329 21329_1   
 7 9130_end             26273 26273_1   
 8 9130_end             27865 27865_1   
 9 9130_end             29073 29073_1   
10 9130_end             30873 30873_1   
# ℹ 166 more rows
# ℹ Use `print(n = ...)` to see more rows
> mhq_samples_forest %>% rename(type = stratum, grts_ranking_draw = grts_address) %>% tail()
# A tibble: 6 × 3
  type    grts_ranking_draw point_code
  <fct>               <int> <chr>     
1 91E0_vo           1167350 1167350_1 
2 91E0_vo           1183570 1183570_1 
3 91E0_vo           1187538 1187538_1 
4 91E0_vo           1189342 1189342_1 
5 91E0_vo           1194934 1194934_1 
6 91E0_vo           1203506 1203506_1 

Wondering whether the used column to join may be at stake (either grts_ranking_draw or point_code), it appears that the number of mismatched sampling units is not the same. But I cannot reproduce the 4 forest cases you find, so we're doing something different I guess.

> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, grts_ranking_draw)) %>% 
    nrow()
[1] 380
> 
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, point_code)) %>% 
    nrow()
[1] 453
> 
> # forest types only
> mhq_terr_samples_grts %>% 
    anti_join(mhq_terr_popunits_grts, join_by(type, point_code)) %>% 
    filter(str_detect(type, "^9")) %>% 
    count(type)
# A tibble: 8 × 2
  type         n
  <fct>    <int>
1 9130_end     3
2 9130_fm      1
3 9160         4
4 91E0_va     21
5 91E0_vc     20
6 91E0_vm     16
7 91E0_vn     19
8 91E0_vo     40