mhq_terr_popunits: some other suspicious values in grts_ranking_draw

[florisvdh]

Using current master; d29cd9b.

Related to #47 (non-existing values in grts_ranking_draw), I started wondering whether other faulty GRTS addresses might exist that do match with the GRTSmaster_habitats data source by mere accident, hence go unnoticed.

Assuming that the cases where grts_ranking_draw != grts_ranking reflect local replacements and consequently the involved cells should be spatially near each other, it was not too difficult to find more suspicious cases.

> grts_mh <- read_GRTSmh()
> 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()
> 
> existing_grts_addresses_suspected <- c(
    1, 22, 54, 73, 82, 86, 105
  )
> 
> mhq_terr_popunits_grts_suspect <- 
    mhq_terr_popunits_grts %>% 
    filter(grts_ranking_draw %in% existing_grts_addresses_suspected) %>% 
    select(
      grts_original = grts_ranking_draw,
      grts_new = grts_ranking
    )
> 
> # making a subset of grts_mh with only the above used GRTS addresses
> grts_mh_subset <- mask(
    grts_mh,
    grts_mh,
    maskvalues = c(
      mhq_terr_popunits_grts_suspect$grts_original,
      mhq_terr_popunits_grts_suspect$grts_new
    ),
    inverse = TRUE
  )
> # making a spatial index of grts_mh_subset before extracting cell center coordinates
> grts_mh_subset_index <- tibble(
    cellid = seq_len(ncell(grts_mh_subset)),
    grts_address = values(grts_mh_subset)[, 1]
  ) %>%
    filter(!is.na(grts_address))  
> 
> # determine cell center coordinates and pairwise distances
> mhq_terr_popunits_grts_suspect %>% 
    inner_join(grts_mh_subset_index, join_by(grts_original == grts_address)) %>% 
    rename(cellid_original = cellid) %>% 
    inner_join(grts_mh_subset_index, join_by(grts_new == grts_address)) %>% 
    rename(cellid_new = cellid) %>% 
    mutate(
      geometry_original_x = xyFromCell(grts_mh_subset, cellid_original)[, "x"],
      geometry_original_y = xyFromCell(grts_mh_subset, cellid_original)[, "y"],
      geometry_new_x = xyFromCell(grts_mh_subset, cellid_new)[, "x"],
      geometry_new_y = xyFromCell(grts_mh_subset, cellid_new)[, "y"]
    ) %>% 
    select(-starts_with("cellid")) %>% 
    st_as_sf(coords = c("geometry_original_x", "geometry_original_y"), crs = "EPSG:31370") %>% 
    as_tibble() %>% 
    rename(geometry_original = geometry) %>% 
    st_as_sf(coords = c("geometry_new_x", "geometry_new_y"), crs = "EPSG:31370") %>% 
    as_tibble() %>% 
    rename(geometry_new = geometry) %>% 
    mutate(
      distance = st_distance(geometry_original, geometry_new, by_element = TRUE) %>% 
        set_units("km")
    ) %>% 
    arrange(desc(distance))
# A tibble: 7 × 5
  grts_original grts_new   geometry_original        geometry_new distance
          <dbl>    <dbl>         <POINT [m]>         <POINT [m]>     [km]
1            82 10578737 (248125.6 195342.1) (37629.59 207374.1)   211.  
2            86 10654897 (185629.6 201422.1) (23933.59 198830.1)   162.  
3            54 10288565 (198109.6 191566.1) (59453.59 221422.1)   142.  
4            22 10082997 (165981.6 199086.1) (65213.59 224302.1)   104.  
5           105 12644145 (122941.6 202670.1) (32861.59 203694.1)    90.1 
6            73 10419637 (98525.59 217486.1) (59005.59 221230.1)    39.7 
7             1   103425 (59933.59 217230.1) (55165.59 218574.1)     4.95