L’écosystème spatial sur R

library(sf)

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

#importer
mtq <- read_sf("data/mtq/martinique.shp")
mtq <- st_read("data/mtq/martinique.shp")

## Reading layer `martinique' from data source `C:\Users\Kim Antunez\Documents\Projets_R\quantilille\lecture\data\mtq\martinique.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 23 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
## Projected CRS: WGS 84 / UTM zone 20N

#exporter
write_sf(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)
st_write(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)

st_crs(mtq)

## Coordinate Reference System:
##   User input: WGS 84 / UTM zone 20N 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 20N",
##     BASEGEOGCRS["WGS 84",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 20N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-63,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     ID["EPSG",32620]]

mtq_4326 <- mtq %>% st_transform(4326)

plot(mtq)

## Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
## all

plot(st_geometry(mtq))

mtq_c <- st_centroid(mtq) 

## Warning in st_centroid.sf(mtq): st_centroid assumes attributes are constant over
## geometries of x

plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]

## Units: [m]
##           [,1]     [,2]      [,3]      [,4]      [,5]
## [1,]     0.000 35297.56  3091.501 12131.617 17136.310
## [2,] 35297.557     0.00 38332.602 25518.913 18605.249
## [3,]  3091.501 38332.60     0.000 15094.702 20226.198
## [4,] 12131.617 25518.91 15094.702     0.000  7177.011
## [5,] 17136.310 18605.25 20226.198  7177.011     0.000

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

mtq_u2 <- mtq %>% 
  group_by(STATUT) %>% 
  summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u2), add=T, lwd=2, border = "red", col=NA)

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

# create a polygon
m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), 
           c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

pts <- st_sample(x = mtq, size = 50)
plot(st_geometry(mtq))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

inter <- st_intersects(mtq, pts)

mtq$nbpts <- sapply(X = inter, FUN = length)
plot(st_geometry(mtq))
# display munucipalities that intersect at least 2 point
plot(st_geometry(mtq[mtq$nbpts>2,]), col = "grey", add=TRUE)
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

mtq_counts <- mtq %>% st_join(st_as_sf(pts)) %>% count(INSEE_COM)
plot(mtq_counts %>% select(n))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)

mtq_v <- st_collection_extract(st_voronoi(x = st_union(mtq_c)))
mtq_v <- st_intersection(mtq_v, st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c)
plot(st_geometry(mtq_v), col='lightblue')

Préparer / récupérer des données

Dans ce premier exemple, les données sont stockées dans un fichier shapefile.

library(sf)
library(dplyr)
# Import de la couche géographique (iris de Paris)
iris.75 <- st_read("data/iris_75.shp", stringsAsFactors = F) 

## Reading layer `iris_75' from data source `C:\Users\Kim Antunez\Documents\Projets_R\quantilille\lecture\data\iris_75.shp' using driver `ESRI Shapefile'
## Simple feature collection with 992 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
## Projected CRS: RGF93 / Lambert-93

Regardons la projection utilisée :

st_crs(iris.75)

## Coordinate Reference System:
##   User input: RGF93 / Lambert-93 
##   wkt:
## PROJCRS["RGF93 / Lambert-93",
##     BASEGEOGCRS["RGF93",
##         DATUM["Reseau Geodesique Francais 1993",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4171]],
##     CONVERSION["Lambert-93",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",46.5,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",3,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",49,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",44,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",700000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",6600000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["France - onshore and offshore, mainland and Corsica."],
##         BBOX[41.15,-9.86,51.56,10.38]],
##     ID["EPSG",2154]]

Dans ce second exemple, il s’agit de données ponctuelles stockées dans un fichier csv avec deux colonnes (latitude et longitude en WGS84). Dans ce cas, on importe le csv puis on convertit le data.frame en data.frame spatial (sf) grâce à la fonction st_as_sf. Il suffit de spécifier le nom des colonnes contenant les coordonnées ainsi que le CRS. Généralement, on utilisera le code epsg de WGS84 (4326).

# Import du dataset  
accidents.2019.paris <- readRDS("data/accidents2019_paris.RDS")
# Transformation en objet sf
accidents.2019.paris <- st_as_sf(accidents.2019.paris,
                                coords = c("long", "lat"),
                                crs = 4326, agr = "constant") %>% 
  st_transform(2154)
plot(st_geometry(accidents.2019.paris))

Regardons rapidement ces deux jeux de données.

Tout d’abord, localisons les personnes accidentées selon la gravité des leurs blessures

plot(st_geometry(iris.75))
# Les personnes non blessées ou blessées légèrement
plot(accidents.2019.paris %>% filter(grav%in%c(1,4)) %>% st_geometry,
     pch = 20, col = "darkgreen", add=TRUE, cex = 0.5)
# Les personnes blessées gravement
plot(accidents.2019.paris %>% filter(grav == 3) %>% st_geometry,
     pch = 20, col = "orange", add=TRUE, cex = 0.5)
# Les personnes tuées
plot(accidents.2019.paris %>% filter(grav == 2) %>% st_geometry,
     pch = 20, col = "red", add=TRUE, cex = 1)

Comptons par iris :

1. le nombre de personnes accidentées (nbacc) ;
2. le nombre de personnes accidentées qui ont été gravement blessées ou tuées (nbacc_blessgravtues).

inter <- st_intersects(iris.75, accidents.2019.paris)
inter_blessgravtues <- st_intersects(iris.75, accidents.2019.paris
                              %>% filter(grav%in%c(2,3)))
iris.75$nbacc <- sapply(X = inter, FUN = length)
iris.75$nbacc_blessgravtues <- sapply(X = inter_blessgravtues, FUN = length)

#Remarque : Il manque 24 accidents
nrow(accidents.2019.paris)-sum(iris.75$nbacc)

## [1] 24

library(osmdata)
# Récupérer les routes principales grâce à osm
bb      <- iris.75 %>% st_transform(4326) %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)
qm      <- add_osm_feature (q, key = 'highway',
                            value = 'motorway', value_exact = FALSE)
qt      <- add_osm_feature (q, key = 'highway',
                            value = 'trunk', value_exact = FALSE)
qp      <- add_osm_feature (q, key = 'highway',
                            value = 'primary', value_exact = FALSE)

motorway<- osmdata_sf(qm)
trunk   <- osmdata_sf(qt)
primary <- osmdata_sf(qp)

roads    <- c(primary,trunk,motorway)$osm_lines %>%
  st_transform(st_crs(iris.75))
roads.geom = st_geometry(roads)

# Récupérer le shape de la seine 
qr <- q %>% 
  add_osm_feature (key = 'waterway') %>% 
  add_osm_feature(key = "name:fr", value = "La Seine")
river <- osmdata_sf(qr)

river.geom <- c(st_geometry(river$osm_lines),
                st_geometry(river$osm_multilines)) %>%
  st_transform(st_crs(iris.75))

# Export road and river layers to shapefile
st_write(roads%>% select(name,osm_id), dsn = "data/osmdata/roadsfull.gpkg")
st_write(roads.geom, dsn = "data/osmdata/road.shp")
st_write(river.geom, dsn = "data/osmdata/river.shp")

# bbox est utilisé pour centrer sur Paris
bb <- st_bbox(iris.75)
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "ivory")
plot(st_geometry(iris.75), col = "ivory", border = "ivory3", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
plot(st_geometry(roads.geom), col="#666666", lwd = 1.2, add = TRUE)
plot(st_geometry(river.geom), col="#87cdde", lwd = 3, add = TRUE)
plot(st_geometry(accidents.2019.paris %>% filter(grav == 3 )) , pch = 20,
     col = "orange", add=TRUE, cex = 1)
plot(st_geometry(accidents.2019.paris %>% filter(grav == 2)) , pch = 20,
     col = "red", add=TRUE, cex = 1)

# Avec les coordonnées présentes dans la base de données
geo_bdd <- accidents.2019.paris %>% 
  filter(catv %in% c("VAE", "EDP à moteur")) %>% slice(1:10)

# Avec banR
# remotes::install_github("joelgombin/banR")
library(banR)
geo_banR <- accidents.2019.paris %>% 
  filter(catv %in% c("VAE", "EDP à moteur")) %>% slice(1:10) %>% 
  geocode_tbl(adresse = voie,code_insee = com) %>%
  select(latitude,longitude) %>%
  st_as_sf(coords = c("longitude", "latitude"),
           crs = 4326, agr = "constant") %>%
  st_transform(2154)

## Writing tempfile to...C:\Users\KIMANT~1\AppData\Local\Temp\RtmpU5gt83\file375c7c7b521.csv

## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.

## If file is larger than 8 MB, it must be splitted
## Size is : 666 bytes

## SuccessOKSuccess: (200) OK

## New names:
## * geometry -> geometry...10
## * geometry -> geometry...13

# Avec tidygeocoder
# install.packages("tidygeocoder")
library(tidygeocoder)

## 
## Attaching package: 'tidygeocoder'

## The following objects are masked from 'package:banR':
## 
##     geocode, reverse_geocode

geo_tidygeocoder <- accidents.2019.paris %>% 
  filter(catv %in% c("VAE", "EDP à moteur")) %>% slice(1:10) %>%
  mutate(addr = paste(voie, ", Paris, France")) %>%  
  geocode(addr,method = "osm") %>% select(lat, long) %>%
  st_as_sf(coords = c("long", "lat"),crs = 4326, agr = "constant") %>%
  st_transform(2154)

## Distances et distance moyenne entre les trois types de géocodage
# Entre celui de la base et celui de banR
st_distance(geo_bdd, geo_banR, by_element = TRUE)

## Units: [m]
##  [1] 262.55263  43.85397  19.52109  68.25663  26.75033 656.61369  50.37249
##  [8] 261.84111 675.63757 675.63757

mean(st_distance(geo_bdd, geo_banR, by_element = TRUE))

## 274.1037 [m]

# Entre celui de la base et celui de tidygeocoder
st_distance(geo_bdd,geo_tidygeocoder, by_element = TRUE)

## Units: [m]
##  [1] 734.059556 861.005957 107.911108 113.374186  23.289554  31.108417
##  [7]   8.484897  43.502021  19.184942  19.184942

mean(st_distance(geo_bdd, geo_tidygeocoder,by_element = TRUE))

## 196.1106 [m]

# Entre celui de la banR et celui de tidygeocoder
st_distance(geo_banR, geo_tidygeocoder, by_element = TRUE)

## Units: [m]
##  [1] 950.52141 820.56843 122.00239 181.62935  47.91376 642.02833  41.90755
##  [8] 305.16959 664.76010 664.76010

mean(st_distance(geo_banR, geo_tidygeocoder, by_element = TRUE))

## 444.1261 [m]

Faire des cartes interactives

Plusieurs solutions existent pour faire des cartes interactives avec R. mapview, leaflet et mapdeck sont les principales. Par simplicité, nous nous concentrons ici sur mapview.

Les cartes interactives ne sont pas forcément très pertinentes pour représenter des informations géostatistiques. En revanche, elles sont utiles pour explorer les bases de données. Voyons un exemple avec mapview concernant les accidents mortels à Paris en 2019.

#remotes::install_github("r-spatial/mapview")
library(mapview)
mapviewOptions(fgb = FALSE) #pour marcher avec le format .Rmd
# construire une carte avec certaines options pour les cercles
# avec mapview la taille des cercles reste constante quel que soit le zoom. 
# grav = 2 : individus tués
individus_tues <- accidents.2019.paris %>%
  filter(grav == 2) %>%
  mutate(age=2019-an_nais) # On ajoute leur âge pour un traitement ultérieur
mapview(individus_tues)

Quand on clique sur un point, la valeur des différentes variables de la base de données apparaissent. Cela peut aider à l’exploration de la base de données.

On customise un peu…

mapview(individus_tues,
        map.types = "Stamen.TonerLite", legend = TRUE,
        cex = 5, col.regions = "#217844", lwd = 0, alpha = 0.9,
        layer.name = 'Individus tués')

On customise encore un peu plus…

Toutefois, ajouter une légende pour la taille des ronds proportionnels ne peut pas être fait facilement.

mapview(individus_tues,
        map.types = "Stamen.TonerLite", legend=TRUE,
        layer.name = 'Individus tués',
        cex="age", zcol="sexe", lwd=0, alpha=0.9
       )

Faire des cartes statiques

Là encore, différents packages R sont utilisés pour faire des cartes statiques :

• ggplot2 est un package très utilisé pour faire tous types de graphiques, et a été adapté spécifiquement aux cartes (geom_sf).
• Le package tmap contient des fonctionnalités avancées basées sur la logique de ggplot2
• mapsf (anciennement cartography) s’appuie sur un langage dit “base R” et permet de faire des représentations cartographiques, basiques comme avancées.

Par simplicité, nous nous concentrons ici sur ggplot2, package très renommé pour tous types de graphiques.

library(ggplot2)
ggplot(accidents.2019.paris) +
  geom_bar(aes(x = catv,group = sexe,fill = sexe))

catv_ol <- accidents.2019.paris %>% st_drop_geometry %>% 
  count(catv) %>% arrange(n) %>% pull(catv)

gg <- accidents.2019.paris %>% 
  mutate(catv_o = factor(catv,levels = catv_ol)) %>%
  filter(catv_o %in% tail(catv_ol, 10))

ggplot()+geom_bar(data = gg,aes(y = catv_o, group = sexe, fill = sexe))+
  scale_fill_brewer("Sexe", palette = "Set1")+
  theme_bw()+
  labs(title = "Nombre d'accidentés par type de véhicule et sexe",
       subtitle = "à Paris en 2019, pour les hommes et les femmes ",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",
       x = "", y = "")

gg = accidents.2019.paris %>% 
  st_drop_geometry %>% 
  filter(catv %in% tail(catv_ol,9)) %>% 
  count(catv,lum,sexe) %>% 
  add_count(catv,wt=n,name="tot") %>% 
  mutate(prop = n/tot)

ggplot(gg)+geom_point(aes(y = lum, x = sexe, color = prop, size = prop))+
  facet_wrap(~catv)+
  scale_color_distiller(palette = "Reds", direction = 1)+
  labs(title = "Part d'accidentés par type de véhicule et éclairage",
       subtitle = "à Paris en 2019. ",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",
       x = "", y = "")

library(tidyr)
catv_sel = c("Bicyclette", "VL seul", "VAE", "VU seul",
             "EDP à moteur", "Scooter < 50 cm3")

gg <- accidents.2019.paris %>% 
  st_drop_geometry %>% 
  filter(catv %in% catv_sel) %>% 
  count(catv,sexe) %>% pivot_wider(names_from = "sexe",values_from="n")

ggplot(gg)+
  geom_segment(aes(x = "Homme",
                   y = `Masculin`,
                   xend = "Femme",
                   yend = `Féminin`,
                   color = catv))+
  geom_text(data=gg %>% filter(catv!="Scooter < 50 cm3"),
            aes(x = "Homme",
                y = `Masculin`,
                label = catv,
                color = catv),hjust="left") +
  geom_text(data=gg %>% filter(catv!="VU seul"),
            aes(x = "Femme",
                y = `Féminin`,
                label = catv,
                color=catv),hjust="right") +
  scale_color_discrete(guide ="none") + 
  theme_bw()+
  labs(title = "Nombre d'accidentés suivant le sexe et le type de véhicule", 
       subtitle = "à Paris en 2019",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",
       x = "", y = "")

ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_sf(data = iris.75 %>%  st_centroid(),
          aes(size= nbacc), colour="#E84923CC", show.legend = 'point') +
  scale_size(name = "Nombre d'accidents",
             breaks = c(1,10,100,200),
             range = c(0,5)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA),
        plot.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Nombre d'accidents de la route à Paris par iris",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",x="",y="")

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

library(RColorBrewer) #pour les couleurs des palettes

# Préparation des données
acc <- iris.75 %>% 
  st_join(accidents.2019.paris) %>% 
  group_by(INSEE_COM, do_union = TRUE) %>% 
  summarize(nb_acc = n(),
            nb_vl = sum(if_else(catv == "VL seul", 1, 0),
                        na.rm = TRUE),
            nb_edp = sum(if_else(catv == "EDP à moteur", 1, 0),
                         na.rm = TRUE),
            nb_velo = sum(if_else(catv == "Bicyclette", 1, 0),
                          na.rm = TRUE)) 

## `summarise()` has grouped output by 'INSEE_COM'. You can override using the `.groups` argument.

# Choix des breaks
# (quintiles de la part des accidents ayant eu lieu à vélo)
bks <- round(quantile(100*acc$nb_velo/acc$nb_acc,
                         na.rm=TRUE,
                         probs=seq(0,1,0.2)),1)

# Intégration dans la base de données
acc <- acc %>% mutate(txaccvelo = 100*nb_velo/nb_acc,
                     txaccvelo_cat = cut(txaccvelo,bks)) 

# Carte
ggplot() +
  geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
  geom_sf(data = acc, aes(fill = txaccvelo_cat)) +
  geom_sf(data = river.geom, colour = "#87cdde",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  scale_fill_brewer(name = "Part (En %)",
                    palette = "Reds",
                    na.value = "grey80") +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(iris.75)[c(1,3)],
           ylim = st_bbox(iris.75)[c(2,4)]) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "ivory",color=NA),
        plot.background = element_rect(fill = "ivory",color=NA)) +
  labs(title = "Part des Accidentés à vélos",
       subtitle = "par arrondissement à Paris en 2019",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",
       x = "", y = "")

Un exemple de traitement géomatique avancé

L’objectif de ce traitement est de compter le nombre d’accidents par tronçons de 100m sur le périphérique parisien.

Pour commencer, nous allons extraire des données OSM le squelette du périphérique. Pour ce faire nous allons filtrer le data.frame roads à partir d’une liste de noms. Nous aurions pu nous servir d’une petite carte interactive pour trouver cette sélection.

library(dplyr)
library(sf)
library(tidygraph) # equivalent dplyr pour les graphes

## 
## Attaching package: 'tidygraph'

## The following object is masked from 'package:stats':
## 
##     filter

library(sfnetworks) # sf + graphes
library(ggplot2)

# A partir de roads, on garde les routes qui correspondent au périphérique
periph_simple <- roads %>% 
  filter(!is.na(name))%>% 
  filter(name %in% c("Boulevard Périphérique Intérieur", "Pont Masséna",
                     "Tunnel Lilas","Pont Amont","Pont Aval")) %>% 
  select(name)

plot(periph_simple %>% st_geometry(),
     col=1:nrow(periph_simple),lwd=4)

C’est un bon début, mais il reste un certains nombre de problèmes :

• Les lignes extraites ne font pas 100m ;
• Il reste les deux voies parallèles au niveau des tunnels

Pour résoudre ces problèmes, nous allons nous servir de la librairie sfnetwork. Elle permet de marier des objets sf avec une description de leur topologie sous forme de graphe et elle s’appuie pour cela sur la librairie tidygraph. Cette description va nous être bien utile pour fusionner toute les lignes qui se touchent en une grande ligne. Commençons par transformer notre data.frame de lignes en un réseau spatial :

# sfnetwork permet de gérer les réseaux géospatiaux.
# On transforme periph_simple en ce type d'objet
net  = as_sfnetwork(st_geometry(periph_simple))
plot(net)

Nous allons maintenant pouvoir supprimer de ce réseau les noeuds inutiles i.e. ceux qui n’ont que deux voisins :

# to_spatial_smooth : on enlève les pseudos-noeuds en 
# préservant la connectivité du réseau
nets = convert(net,to_spatial_smooth)
plot(nets)

Cela commence à ressembler à quelque chose. Mais il reste encore quelques liens isolés :

plot(nets  %>% 
       activate(edges) %>%
       filter(st_length(x)<units::as_units(1000,"m"))
     )

Pour les supprimer, nous allons calculer la longueur de chaque lien et ne conserver que ceux qui nous permettent de construire une voie continue autour du périphérique :

# On calcule la longueur des edges, on trie par ordre décroissant
nets = nets %>% 
  activate(edges) %>% #on travaille sur les liens (et non les noeuds)
  mutate(length=st_length(x)) %>% 
  arrange(desc(length)) %>% 
  mutate(lid=1:n()) %>% 
  filter(lid %in% c(1,2,3,5)) #je ne garde que les grands accès
nets = nets %>% 
  activate(nodes) %>% 
  mutate(deg=centrality_degree()) %>% 
  filter(deg!=0)
nets

## # A sfnetwork with 4 nodes and 4 edges
## #
## # CRS:  RGF93 / Lambert-93 
## #
## # A directed simple graph with 1 component with spatially explicit edges
## #

## Registered S3 method overwritten by 'cli':
##   method     from         
##   print.boxx spatstat.geom

## # Node Data:     4 x 3 (active)
## # Geometry type: POINT
## # Dimension:     XY
## # Bounding box:  xmin: 647435.4 ymin: 6858209 xmax: 655197.6 ymax: 6864941
##                    x .tidygraph_node_index   deg
##          <POINT [m]>                 <int> <dbl>
## 1 (649843.4 6858209)                    49     1
## 2 (649423.8 6858362)                    94     1
## 3 (655197.6 6858709)                   114     1
## 4 (647435.4 6864941)                   137     1
## #
## # Edge Data:     4 x 6
## # Geometry type: LINESTRING
## # Dimension:     XY
## # Bounding box:  xmin: 645125.9 ymin: 6857527 xmax: 657088 ymax: 6866979
##    from    to                                  x .tidygraph_edge~   length   lid
##   <int> <int>                   <LINESTRING [m]> <list>                [m] <int>
## 1     4     3 (647435.4 6864941, 647451.5 68649~ <dbl [77]>       18283.7~     1
## 2     2     4 (649423.8 6858362, 649345.7 68583~ <dbl [57]>       10334.9~     2
## 3     3     1 (655197.6 6858709, 655036.5 68586~ <dbl [19]>        5927.2~     3
## # ... with 1 more row

plot(nets)

Nous pouvons maintenant récupérer uniquement la géométrie des liens du réseau, et construire une ligne unique couvrant tout le périphérique :

#On transforme le tout en lignes
lines.geom = nets %>% activate(edges) %>% st_geometry()
points_ordered = lines.geom[c(2,1,3)] %>% st_cast("POINT")
points_ordered = c(points_ordered,points_ordered[1])
line.geom = points_ordered %>% st_combine() %>% st_cast("LINESTRING")
line = st_as_sf(line.geom,id=1)
plot(line)

Il reste a découper cette longue ligne en tronçons de 500m. Pour cela, nous allons commencer par créer un ensemble de points distants de 500m le long de la ligne et uniformiser l’échantillonnage de la ligne :

# On prend une ligne et on met un point tous les 100 mètres
points_eqd = line.geom  %>% st_line_sample(density = 1/10) %>%
  st_cast("POINT")
lines_eqd = line.geom  %>% st_line_sample(density = 1/10) %>%
  st_cast("LINESTRING")
split_points = points_eqd[seq(1,length(points_eqd),by=50)]
plot(points_eqd,pch=20,cex=0.2)

Ceci va nous permettre de découper la ligne en tronçons de taille identique :

# Je découpe ma ligne avec mes points
troncons.col = lwgeom::st_split(lines_eqd,split_points)
troncons.geom = st_collection_extract(troncons.col,type = "LINESTRING")
troncons = st_sf(troncons.geom,id=1:length(troncons.geom))
plot(troncons)

st_length(troncons)

## Units: [m]
##  [1] 480.0390 500.0405 500.0396 500.0283 500.0361 500.0372 500.0316 500.0050
##  [9] 500.0274 499.9622 499.9695 499.9866 499.9942 499.9227 500.0254 499.9795
## [17] 499.9905 500.0196 500.0173 499.9714 500.0056 500.0218 500.0266 500.0263
## [25] 499.9816 500.0230 500.0182 500.0366 500.0180 500.0406 500.0405 500.0405
## [33] 500.0406 500.0332 500.0399 500.0355 500.0244 500.0359 499.9975 499.9851
## [41] 500.0402 500.0391 500.0127 500.0176 499.9421 500.0071 499.9993 500.0339
## [49] 500.0400 500.0180 499.9863 500.0232 500.0357 500.0337 499.9849 500.0223
## [57] 499.9858 500.0151 500.0281 500.0301 499.9669 500.0406 499.9588 499.9910
## [65] 499.8952 499.9658 500.0035 499.9061 500.0388 500.0317

Le plus dur est fait. Il ne reste plus qu’à compter et à faire une carte :). Notez le endCapStyle sur le st_buffer pour ne pas rajouter de marge aux deux extrémités des tronçons :

# On fait la jointure avec les accidents
periph_count= troncons  %>% st_buffer(100,endCapStyle = 'FLAT') %>% 
  st_join(accidents.2019.paris %>% filter(!duplicated(Num_Acc))) %>% 
  filter(!is.na(Num_Acc)) %>% 
  count(id)

# On met la bonne géométrie à periph_count
st_geometry(periph_count)=troncons.geom[match(periph_count$id,troncons$id)]

# On fait une carte
ggplot(periph_count) + 
  #ggspatial::annotation_map_tile(zoom=13,type="stamenbw") + 
  geom_sf(data=roads.geom, colour = "#666666",size=0.5)+
  geom_sf(aes(color=n),size=3)+
  scale_color_distiller("",palette = "Reds",direction=1)+
  coord_sf(crs = 2154, datum = NA) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "white",color=NA),
        plot.background = element_rect(fill = "white",color=NA)) +
  labs(title = "Nombre de personnes accidentées sur le périphérique",
       subtitle = "en 2019, par portion de 500m",
       caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",x="",y="")

Crédits et reproductibilité

Présentation faite grâce au package rmdformats.

Elle s’inspire, ainsi que son tutoriel, d’une précédente formation donnée par les mêmes auteurs avec Timothée Giraud.

Partage de la configuration de R et des packages utilisés :

sessionInfo()

## R version 4.0.5 (2021-03-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19041)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
## [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
## [5] LC_TIME=French_France.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] sfnetworks_0.5.2   tidygraph_1.2.0    RColorBrewer_1.1-2 tidyr_1.1.3       
##  [5] ggplot2_3.3.3      mapview_2.9.9      tidygeocoder_1.0.3 banR_0.2.2        
##  [9] dplyr_1.0.5        sf_0.9-8          
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152            spatstat.sparse_2.0-0   satellite_1.0.2        
##  [4] webshot_0.5.2           httr_1.4.2              tools_4.0.5            
##  [7] bslib_0.2.4             utf8_1.2.1              R6_2.5.0               
## [10] rpart_4.1-15            KernSmooth_2.23-18      mgcv_1.8-34            
## [13] DBI_1.1.1               colorspace_2.0-0        raster_3.4-10          
## [16] withr_2.4.2             sp_1.4-5                tidyselect_1.1.1       
## [19] leaflet_2.0.4.1         curl_4.3.1              compiler_4.0.5         
## [22] cli_2.5.0               leafem_0.1.6            labeling_0.4.2         
## [25] bookdown_0.22.3         sass_0.3.1              scales_1.1.1           
## [28] spatstat.data_2.1-0     classInt_0.4-3          readr_1.4.0            
## [31] proxy_0.4-25            goftest_1.2-2           spatstat_2.1-0         
## [34] systemfonts_1.0.1       stringr_1.4.0           digest_0.6.27          
## [37] spatstat.utils_2.1-0    rmarkdown_2.9.1         svglite_2.0.0          
## [40] base64enc_0.1-3         pkgconfig_2.0.3         htmltools_0.5.1.1      
## [43] highr_0.9               htmlwidgets_1.5.3       rlang_0.4.11           
## [46] rstudioapi_0.13         jquerylib_0.1.4         farver_2.1.0           
## [49] generics_0.1.0          jsonlite_1.7.2          crosstalk_1.1.1        
## [52] magrittr_2.0.1          spatstat.linnet_2.1-1   Matrix_1.3-2           
## [55] Rcpp_1.0.6              munsell_0.5.0           fansi_0.4.2            
## [58] abind_1.4-5             lifecycle_1.0.0         stringi_1.6.2          
## [61] yaml_2.2.1              grid_4.0.5              crayon_1.4.1           
## [64] deldir_0.2-10           lattice_0.20-41         splines_4.0.5          
## [67] tensor_1.5              hms_1.0.0               leafpop_0.1.0          
## [70] knitr_1.33              pillar_1.6.0            igraph_1.2.6           
## [73] uuid_0.1-4              spatstat.geom_2.1-0     codetools_0.2-18       
## [76] stats4_4.0.5            glue_1.4.2              evaluate_0.14          
## [79] leaflet.providers_1.9.0 png_0.1-7               vctrs_0.3.8            
## [82] rmdformats_1.0.2        gtable_0.3.0            purrr_0.3.4            
## [85] spatstat.core_2.1-2     polyclip_1.10-0         assertthat_0.2.1       
## [88] xfun_0.24               mime_0.10               lwgeom_0.2-6           
## [91] e1071_1.7-6             class_7.3-18            tibble_3.1.1           
## [94] sfheaders_0.4.0         units_0.7-1             ellipsis_0.3.2         
## [97] brew_1.0-6