Zones tampons avec attributs du plus proche voisin

Il y a quelques semaines, Christina Buelow a demandé sur Twitter comment elle pouvait créer les polygones des eaux côtières en partant des polygones terrestres, et en s’assurant de ne pas avoir de chevauchement de zone tampon. Cela permettrait de savoir quel secteur terrestre est le plus proche de n’importe quel endroit dans la mer. Explorons quelle solution est ressortie des discussions.

Créer une zone tampon autour des polygones adjacents

Comme toujours lorsqu’on essaie de trouver une solution à un problème spécifique, on commence par un exemple reproductible, communément appelé “reprex” chez les utilisateurs R. Cela signifie que nous créons le code minimal qui permet à quiconque d’exécuter le code, de faire face au même problème et de montrer tout ce que nous avons testé. Lors de la discussion sur Twitter, je n’ai pas eu assez de temps pour trouver une solution, mais j’ai préparé un gist avec mes essais, le point que j’avais atteint et, plus important encore, ce que je comptais réaliser après.

Tout d’abord, chargeons quelques packages

library(sf)
library(ggplot2)
library(dplyr)

my_blue <- "#1e73be"

Ensuite, téléchargeons des données

# Define where to save the dataset
extraWD <- "."
# Get some data available to anyone
if (!file.exists(file.path(extraWD, "departement.zip"))) {
  githubURL <- "https://github.com/statnmap/blog_tips/raw/master/2018-07-14-introduction-to-mapping-with-sf-and-co/data/departement.zip"
  download.file(githubURL, file.path(extraWD, "departement.zip"))
  unzip(file.path(extraWD, "departement.zip"), exdir = extraWD)
}

Je réduis le jeu de données pour mieux voir ce que je fais. Ici, je filtre sur la région Bretagne.
Comme vous pouvez le voir, il y a 4 départements dans cette région, donc 4 entités spatiales.

departements_l93 <- read_sf(dsn = extraWD, layer = "DEPARTEMENT")

# Reduce dataset
bretagne_l93 <- departements_l93 %>%
  filter(NOM_REG == "BRETAGNE")

bretagne_l93
## Simple feature collection with 4 features and 11 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 99217.1 ymin: 6704195 xmax: 400572.3 ymax: 6881964
## projected CRS:  RGF93_Lambert_93
## # A tibble: 4 x 12
##   ID_GEOFLA CODE_DEPT NOM_DEPT CODE_CHF NOM_CHF X_CHF_LIEU Y_CHF_LIEU X_CENTROID
## * <chr>     <chr>     <chr>    <chr>    <chr>        <int>      <int>      <int>
## 1 DEPARTEM… 35        ILLE-ET… 238      RENNES      351870    6789646     326237
## 2 DEPARTEM… 22        COTES-D… 278      SAINT-…     274947    6839206     260216
## 3 DEPARTEM… 29        FINISTE… 232      QUIMPER     171326    6789950     115669
## 4 DEPARTEM… 56        MORBIHAN 260      VANNES      267889    6744074     256790
## # … with 4 more variables: Y_CENTROID <int>, CODE_REG <chr>, NOM_REG <chr>,
## #   geometry <MULTIPOLYGON [m]>

Maintenant, je peux créer un graphique avec une zone tampon autour des polygones séparés.

# Create buffer
bretagne_buffer_l93 <- bretagne_l93 %>% 
  st_buffer(
  dist = 
    units::set_units(10, km)
) 

ggplot() +
  geom_sf(data = bretagne_buffer_l93, fill = my_blue, alpha = 0.2, colour = my_blue) +
  geom_sf(data = bretagne_l93, fill = NA) +
  theme_bw()

Vous pouvez voir qu’il y a des zones de chevauchement.
C’est là que les polygones de Voronoï vont venir nous aider. À partir d’ici, je vais prendre quelques idées que Christina a recueillies lors de ses discussions sur Twitter. Son code pour les zones tampons des zones côtières est ici.

# Create a merged region entity
bretagne_union_l93 <- bretagne_l93 %>% 
  summarise()

ggplot() + 
  geom_sf(data = bretagne_union_l93, colour = my_blue, fill = NA)

# Create a doughnut for regional seas areas, 30km off coasts
bretagne_seas_l93 <- bretagne_union_l93 %>% 
  st_buffer(
    dist = units::set_units(30, km)
  ) %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_seas_l93, colour = my_blue, fill = NA)

# Remove inside terrestrial parts
bretagne_donut_l93 <- bretagne_seas_l93 %>% 
  st_difference(bretagne_union_l93) %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_donut_l93, colour = my_blue, fill = my_blue, alpha = 0.2)

Divisons cette zone en plusieurs polygones

Nous utilisons le trait de côte comme base.

# density <- units::set_units(0.1, 1/km)

# First merge everything and transform as lines
bretagne_points_l93 <- bretagne_union_l93 %>% 
  # transform polygons as multi-lines
  st_cast("MULTILINESTRING") %>% 
  # transform multi-lines as separate unique lines
  st_cast("LINESTRING") %>% 
  # transform as points, 0.1 per km (= 1 every 10 km)
  # Choose density according to precision needed
  st_line_sample(density = units::set_units(0.1, 1/km)) %>% 
  st_sf() #%>% 
  # remove empty geometry (small islands where sample is not possible with this density)
  # filter(!st_is_empty(geometry))
  
ggplot() + 
  geom_sf(data = bretagne_points_l93, colour = my_blue, fill = NA)

# Create voronoi polygons around points
bretagne_voronoi_l93 <- bretagne_points_l93 %>% 
  st_combine() %>% 
  st_voronoi() %>% 
  st_cast()

ggplot() + 
  geom_sf(data = bretagne_voronoi_l93, colour = my_blue, fill = NA)

Jointure des polygones de voronoi et des polygones type “doughnut”

bretagne_voronoi_donut_l93 <- bretagne_voronoi_l93 %>% 
  st_intersection(bretagne_donut_l93) %>% 
  st_cast() %>% 
  st_sf()

ggplot() + 
  geom_sf(data = bretagne_voronoi_donut_l93, colour = my_blue, fill = NA)

Joindre les polygones de voronoi-donut au jeu de données original pour l’intersection avec les régions

Nous réalisons une sorte de left_join (jointure à gauche) basée sur des attributs spatiaux en utilisant st_join(). Cela permettra de conserver la géométrie originale de nos polygones de Voronoï, tout en récupérant les variables du jeu de données terrestres jointes lorsqu’il y a intersection.

# Equivalent to left_join
bretagne_vd_region_l93 <- bretagne_voronoi_donut_l93 %>% 
  st_join(bretagne_l93)

ggplot() + 
  geom_sf(data = bretagne_vd_region_l93,
          aes(colour = NOM_DEPT, fill = NOM_DEPT),
          alpha = 0.25)

# Interactive view
mapview::mapview(bretagne_vd_region_l93, zcol = "NOM_DEPT")