`mapvotr` alternative to Etalab shapes, using `mapvotr` directly

# 1. Get the data

```r
library(arrow)
library(happign)
library(tidyverse)

# code based on
# https://github.com/InseeFrLab/mapvotr/blob/main/data-raw/address_sample.R

city_id <- "59350"

# https://paul-carteron.github.io/happign/reference/get_wfs.html
# (slightly changed since {mapvotr} came out)
contours_com <- happign::get_wfs(x = NULL,
layer = "BDTOPO_V3:commune",
# apikey = "topographie",
ecql_filter = str_c("code_insee = ", city_id))

contours_com_sample <- contours_com %>%
dplyr::filter(code_insee %in% city_id) %>%
dplyr::select(code_insee, geometry)

readr::write_rds(contours_com_sample, "contours_com_sample-lille.rds")

addresses_sample <- "table-adresses-reu.parquet" %>%
arrow::read_parquet() %>%
dplyr::filter(code_commune_ref %in% city_id) %>%
dplyr::select(id_brut_bv_reu,
code_commune_ref,
longitude,
latitude,
geo_type,
geo_score,
nb_adresses)

readr::write_rds(addresses_sample, "addresses_sample-lille.rds")
```

# 2. Get contours

```r
library(mapvotr)
library(sf)
library(tidyverse)

# code based on
# https://inseefrlab.github.io/mapvotr/articles/mapvotr.html

contours_com_sample <- readr::read_rds("contours_com_sample-lille.rds")
stopifnot(sf::st_is_valid(contours_com_sample))
sf::st_crs(contours_com_sample) # "WGS 84"

addresses_sample <- readr::read_rds("addresses_sample-lille.rds")
glimpse(addresses_sample)

# 1. prepare addresses
prep_adr <- mapvotr::prepare_address(
address = addresses_sample,
contours_com = contours_com_sample,
var_cog1 = "code_commune_ref",
var_cog2 = "code_insee",
var_bv1 = "id_brut_bv_reu",
path_log = NULL)

# 2. build contours

#### Hypothèses de modélisation
MIN_POINT_COM <- 50
MIN_ADDRESS_BV <- 15
MIN_ADDRESS_SHOOT <- 5

# Création des contours de Douarnenez
cog <- "59350"
lcontours <- mapvotr::create_contours(
prep_adr,
cog,
min_points_com = MIN_POINT_COM,
min_address_bv = MIN_ADDRESS_BV,
min_address_shoot = MIN_ADDRESS_SHOOT,
var_cog1 = "code_commune_ref",
var_cog2 = "code_insee",
var_bv1 = "id_brut_bv_reu",
var_geo_score = "geo_score",
var_nbaddress = "nb_adresses",
path_log = NULL)

stopifnot(sf::st_is_valid(lcontours$contours))
stopifnot(sf::st_is_valid(lcontours$contours_simplified))

sf::st_crs(lcontours$contours) # "RGF93 v1 / Lambert-93"
sf::st_crs(lcontours$contours_simplified) # "RGF93 v1 / Lambert-93"

ggplot(lcontours$contours) +
geom_sf(color = "red") +
geom_sf(data = lcontours$contours_simplified,
fill = NA, color = "blue")

readr::write_rds(lcontours, "bv-lille-mapvotr.rds")
```

# 3. Test areal interpolation and compare

Comparison with aggregated file shown to have an issue in #3

```r
library(areal)
library(sf)
library(tidyverse)

lcontours <- readr::read_rds("bv-lille-mapvotr.rds")

# compare to aggregated file (see issue #3)
# https://github.com/briatte/selection-bv/issues/3
bv_cedricr <- readr::read_rds("https://f.briatte.org/temp/bv-lille.rds")
sf::st_crs(bv_cedricr)
table(sf::st_is_valid(bv_cedricr)) # 3 FALSE, 123 TRUE
# this flips the validity?! 123 FALSE
table(sf::st_is_valid(sf::st_set_crs(bv_cedricr, "WGS84")))
# back to 3 FALSE 123 TRUE
table(sf::st_is_valid(sf::st_transform(sf::st_set_crs(bv_cedricr, "WGS84"), 2154)))

# completely different bbox
sf::st_bbox(lcontours$contours_simplified)
sf::st_bbox(bv_cedricr)

# strictly equivalent when plotting, though
ggplot(lcontours$contours) +
geom_sf(data = lcontours$contours_simplified,
fill = NA, color = "blue", size = 3) +
geom_sf(data = sf::st_set_crs(bv_cedricr, "WGS84"),
fill = NA, color = "yellow", size = 1.5)

# now onto the areal interpolation test

bv <- lcontours$contours_simplified
sf::st_crs(bv) # "RGF93 v1 / Lambert-93"

iris <- readr::read_rds("https://f.briatte.org/temp/iris-lille.rds")
sf::st_crs(iris) # "RGF93 Lambert 93"

# looks like we're good to go
areal::ar_validate(iris, bv, varList = "test_variable", verbose = TRUE)

# everything's TRUE, interpolation should work

interp1 <- areal::aw_interpolate(bv, source = iris,
tid = "id_brut_bv_reu", sid = "CODE_IRIS",
extensive = "test_variable",
weight = "total", output = "tibble")

# ... and it does

interp1
# # A tibble: 126 × 2
# id_brut_bv_reu test_variable
#
# 1 59350_1 1.13
# 2 59350_10 1.29
# 3 59350_101 3.53
# 4 59350_102 1.01
# 5 59350_104 0.700
# 6 59350_105 0.419
# 7 59350_106 1.13
# 8 59350_107 2.29
# 9 59350_108 0.528
# 10 59350_109 0.751

# compare to results with aggregated file
bv4 <- st_make_valid(bv_cedricr) %>%
sf::st_set_crs("WGS84") %>%
sf::st_transform(sf::st_crs(iris))

# .. and NOW this works
interp2 <- areal::aw_interpolate(bv4, source = iris,
tid = "codeBureauVote", sid = "CODE_IRIS",
extensive = "test_variable",
weight = "total", output = "tibble") %>%
select(codeBureauVote, test_variable)

interp2
# # A tibble: 126 × 2
# codeBureauVote test_variable
#
# 1 59350_0001 1.13
# 2 59350_0010 1.29
# 3 59350_0101 3.53
# 4 59350_0102 1.01
# 5 59350_0104 0.700
# 6 59350_0105 0.419
# 7 59350_0106 1.13
# 8 59350_0107 2.29
# 9 59350_0108 0.528
# 10 59350_0109 0.751

# very very minor differences
interp1 %>%
mutate(id = str_remove(id_brut_bv_reu, "59350_") %>%
str_pad(width = 4, pad = "0") %>%
str_c("59350_", .)) %>%
select(-id_brut_bv_reu) %>%
full_join(interp2, by = c("id" = "codeBureauVote")) %>%
mutate(diff = test_variable.x - test_variable.y) %>%
filter(diff > 0.001)
```