library(tidyverse)
library(sf)
library(spatstat.linnet)

Dataset

Le jeu de données est issu du programme Dynarep, en open access (licence Open Database License) et accessible sur https://www.fabriquenumeriquedupasse.fr/explore/dataset/dynarep-les-voies-dapres-la-carte-detat-major-au-1-320e-1854/.

network_data <- read_sf("dynarep-les-voies-dapres-la-carte-detat-major-au-1-320e-1854.shp") %>%
  st_transform(x = ., crs = 2154) %>%
  st_cast(x = ., to = "LINESTRING")

plot(network_data$geometry)

Import d’un jeu de données sur les villes en cours d’étude:

towns_data <- read_sf("communes_edition_ign_2016.gpkg") %>%
  st_centroid(x = .)

# snaping towns points on road network
snaping_data <- maptools::snapPointsToLines(points = as(object = towns_data, Class = "Spatial"),
                                            lines = as(network_data, "Spatial"),
                                            maxDist = 5000) %>%
  st_as_sf(x = .)

# visu
plot(network_data$geometry)
plot(snaping_data$geometry, cex = 0.75, pch = 24, bg = "darkgreen", col = "darkgreen", add = TRUE)

Identification de problèmes de topologie

# network as linnet object
network <- as.linnet(X = as(object = network_data %>% select(geometry), Class = "Spatial"))
# towns as lpp object
lpptowns <- lpp(X = as.ppp(X = as(object = snaping_data %>% select(geometry), Class = "Spatial")), L = network)
# extracting all vertices of road network as points
xy_roads <- tibble(x = network$vertices$x, y = network$vertices$y) %>%
  st_as_sf(x = ., coords = c("x", "y")) %>%
  st_set_crs(x = ., value = 2154)
# creating vertices as lpp object
xy_roadslpp <- lpp(X = as.ppp(X = as(object = xy_roads, Class = "Spatial")), L = network)

# distance matrix calculation between all towns and all vertices of road network
distances <- crossdist.lpp(X = xy_roadslpp, Y = lpptowns)
distances <- distances %>%
  as_tibble() %>%
  rowid_to_column() %>%
  bind_cols(xy_roads)


# ugly visu of vertices bugs: all vertices which are not link with the road network
plot(network_data$geometry)
plot(distances %>% filter(V1 == 'Inf') %>% select(geometry), cex = 0.25, col = 'blue', add = TRUE)

Export:

st_write(obj = distances %>% filter(V1 == 'Inf') %>% select(geometry), dsn = "dynarep_voies_detections_bugs_potentiels_.gpkg")
## Writing layer `dynarep_voies_detections_bugs_potentiels_' to data source 
##   `dynarep_voies_detections_bugs_potentiels_.gpkg' using driver `GPKG'
## Writing 5366 features with 0 fields and geometry type Point.