Interactions potentielles entre les lieux : application du modèle gravitaire aux villes du nord de la France en 1836
Julie Gravier
21 décembre 2017
Licence du post CC BY-NC 2.0 FR : https://creativecommons.org/licenses/by-nc/2.0/fr/
Objectifs
Implémentation du modèle gravitaire pour l’appliquer à des périodes anciennes pour lesquelles les données sur les flux entre les lieux sont soit inexistantes, soit non exhaustives d’un point de vue spatial. Le test est ici réalisé sur les villes du nord de la France avant la Révolution industrielle (en 1836). Le modèle gravitaire généralisé par Pareto s’écrit : \[F_{ij} = k . \frac{P_i*P_j}{D_{ij}^\alpha}\] On a donc le flux entre un lieu i et un lieu j (\(F_{ij}\)) est égal à \(k\) multiplé par le poids de i fois le poids de j (\(P{i}*P{j}\)) divisé par la distance entre i et j (\(D_{ij}\)) à la puissance \(\alpha\). Dans la bibliographie relative au modèle gravitaire, les poids sont également appelés les masses des lieux (notés \(M{i}\) et \(M{j}\)). Le poids est relatif à la taille des lieux, le plus souvent mesurée par la population recensée. Quant à la distance, il s’agit généralement d’une distance euclidienne.
Dans ce modèle, \(k\) et \(\alpha\) sont des paramètres, avec \(k\) correspondant à une constante (de proportionnalité) permettant d’ajuster le modèle selon les types de flux (quand ils sont donc connus) et \(\alpha\) correspondant au frein de la distance. Ainsi, le modèle testé pour évaluer des flux potentiels est le suivant : \[F_{ij} = \frac{P_i*P_j}{D_{ij}^\alpha}\]
De l’intérêt des données
Les données ont été construites à partir des cartes d’État-major disponibles sur le géoportail de l’IGN et des données de recensement de population à échelle communale en 1836 (en ligne : projet Cassini du Ldh-EHESS). Plus précisément, les lieux de peuplement que sont les “villes” ont été repérées sur les cartes. On a ainsi pris en compte les lieux désignés comme tels par les acteurs de l’époque.
L’intérêt de tester le modèle sur ces données est triple :
1) on possède des données sur les populations ;
2) on se situe avant la Révolution industrielle, donc avant l’arrivée du chemin de fer. Ainsi, les modes de déplacements, et dès lors les relations distance-temps, sont semblables aux époques antérieures que l’on étudie par ailleurs (de l’Antiquité à l’époque moderne). Si le modèle apparaît relativement satisfaisant, il pourra ainsi être appliqué sur des données plus anciennes ;
3) sachant que l’on a pour la même époque (1833), des données sur des réseaux existants (à savoir le réseau des routes de Poste grâce à la base de données réalisée par Anne Bretagnolle et Nicolas Verdier), on pourra faire des comparaisons entre les résultats du modèle selon différentes valeurs de paramètre \(\alpha\) et le réseau observé. La comparaison nous permettra de “valider” le modèle, disons tout du moins, d’évaluer des valeurs de paramètres efficaces dans un contexte préindustriel.
Format des données en entrée
En entrée, nous avons d’une part un tableau (en .csv) avec les villes étudiées (i) et les données sur les populations. Et nous avons d’autre part les localisations des villes i en format .shp.
Il aurait été plus simple pour les traitements d’avoir un seul tableau avec les coordonnées géographiques des villes en tant qu’attribut. Toutefois, l’enregistrement n’a pas été réalisé de cette manière au départ (en particulier car les localisations ont été récupérés directement par calcul dans un SIG du centroïde des communes car nous avions anciennement réalisé d’autres traitement en SIG demandant de le faire). Par ailleurs, ces formats de données sont souvent utilisés dans les chaînes d’enregistrement et de traitements de données des archéologues. Ils peuvent ainsi, selon nous, s’approprier et appliquer le code suivant plus facilement.
Traitements
Packages nécessaires
library(sp)
library(rgdal)## rgdal: version: 1.2-7, (SVN revision 660)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/Paris/Documents/R/win-library/3.4/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: C:/Users/Paris/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-4library(Matrix)
library(spdep)
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3library(cartography)
library(igraph)##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionlibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(reshape2)Imports des données
# Le semis étudié
Semis <- readOGR(dsn = "Data_Villes_PasBourgs.shp",
layer = "Data_Villes_PasBourgs",
encoding = "utf8",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "Data_Villes_PasBourgs.shp", layer: "Data_Villes_PasBourgs"
## with 72 features
## It has 30 fields
## Integer64 fields read as strings: ID_GEOFLA SUPERFICIE Tquant_Ran Tquant_POP Tquant_P_1 Tquant_P_2 Tquant_P_3# Un fond géographique pour visualiser les données
Fond <- readOGR(dsn = "G:/These/These/Sujet/3-Traitements/2-SIG/20-SurfacesVilles/1-BDAB&ETJM_These/ModeleGravitaire/PorjetR/Fond_R.shp",
layer = "Fond_R",
encoding = "utf8",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "G:/These/These/Sujet/3-Traitements/2-SIG/20-SurfacesVilles/1-BDAB&ETJM_These/ModeleGravitaire/PorjetR/Fond_R.shp", layer: "Fond_R"
## with 15 features
## It has 1 fields
## Integer64 fields read as strings: id# Les données sur les populations
PopulationDonneesEtmj <- read.csv("Data_R.csv", header = T, sep = ";", dec = ",", stringsAsFactors = F)Vérification des projections géographiques
proj4string(Semis)## [1] "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"proj4string(Fond)## [1] "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs"Qu’a-t-on dans le jeu de données ?
head(PopulationDonneesEtmj)## INSEE NOM Etjm3km2 Etjm3 Rang Sup30ha TypeEtmj POP1793
## 1 80001 ABBEVILLE 3.079000 307.9000 5 307,9 Ville 18125
## 2 80013 AIRAINES 0.467146 46.7146 53 46,7146 Bourg 1490
## 3 62014 AIRE-SUR-LA-LYS 0.482772 48.2772 50 48,2772 Ville 8500
## 4 80016 ALBERT 0.539511 53.9511 44 53,9511 Ville 2000
## 5 80021 AMIENS 3.345310 334.5310 4 334,531 Ville 40000
## 6 2018 ANIZY-LE-CHÂTEAU 0.303141 30.3141 90 30,3141 Bourg 1017
## POP1821 POP1836 POP1866
## 1 18654 18247 19385
## 2 1857 1935 2270
## 3 8713 8717 8803
## 4 2299 2542 4019
## 5 41107 46129 61063
## 6 1043 1073 1116On recode les ID et on sélectionne les villes (hors Paris)
PopulationDonneesEtmj$INSEE <- ifelse(test = nchar(as.character(PopulationDonneesEtmj$INSEE)) == 4,
yes = paste("0", PopulationDonneesEtmj$INSEE, sep=""),
no = PopulationDonneesEtmj$INSEE)
# On vire Paris qui est beaucoup trop exeptionnel
PopulationDonneesEtmj <- PopulationDonneesEtmj %>%
filter(TypeEtmj == "Ville") %>%
filter(NOM != "PARIS")## Warning: package 'bindrcpp' was built under R version 3.4.3Création des deux matrices
Matrice des poids de i fois les poids de j (\(P_{i}*P_{j}\))
On va tout d’abord réordonner les lignes pour être certain d’avoir les mêmes matrices finales (celle des poids et celle des distances)
# Réordonnancement des lignes
PopulationDonneesEtmj <- PopulationDonneesEtmj[order(PopulationDonneesEtmj$INSEE),]
# création de la matrice
MatricePoids <- matrix(nrow = nrow(PopulationDonneesEtmj), ncol = nrow(PopulationDonneesEtmj), byrow = TRUE,
dimnames = list(PopulationDonneesEtmj$INSEE, PopulationDonneesEtmj$INSEE))
head(MatricePoids)## 02173 02304 02361 02408 02468 02691 02722 02789 02810 59017 59036
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 59043 59053 59067 59092 59094 59122 59135 59136 59153 59178 59183
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 59273 59295 59331 59350 59392 59481 59512 59571 59599 59606 60057
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 60104 60143 60157 60159 60175 60176 60471 60509 60612 62014 62038
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 62041 62119 62160 62193 62268 62318 62447 62588 62765 62767 62770
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 76255 77153 80001 80016 80021 80136 80212 80253 80410 80561 80585
## 02173 NA NA NA NA NA NA NA NA NA NA NA
## 02304 NA NA NA NA NA NA NA NA NA NA NA
## 02361 NA NA NA NA NA NA NA NA NA NA NA
## 02408 NA NA NA NA NA NA NA NA NA NA NA
## 02468 NA NA NA NA NA NA NA NA NA NA NA
## 02691 NA NA NA NA NA NA NA NA NA NA NA
## 80620 80685 80721 93066 95052
## 02173 NA NA NA NA NA
## 02304 NA NA NA NA NA
## 02361 NA NA NA NA NA
## 02408 NA NA NA NA NA
## 02468 NA NA NA NA NA
## 02691 NA NA NA NA NA# Boucle permettant de remplir la matrice
for(i in 1:nrow(PopulationDonneesEtmj)){
for (j in 1:nrow(PopulationDonneesEtmj))
if (i != j){
MatricePoids[i,j] <- PopulationDonneesEtmj$POP1836[i]*PopulationDonneesEtmj$POP1836[j];
}
}## Warning in PopulationDonneesEtmj$POP1836[i] * PopulationDonneesEtmj
## $POP1836[j]: NA produit par débordement d'entier par le haut
## Warning in PopulationDonneesEtmj$POP1836[i] * PopulationDonneesEtmj
## $POP1836[j]: NA produit par débordement d'entier par le hauthead(MatricePoids)## 02173 02304 02361 02408 02468 02691 02722
## 02173 NA 11884433 14529403 36895090 7558338 92215310 37764792
## 02304 11884433 NA 8591891 21817730 4469586 54531070 22332024
## 02361 14529403 8591891 NA 26673430 5464326 66667370 27302184
## 02408 36895090 21817730 26673430 NA 13875780 169291100 69329520
## 02468 7558338 4469586 5464326 13875780 NA 34681020 14202864
## 02691 92215310 54531070 66667370 169291100 34681020 NA 173281680
## 02789 02810 59017 59036 59043 59053 59067
## 02173 11525793 12068236 29193296 13583490 44431013 7396950 26754544
## 02304 6815721 7136492 17263312 8032530 26274061 4374150 15821168
## 02361 8332611 8724772 21105392 9820230 32121551 5347650 19342288
## 02408 21159330 22155160 53593760 24936900 81567530 13579500 49116640
## 02468 4334706 4538712 10979232 5108580 16709946 2781900 10062048
## 02691 52885470 55374440 133951840 62327100 203869270 33940500 122761760
## 59092 59094 59122 59135 59136 59153 59178
## 02173 5146484 11328541 80003618 20151085 26965245 23746451 85952559
## 02304 3043348 6699077 47309746 11916245 15945765 14042347 50827623
## 02361 3720668 8190007 57838886 14568295 19494615 17167577 62139693
## 02408 9448040 20797210 146872580 36993850 49503450 43594310 157793790
## 02468 1935528 4260522 30088356 7578570 10141290 8930742 32325678
## 02691 23614360 51980390 367092220 92462150 123728550 108959290 394388610
## 59183 59273 59295 59331 59350 59392 59481
## 02173 106731264 20361786 34402542 16492957 322798415 28525329 14708723
## 02304 63115008 12040842 20343774 9753029 190885255 16868313 8697931
## 02361 77161728 14720622 24871434 11923639 233368205 20622483 10633721
## 02408 195939840 37380660 63157020 30278170 592601150 52367490 27002630
## 02468 40140288 7657812 12938364 6202794 121400430 10728018 5531766
## 02691 489730560 93428940 157854180 75677030 1481142850 130886910 67490170
## 59512 59571 59599 59606 60057 60104 60143
## 02173 87216765 22401551 89507578 87414017 58646606 10826445 4989579
## 02304 51575205 13247047 52929866 51691849 34680382 6402165 2950563
## 02361 63053655 16195277 64709806 63196259 42398762 7827015 3607233
## 02408 160114650 41125310 164320180 160476770 107664860 19875450 9159990
## 02468 32801130 8424942 33662676 32875314 22056252 4071690 1876518
## 02691 400189350 102788290 410700620 401094430 269096740 49676550 22894410
## 60157 60159 60175 60176 60471 60509 60612
## 02173 14502505 39876285 7576270 11575106 26651435 11126806 22486728
## 02304 8575985 23580645 4480190 6844882 15760195 6579782 13297416
## 02361 10484635 28828695 5477290 8368262 19267745 8044162 16256856
## 02408 26624050 73205850 13908700 21249860 48927350 20426860 41281680
## 02468 5454210 14996970 2849340 4353252 10023270 4184652 8456976
## 02691 66543950 182970150 34763300 53111740 122288650 51054740 103179120
## 62014 62038 62041 62119 62160 62193 62268
## 02173 39078311 9638450 105283255 30506815 115356556 48707795 12319284
## 02304 23108767 5699650 62258735 18040055 68215532 28803115 7284948
## 02361 28251797 6968150 76114885 22055005 83397412 35213465 8906268
## 02408 71740910 17694500 193281550 56005150 211774360 89418950 22616040
## 02468 14696862 3624900 39595710 11473230 43384152 18318390 4633128
## 02691 179308690 44225500 483086450 139978850 529307240 223493050 56526360
## 62318 62447 62588 62765 62767 62770 76255
## 02173 8109747 15493248 17335761 85320456 15475316 10140546 16761937
## 02304 4795659 9161856 10251417 50453832 9151252 5996562 9912089
## 02361 5862969 11200896 12532947 61682712 11187932 7331142 12118099
## 02408 14888070 28442880 31825410 156633360 28409960 18616260 30771970
## 02468 3049974 5826816 6519762 32087952 5820072 3813732 6303954
## 02691 37211130 71089920 79544190 391488240 71007640 46529340 76911230
## 77153 80001 80016 80021 80136 80212 80253
## 02173 7809386 81801301 11395786 206796307 6298615 11812705 17537496
## 02304 4618042 48372797 6738842 122287979 3724655 6985385 10370712
## 02361 5645822 59138527 8238622 149504089 4553605 8540035 12678792
## 02408 14336660 150172810 20920660 379641670 11563150 21686050 32195760
## 02468 2937012 30764442 4285812 77773494 2368830 4442610 6595632
## 02691 35832940 375340790 52288940 948873530 28900850 54201950 80469840
## 80410 80561 80585 80620 80685 80721 93066
## 02173 9795355 16990570 7643515 18465477 16452610 14726655 41835356
## 02304 5792435 10047290 4519955 10919469 9729170 8708535 24739132
## 02361 7081585 12283390 5525905 13349679 11894470 10646685 30245012
## 02408 17982550 31191700 14032150 33899370 30204100 27035550 76802360
## 02468 3683910 6389940 2874630 6944634 6187620 5538510 15733752
## 02691 44945450 77960300 35071850 84727830 75491900 67572450 191959240
## 95052
## 02173 8401142
## 02304 4967974
## 02361 6073634
## 02408 15423020
## 02468 3159564
## 02691 38548180Matrice des distances (\(D_{ij}\))
À partir du shapefile on va créer une matrice des distance entre les villes
# récupération des coordonnées (si en entrée on a un polygone, la fonction renvoie les centroïdes)
coordVilles <- as.data.frame(coordinates(Semis))
# On enrichit le tableau
colnames(coordVilles) <- c("CoordX", "CoordY")
coordVilles$INSEE <- Semis@data$INSEE_COM
coordVilles <- coordVilles %>%
filter(INSEE != "75101") # on vire Paris également
coordVilles <- coordVilles[order(coordVilles$INSEE),] # Réordonnancement des lignes pour avoir les mêmes matrices
row.names(coordVilles) <- coordVilles$INSEE
DistVilles <- dist(coordVilles, method = "euclidean") # calcul des distances entre les villes (euclidiennes)
matrice_DistVilles <- as.matrix(DistVilles) # transformation en matrice
matrice_DistVilles[matrice_DistVilles == 0] <- NA
head(matrice_DistVilles)## 02173 02304 02361 02408 02468 02691 02722
## 02173 NA 12049.41 42440.59 29656.90 41851.15 25663.82 27867.23
## 02304 12049.41 NA 31500.37 21220.82 29887.65 21201.23 32112.26
## 02361 42440.59 31500.37 NA 36508.03 19604.07 25550.67 61464.90
## 02408 29656.90 21220.82 36508.03 NA 22304.21 39743.80 30040.34
## 02468 41851.15 29887.65 19604.07 22304.21 NA 36918.89 51719.11
## 02691 25663.82 21201.23 25550.67 39743.80 36918.89 NA 52197.21
## 02789 02810 59017 59036 59043 59053 59067
## 02173 56324.60 43278.60 134374.9 94796.38 140100.1 103064.86 170197.0
## 02304 44275.19 51483.54 132189.1 86191.78 138640.1 95224.01 169539.2
## 02361 22661.19 82487.20 117826.6 65804.85 125651.1 73057.82 157166.9
## 02408 36975.01 52812.24 147199.7 86768.56 154493.0 99527.44 186434.1
## 02468 15280.40 74297.49 135298.8 71537.66 143358.1 83295.75 175522.3
## 02691 46632.65 68918.59 113219.1 79386.56 119452.8 84050.20 149720.9
## 59092 59094 59122 59135 59136 59153 59178
## 02173 93269.70 173737.5 83665.81 151892.9 80861.80 113698.65 102657.18
## 02304 89011.07 174189.8 80372.32 151519.4 74999.96 107149.96 99938.05
## 02361 74539.14 164425.4 70049.44 140561.2 61133.63 85323.92 87277.42
## 02408 100016.30 192053.3 91944.81 168520.7 81578.81 115131.46 113521.69
## 02468 88431.23 182699.3 82829.67 158526.5 70471.14 99076.81 102864.80
## 02691 74283.56 154638.4 66966.13 132160.6 65431.51 92263.99 83177.80
## 59183 59273 59295 59331 59350 59392 59481
## 02173 178584.3 180778.0 144254.3 86687.91 126939.4 107676.67 95219.49
## 02304 178160.9 181281.7 143907.5 79595.28 123860.3 98999.85 88534.72
## 02361 166091.6 171431.4 133394.2 62380.91 107972.1 75523.77 69442.39
## 02408 195314.1 199223.0 160840.7 83944.30 137682.4 100995.28 95042.73
## 02468 184554.4 189798.4 151144.9 70781.15 124955.5 83885.31 80743.89
## 02691 158262.4 161646.4 124782.8 70734.46 105417.1 89388.27 76564.92
## 59512 59571 59599 59606 60057 60104 60143
## 02173 132065.0 86819.74 135916.6 102574.81 102452.4 88589.75 119396.9
## 02304 128047.9 81013.06 131991.4 96641.24 112338.4 97034.76 129977.3
## 02361 109737.3 65364.83 113644.8 77424.56 135256.2 115834.79 155059.1
## 02408 140710.0 88534.05 144833.5 105268.36 126120.2 112349.13 142621.3
## 02468 126539.5 76338.41 130610.6 90852.75 138741.6 121617.44 157182.3
## 02691 109738.1 69589.23 113478.7 82312.05 113114.4 94422.42 132633.1
## 60157 60159 60175 60176 60471 60509 60612
## 02173 86792.72 68320.25 88677.12 75721.13 60442.60 81194.22 87105.28
## 02304 96015.39 75026.24 97857.93 82554.46 64527.08 89721.93 95859.15
## 02361 119951.04 97503.05 123534.28 107284.50 81318.21 114515.50 121841.99
## 02408 107266.22 82373.72 106604.35 86392.99 73071.03 98146.45 102998.93
## 02468 120986.11 95745.82 122208.24 102467.55 81882.93 113106.12 119219.91
## 02691 99893.12 82047.96 104723.83 93221.29 67962.91 96604.30 104235.87
## 62014 62038 62041 62119 62160 62193 62268
## 02173 141645.0 175942.2 101019.60 124671.7 179104.2 186386.1 165241.4
## 02304 142404.6 177551.6 101314.89 124556.6 182756.5 188111.2 168378.5
## 02361 134916.1 170738.0 96273.67 116044.3 181224.1 181277.9 165945.0
## 02408 160046.2 196162.9 117423.93 141183.2 202441.3 206869.1 187707.9
## 02468 152179.8 188614.8 111198.51 132784.6 198175.3 199290.2 182884.2
## 02691 123964.0 158655.7 85416.97 106568.5 165355.6 169112.6 150819.0
## 62318 62447 62588 62765 62767 62770 76255
## 02173 162312.6 133626.2 152610.8 157106.8 121969.6 136389.5 155887.2
## 02304 166864.2 137624.6 157084.6 157968.8 124351.4 136540.6 163460.4
## 02361 168383.9 139049.6 158806.6 150003.5 122370.2 127896.6 174894.2
## 02408 186636.4 156737.2 176684.3 175860.1 142505.6 153589.0 182120.0
## 02468 184304.9 154338.4 174478.5 167660.5 138022.3 145060.0 186451.0
## 02691 150790.1 122008.8 141232.5 139139.9 107935.4 118123.7 153375.5
## 77153 80001 80016 80021 80136 80212 80253
## 02173 105454.8 137725.8 97859.36 107340.8 93554.07 99646.76 116973.3
## 02304 112934.8 143813.7 100737.22 113014.2 96266.53 104105.42 120560.6
## 02361 138076.6 152098.0 105329.72 123454.2 101512.75 112536.12 123883.6
## 02408 116378.7 161928.3 115287.83 128890.3 110001.59 119097.00 137288.1
## 02468 133057.2 164150.4 115327.40 132661.7 110543.07 121604.31 136345.4
## 02691 123347.5 132638.4 91268.61 105183.5 87979.09 96150.30 108218.3
## 80410 80561 80585 80620 80685 80721 93066
## 02173 80319.89 91438.44 82998.20 88123.04 84834.75 152237.0 134431.6
## 02304 81542.43 97327.57 85651.60 89365.69 88945.81 158530.0 142670.3
## 02361 88971.12 112546.00 94854.71 92905.13 101341.40 166507.5 168005.3
## 02408 89839.90 109431.48 95946.86 101034.81 99706.19 177147.9 147316.8
## 02468 92523.06 116930.22 99544.66 100692.18 105431.52 179169.7 163969.2
## 02691 79777.80 95740.66 82848.56 82338.88 87279.60 146728.0 151788.9
## 95052
## 02173 126201.9
## 02304 134368.0
## 02361 157706.1
## 02408 141632.7
## 02468 156200.9
## 02691 140498.6Au cas où on a un problème de correspondance entre les matrices (pas les mêmes villes par exemple), on réalise un test pour vérifier qu’il n’y a pas de soucis dans nos données.
TestIdentifiants <- inner_join(PopulationDonneesEtmj, coordVilles, "INSEE")
nrow(PopulationDonneesEtmj)## [1] 71nrow(TestIdentifiants) # doit être égal au nbre de ligne du tableau en entrée sur les populations ou les coordonnées## [1] 71Calcul des flux potentiels (\(F_{ij}\))
On prend ici une valeur d’\(\alpha\) égale à 2, ce qui est relativement traditionnel (cf. modèle de Stewart). L’objectif étant bien entendu (mais pas ici) de tester différentes valeurs d’\(\alpha\) car c’est un paramètre.
# Voici la différence entre alpha = 1 et alpha = 2
curve(1/x, from=1, to=50, xlab="Distance", ylab="Fij", col = "blue")
curve(1/(x^2), from=1, to=50, xlab="Distance", ylab="Fij", add = T, col = "red")
# Matrice des flux potentiels
FluxPotentielsDijCarre <- MatricePoids/(matrice_DistVilles^2)Transformation de la matrice en tableau à trois entrées :
- l’origine du flux (i),
- la destination du flux (j),
- le flux potentiel
OrigineDestDij <- melt(FluxPotentielsDijCarre, id = c("INSEE_Dep", "INSEE_Arr"), measure = "FluxPotentielDij")
# On va virer les NA (= ce sont les i et j identiques)
OrigineDestDij <- OrigineDestDij %>%
filter(!is.na(value))
# on renomme les colonnes et réorganise le tableau
colnames(OrigineDestDij) <- c("Orig", "Dest", "Fij")
OrigineDestDij <- OrigineDestDij[order(OrigineDestDij$Orig),]
head(OrigineDestDij)## Orig Dest Fij
## 71 2173 2304 0.081855312
## 141 2173 2361 0.008066497
## 211 2173 2408 0.041948563
## 281 2173 2468 0.004315306
## 351 2173 2691 0.140010390
## 421 2173 2722 0.048629453Création d’un graphe des flux potentiels
On crée un graphe à partir du tableau des flux origine-destination. Du fait que les flux sont symétriques, on ne va pas garder la totalité des \(F_{ij}\), mais on ne va sélectionner qu’une partie de la matrice.
Graphe_Dij <- graph_from_data_frame(OrigineDestDij, directed = F)
summary(Graphe_Dij)## IGRAPH 1b10ee5 UN-- 71 4968 --
## + attr: name (v/c), Fij (e/n)# On simplifie le graphe <=> on ne sélectionne qu'une partie (inférieure ou supérieure) de la matrice des flux
Graphe_Dij <- simplify(Graphe_Dij)
summary(Graphe_Dij)## IGRAPH 1b12c31 UN-- 71 2484 --
## + attr: name (v/c)On va ensuite enrichir le graphe de données, sur les liens et sur les noeuds.
# Enrichissement des liens
doublonFij <- !duplicated(OrigineDestDij$Fij)
sum(doublonFij, na.rm=FALSE) # doit être égale au nombre de liens du graphe## [1] 2484OrigineDestDij <- OrigineDestDij[doublonFij,]
E(Graphe_Dij)$Fij <- OrigineDestDij$Fij
# Enrichissement des noeuds
V(Graphe_Dij)$POP <- PopulationDonneesEtmj$POP1836
V(Graphe_Dij)$NOM <- PopulationDonneesEtmj$NOM
summary(Graphe_Dij)## IGRAPH 1b12c31 UN-- 71 2484 --
## + attr: name (v/c), POP (v/n), NOM (v/c), Fij (e/n)Pour représenter le graphe, on a besoin que chaque noeud soit géocodé.
# On récupère les coordonnées des points d'après le tableau précédent (coordVilles)
coordsGraphe <- as.matrix(coordVilles[,1:2])
Graphe_Dij$layout <- coordsGrapheVisualisation du graphe
On peut maintenant visualiser les flux potentiels entre les villes !
plot(Fond, border = "grey")
plot.igraph(Graphe_Dij, rescale = F,
vertex.size = 2,
vertex.label = V(Graphe_Dij)$NOM, vertex.label.cex = 0.7,
edge.curved = 0.4,
edge.color = "red", add=T)
En toute logique, on ne voit pas grand chose car il y a des flux potentiels entre toutes les villes. On va donc sélectionner une partie seulement des flux à représenter. Pour cela, on décide de représenter que les 10% des flux les plus importants (en volume).
# Calcul des déciles
DecilesFlux <- quantile(E(Graphe_Dij)$Fij, probs = seq(0,1,0.1))
DecilesFlux <- as.data.frame(DecilesFlux)
# Sélection d'une partie du graphe
GraphePrincipalDij <- delete_edges(graph = Graphe_Dij,
edges = E(Graphe_Dij)[E(Graphe_Dij)$Fij < DecilesFlux[10,1]])
summary(GraphePrincipalDij)## IGRAPH 1bee323 UN-- 71 249 --
## + attr: layout (g/n), name (v/c), POP (v/n), NOM (v/c), Fij (e/n)Visualisation
plot(Fond, border = "grey")
plot.igraph(GraphePrincipalDij, rescale = F,
vertex.size = 1,
vertex.label = V(GraphePrincipalDij)$NOM, vertex.label.cex = 0.5,
vertex.color = "grey50",
vertex.label.color = "black",
edge.curved = 0.2,
edge.width = E(GraphePrincipalDij)$Fij,
edge.color = "darkorange", add=T)
On peut visualiser les flux différemment, en particulier sous forme de heatmap. On récupère sous forme matricielle le deuxième graphe (on aurait également pu travailler directement sur le tableau origine/destination en ne sélectionnant que les 10% des flux majeurs et en le transformant en matrice), puis on le représente en heatmap.
# Récupération de la matrice d'adjacence du graphe composé des flux majeurs
MatGrapheMaj <- as_adjacency_matrix(GraphePrincipalDij,
attr = "Fij", sparse = F) # les valeurs de la matrice sont celles des flux
row.names(MatGrapheMaj) <- V(GraphePrincipalDij)$NOM
colnames(MatGrapheMaj) <- V(GraphePrincipalDij)$NOM
# Représentation sous forme de heatmap
PalettePerso <- colorRampPalette(c("gold", "red"))
heatmap(MatGrapheMaj, Rowv = NA, Colv = NA,
col = PalettePerso(100),
scale = "none")
On observe que les flux principaux sont entre Lille, Roubaix et Tourcoing. Les flux potentiels entre ces trois villes sont tellement important qu’on ne visualise pas bien tous les autres.
Si l’objectif est de visualiser quelles sont les villes qui sont les plus connectées, on peut également sélectionner la matrice d’adjacence du 2e graphe mais en regardant les flux en présence/absence.
# Récupération de la matrice d'adjacence du graphe composé des flux majeurs
MatGrapheMaj_pres <- as_adjacency_matrix(GraphePrincipalDij, sparse = F)
row.names(MatGrapheMaj_pres) <- V(GraphePrincipalDij)$NOM
colnames(MatGrapheMaj_pres) <- V(GraphePrincipalDij)$NOM
# Heatmap présence/absence de flux majeurs
heatmap(MatGrapheMaj_pres, Rowv = NA, Colv = NA,
col = PalettePerso(100),
scale = "none")
# Heatmap avec réordonnancement des lignes et colonnes
heatmap(MatGrapheMaj_pres,
col = PalettePerso(100),
scale = "none")
On observe ici que les villes les mieux connectées sont Amiens, Lille, Boulogne-sur-Mer, etc. D’autres villes sont moyennement connectées (plus ou moins celles comprises entre Bergues et Beauvais sur le heatmap). D’autres sont très peu connectées, telles que Noyon, Breteuil, Clermont. Enfin, certaines villes n’ont aucune connexion (quand on regarde les 10% des flux les plus importants), telles que Guise ou Ham.