Sélection de flux avec R : le package flows

L’analyse et la représentation des flux supposent le plus souvent une sélection de ces derniers. La méthode classique en analyse spatiale des flux dominants (voir ce billet) et les méthodes des flux majeurs (garder x premiers flux par lieu, les x flux les plus install.packages(« flows », dependencies = TRUE)importants en volume, les flux représentant plus de x % d’un volume donné etc.) forment les deux grandes familles de sélection couramment employées.

Le package flows, développé en juin-juillet au sein de l’UMS RIATE avec mes collègues Timothée Giraud (voir le blog rgeomatic) et Marianne Guérois permet de choisir ses critères de sélection, de calculer des indicateurs statistiques et, pour les flux dominants, de visualiser la sélection effectuée. Il est disponible sur le CRAN donc s’installe normalement avec un install.packages(« flows », dependencies = TRUE).

Les fonctions de sélection et d’analyse prennent en entrée une matrice carrée origine – destination. Si les données de départ sont de type origine – destination – intensité, la fonction prepflows permet de passer en format carré.

La fonction statmat fournit pour une matrice origine destination un certain nombre d’éléments utiles pour connaître les données de départ mais également justifier statistiquement la méthode ou le seuil retenu. On obtient notamment le volume total de flux et les indicateurs statistiques de base sur ces flux (quartiles, moyenne, etc.), la densité, le nombre, la taille et la composition des composantes connexes, le degré pondéré et non pondéré des sommets. Des sorties graphiques permettent d’étudier la distribution des degrés (pondérés et non pondérés) et des liens (boîte à moustache et courbe de Lorenz sur l’intensité des flux).

En ce qui concerne la sélection, trois familles de méthode sont implémentées :

  • flux dominants classiques ;
  • critères locaux (pour chaque sommet, garder les flux obéissant à tel critère) ;
  • critères globaux (considérant l’ensemble des flux, garder ceux obéissant à tel critère).

La fonction compmat permet de comparer deux matrices origine-destination de même dimension, notamment pour connaître la perte d’information consécutive à une sélection.

Enfin, la représentation graphique est limitée aux résultats issus de la méthode des flux dominants dans la mesure où d’autres packages existent pour visualiser des graphes valués ou réaliser des cartes de flux (un package de cartographie développé par le déjà cité Timothée Giraud et son collègue Nicolas Lambert devrait être disponible très prochainement). Elle permet d’obtenir des figures de ce type :

exemplemap

Exceptionnellement, je ne fournis pas de script : je me permets de renvoyer à la vignette qui accompagne le package et qui explique tant les méthodes que l’utilisation des fonctions.

Dernier point : n’hésitez pas à envoyer vos remarques (par mail ou en commentaire), demandes et critiques, pour que nous puissions améliorer ce package.

R et igraph – revue de détail

Suite du billet consacré à statnet il y a peu ici même.

L’objectif du billet est le même que pour celui signalé une ligne au dessus : importer et effectuer une série de mesures sur le graphe suivant. Téléchargement des fichiers de données et du script commenté ici (encodage en utf-8).
td1plot

Importer ses données

L’importation des données avec igraph est plus simple dans la mesure où les liens et les attributs des sommets peuvent être importés simultanément à l’aide de la commande graph.data.frame – et la gestion des sommets isolés ne pose donc aucun problème particulier. En ce qui concerne les liens, igraph considère que les deux premières colonnes correspondent aux origines et aux destinations, toutes les éventuelles colonnes suivantes seront considérées comme des attributs des liens.

Après transformation, taper le nom de l’objet créé permet d’obtenir à peu près les mêmes informations qu’avec statnet mais sous une forme un peu moins lisible.

library(igraph)

liens <- read.table("liens2.txt", header=TRUE, sep=";")
attributs <- read.table("attributs.txt", header=TRUE, sep=";")

g <- graph.data.frame(liens, vertices = attributs)
g
IGRAPH DN-- 7 8 -- 
+ attr: name (v/c), ANC (v/n), TEAM (v/c), I (e/n)

g est donc un objet igraph de type directed network (DN), comprenant 7 sommets et 8 liens. Il existe 3 attributs pour les sommets (attr v), deux de type character (name et TEAM), un de type numeric (ANC (v/n)) ; un attribut numérique I pour le lien.

Filtrer les données

Sélectionner une partie du graphe en fonction d’un attribut (des sommets et/ou des liens) est d’une grande simplicité une fois comprise la logique de la syntaxe.

Si je souhaite supprimer les liens d’une intensité inférieure à 2

g2 <- delete.edges(g,E(g)[I<2])

soit, si l’on décompose, la fonction delete.edges (nom du graphe, E comme edge (nom du graphe) [ attribut et condition]). Si l’on veut supprimer des sommets, on remplace par delete.vertices et E devient V comme ici où je supprime les sommets de degré inférieur à 2.

g22 <- delete.vertices(g, V(g)[degree(g) < 2])

Si l’on supprime des sommets, les liens adjacents disparaissent également. Si l’on supprime des liens, par défaut, tous les sommets du graphe restent présents.

Mesurer

La recherche de la densité (graph.density) donne – et c’est heureux – le même résultat que statnet, soit 0.19 (8/42).

La recherche des composantes connexes (clusters) renvoie une liste de 3 objets, l’appartenance des sommets (les sommets abcd sont dans la composante connexe numérotée 1), la taille des composantes connexes et leur nombre.

$membership
[1] 1 1 1 1 2 2 3
$csize
[1] 4 2 1
$no
[1] 3

Contrairement à statnet, igraph considère par défaut que la connexité faible est recherchée. Pour la connexité forte, il faudra utiliser l’option mode = « strong ».

La recherche des cliques (fonction cliques) est ici efficace mais, par défaut, igraph fournit toutes les cliques de taille 1, 2 etc. Il est donc conseillé d’indiquer un min. La direction des liens n’est pas et ne peut pas être prise en compte. Les commandes dyad.census et triad.census donnent exactement les mêmes résultats que statnet mais la présentation est moins claire ici pour cette dernière fonction car les codes MAN ne sont pas rappelés. Le résultat reproduit ci-dessous n’est pas très simple à interpréter.

 triad.census(g) 
 [1] 11 12  8  0  1  0  1  0  1  0  0  0  1  0  0  0

La commande degree fournit par défaut est la somme des degrés entrants et sortants. La normalisation (normalized = TRUE) se fait ici par le degré maximal possible. Si je calcule les degrés entrants normalisés, j’obtiens le résultat suivant :

 degree(g, mode="in", normalized = TRUE) 
    a         b      c     d      e       f         g 
0.1666667 0.3333333 0.50 0.000 0.1666667 0.1666667 0.00

Le degré entrant de c est 3, il y a 7 sommets donc le maximum possible est de 6 (pas de boucle) : c a donc un degré entrant normalisé de 0.5. La fonction degree.distribution donne la proportion de sommets ayant un degré zéro, un, deux etc. Mais si on visualise avec la fonction plot(degree.distribution(g)), l’axe des abscisses est faux (il commence à un au lieu de zéro). Et si vous faites partie des adeptes du scale-free, le fameux plot log log s’obtient aisément (voir le script) mais il n’a évidemment aucun intérêt pour un réseau aussi modeste.

La fonction closeness appliquée à un graphe non connexe donne les résultats aberrants suivants :

      a          b          c          d          e          f          g 
0.03333333 0.03333333 0.02380952 0.04000000 0.02777778 0.02777778 0.02380952

Il n’est par définition pas possible de calculer une centralité de proximité pour un sommet isolé, igraph y parvient… La solution permettant d’obtenir un résultat cohérent est de décomposer le graphe puis de la mesurer sur chacune des composantes considérée individuellement. L’option mode= »all » permet d’ignorer la direction des liens.

g2 <- decompose.graph(g)
plot(g2[[1]])
sapply(g2, closeness, mode="all")
[[1]]
        a         b         c         d 
0.2500000 0.3333333 0.3333333 0.2500000 
[[2]]
e f 
1 1 
[[3]]
  g 
NaN

Le résultat est plus crédible (NaN soit Not a Number pour g effectivement). 1 pour e et f, c’est logique. Mais pourquoi 0.25 pour a et d, 0.33 pour b et c. Petit rappel : pour mesurer la closeness, on divise 1 par la somme des plus courts chemins entre un sommet et tous les autres. Pour a et d, on a donc 1/(1+1+2) ; pour b et c, 1/(1+1+1). À retenir en résumé : faire attention avec la fonction closeness qui peut facilement donner des résultats absurdes…

La fonction betweenness donne les mêmes résultats que statnet mais igraph propose également de mesurer l’intermédiarité des liens, c’est-à-dire le nombre de plus courts chemins empruntant un lien donné. Par définition, cette intermédiarité des liens sera toujours supérieure ou égale à 1 (le lien entre a et b est placé sur le plus court chemin entre ces deux points). Les résultats sont présentés dans l’ordre des liens du fichier de départ. Soit le résultat suivant (les codes des liens ont été rajoutés manuellement a posteriori) :

edge.betweenness(g, directed=TRUE) 
   ef fe  ab  ba  ac db dc bc
[1] 1  1   1  2   1  2   1   1
edge.betweenness(g, directed=FALSE)
    ef  fe  ab        ba         ac      db        dc         bc
[1] 0.5 0.5 0.8333333 0.8333333 1.3333333 1.6666667 1.3333333 1.0

Si on prend en compte l’orientation des liens, tous ont un minimum de 1 et le lien ba atteint 2 : il se situe en effet sur les plus courts chemins ba et dba.  Si on ne la prend pas en compte, les résultats sont difficiles à comprendre. Pourquoi des résultats différents pour les liens db, dc, ab et ba ? Tout simplement parce qu’il y a trop de liens pour l’opération demandée : si je ne tiens pas compte de la direction, mon graphe ne comprend que 6 liens.  Donner des résultats pour 7 oblige à des pondérations étranges… La solution la plus pertinente est de transformer d’abord le réseau en réseau non orienté (as.undirected) puis de mesurer l’intermédiarité des liens. On obtient alors un résultat cohérent :

[1] 1.5 1.5 1.0 1.5 1.5 1.0

1 correspond ici aux liens ef et bc et les quatre autres liens obtiennent un score identique. Chacun est placé sur le plus court chemin entre ses sommets terminaux et est l’un des deux plus courts chemins entre a et d : tous ont donc 1.5.

La fonction transitivity (nombre de triangles fermés parmi les voisins d’un sommet) ne permet pas de prendre en compte l’orientation des liens : il est donc nécessaire pour obtenir des résultats corrects de prendre le graphe non orienté. On obtient alors les résultats suivants

transitivity(g3, "local")
[1] 1.0000000 0.6666667 0.6666667 1.0000000       NaN       NaN       NaN

a et d sont situés à la pointe dans des triangles fermés donc la transitivité est de 1 ; efg ne font pas partie de triades donc le résultat n’est pas un nombre (NaN = Not a Numer). b et c sont voisins et tous deux liés à a et à d ; comme il manque 1 lien entre a et d, ils obtiennent un score de 2/3 (2 triangles fermés sont présents sur 3 possibles). La recherche des plus chemins, par défaut, ne tient pas compte non plus de l’orientation des liens. La fonction centralization donne des résultats corrects pour les graphes k-réguliers (score de 0) et non valides pour les graphes en étoile (voir le script où un graphe étoilé de 10 sommets obtient un score de centralisation par le degré de 0.8 seulement). En ce qui concerne les graphes bipartis, même remarque que pour statnet : igraph permet d’en définir, d’en représenter mais ne permet pas de les analyser – lorsqu’une mesure est lancée, il est considéré comme un one-mode.

Pour résumer : il est beaucoup plus facile d’importer ses données et de les filtrer avec igraph, il est aussi beaucoup plus facile d’obtenir des résultats aberrants si l’on n’est pas attentif… À venir très prochainement, un billet détaillé consacré au package bipartite.

R et statnet – revue de détail

Suite aux cours et aux TD donnés à l’ENSAI, j’ai découvert de nouveaux (et agaçants) problèmes liés à l’analyse de réseau avec R. De plus, la nécessité de commenter les résultats m’est clairement apparue. Ce billet vise donc à détailler au maximum les étapes et les résultats produits par le package statnet – un billet consacré à igraph est en cours de préparation. Le script et les fichiers de données sont téléchargeables ici.

Rappel : statnet (et igraph) permettent d’analyser des graphes non valués et one-mode. Si statnet permet de mener des modèles type ERGM, igraph est plus simple d’utilisation. Pour l’analyse de réseau biparti avec R, oubliez statnet et igraph (qui proposent pourtant des options bipartite=TRUE) et voyez plutôt du côté du package bipartite qui est, à ma connaissance, le seul à donner des résultats corrects pour ce type de réseau (voir la démonstration dans le script).

1. Liens, attributs et isolés : sur le format et l’importation des données

Le graphe fictif qui va servir d’exemple tout au long de ce billet représente les interactions autour d’une machine à café dans une petite entreprise : la population est constituée de 7 individus répartis en 2 équipes et de niveaux d’ancienneté différents. Un lien signale qu’une personne a payé un café à une autre. Le graphe est simple : pas de boucle (on ne prend pas en compte les cafés que s’offre une personne) ; pas de liens de nature différente. Il est orienté (une personne offre à une autre) et valué (combien de cafés ont été offerts par une personne à une autre dans la semaine ?).td1plot

Les sommets sont définis par deux attributs : un attribut quantitatif relatif à l’ancienneté (ANC – taille du sommet), un attribut qualitatif relatif à l’équipe (TEAM – couleur du sommet). Les liens sont définis par une direction et une intensité (INT). La valeur des attributs est donnée dans le tableau ci-dessous.

td11

Un rappel important : statnet (comme igraph) ne prend pas en compte de l’intensité des liens pour les calculs. Elle peut juste être utilisée pour la représentation graphique.

Premier problème : statnet génère un graphe à partir d’une liste de liens. Or, le sommet g est isolé. Une solution est de créer dans son fichier de départ un lien fictif g-g puis de préciser lors de l’importation qu’il n’y a pas de boucle.

L’importation des données se fait en deux temps (liens, attributs des sommets). La transformation en objet network, manipulable par statnet, nécessite de ne prendre en compte que les deux premières colonnes (origine – destination). Taper ensuite le nom du graphe permet d’en connaître les caractéristiques. De façon étrange, alors qu’il n’y a pas de boucle, statnet indique malgré tout 9 liens et non 8. Pourtant, le calcul de la densité (cf infra) donne un résultat correct – seuls les 8 liens réellement présents sont donc pris en compte dans les calculs.

library(statnet)
#import fichier
 liens <- read.table("liens.txt", header=TRUE, sep=";")
 attributs <- read.table("attributs.txt", header=TRUE, sep=";")
g <- as.network(as.matrix(liens[,-3]), loops = FALSE)
g
 Network attributes:
 vertices = 7
 directed = TRUE
 hyper = FALSE
 loops = FALSE
 multiple = FALSE
 bipartite = FALSE
 total edges= 9
 missing edges= 0
 non-missing edges= 9
Vertex attribute names:
 vertex.names
No edge attributes

Pour l’instant, les sommets ont un attribut (leur nom) et les liens n’en ont aucun. Les lignes suivantes permettent de transférer les données des data.frame vers l’objet network.

g %v% "ANC" <- attributs$ANC
g %v% "TEAM" <- as.numeric(attributs$TEAM)
g %e% "INT" <- liens$I

Si l’on tape à nouveau g, on observe que les attributs sont désormais plus nombreux.

 Vertex attribute names:
 ANC TEAM vertex.names
Edge attribute names:
 INT

Il est sans aucun doute possible de filtrer le graphe selon les attributs des liens et des sommets (par exemple, extraire le graphe des relations entre membres de l’équipe A) mais si l’opération est facile avec igraph, elle ne l’est pas avec statnet (i.e. je n’ai pas trouvé comment faire, si ce n’est en filtrant les données en amont, avant la transformation en graphe…).

Pour obtenir la figure ci-dessus, deux options sont possibles : appeler les données des data.frame, appeler les attributs de l’objet network. Le résultat est heureusement le même.

#à partir des fichiers de départ
gplot(g, gmode="digraph", displaylabels=TRUE, vertex.col=nodeColors, 
      vertex.cex= log(attributs$ANC), edge.lwd = liens$I*2)

#à partir des attributs de l'objet network
gplot(g, gmode="digraph", displaylabels=TRUE, vertex.col=nodeColor, 
      vertex.cex= log(g %v% "ANC"), edge.lwd = g %e% "INT")

2. Quelques mesures et quelques précisions

Le calcul de la densité (fonction gden) donne comme résultat 0.1904762, ce qui correspond bien à 8 (nombre de liens présents) divisé par 7*6 (nombre de liens possibles) – pour mémoire, dans un graphe simple orienté de n sommets, le nombre maximum de liens est égal à n*(n-1).

Le nombre de composantes connexes s’obtient avec la fonction components et donne ici comme résultat 5. Pourquoi 5 ? Par défaut, statnet considère que, dans un graphe orienté, la connexité recherchée est la forte : il doit exister un chemin (une suite de liens orientés) entre toute paire de sommets au sein d’une composante connexe. La composition des 5 composantes connexes est obtenue avec la commande component.dist. Si l’on souhaite les composantes faiblement connexes (au nombre de 3), il suffit de modifier l’option.

components(g)  #   5 composantes fortement connexes
cs <- component.dist(g) 
cs$membership #  et leur composition

components(g, connected="weak") #  3 composantes faiblement connexes
cw <- component.dist(g,connected="weak")
cw$membership #composition
cw$csize      #et taille

La liste des plus chemins – retour de la commande geodist – comprend deux matrices : $counts donne le nombre de plus courts chemins entre paires de sommets (il peut très bien y avoir plusieurs plus courts chemins entre deux sommets) ; $gdist donne la longueur en liens des plus courts chemins. La commande dyad.census donne un résultat facilement interprétable (nombre de liens mutuels, asymétriques et nuls). La commande triad.census donne en sortie le code MAN (mutuel, asymétrique, nul) des triades et le nombre de triades respectives. Dans le graphe fictif étudié, il n’y a aucune triade complète, et donc 0 sous le code 300.

Il est déconseillé d’utiliser la commande clique.census : même avec un graphe aussi petit, et sur un ordinateur récent et plutôt puissant, RStudio plante…

En ce qui concerne le degré, la commande degree fournit par défaut la somme des degrés entrants et sortants par sommet. Si l’on choisit de normaliser, la somme des degrés des sommets sera égale à 1 (le mode de normalisation est différent dans igraph). Le graphe étant non connexe, la commande closeness donne un résultat égal à zéro pour tous les sommets (igraph donne lui des résultats aberrants). La centralité d’intermédiarité (betweenness) donne ici 0 pour tous les sommets sauf pour b qui obtient 1 : il y a en effet un seul plus court chemin (dba) passant par b. Si l’on ne prend pas en compte la direction des liens, tous les autres sommets restent à zéro et b descend à 0.5. Pourquoi ? Il y a deux plus courts chemins entre d et a (dba et dca), et 1 sur 2 passe par b. Que c reste à zéro est par contre difficilement explicable… (si vous avez une idée, n’hésitez pas à la soumettre).

Enfin, et je termine ici pour l’examen détaillé de certaines fonctions de statnet, la mesure de la centralisation du graphe est correcte dans statnet – on obtient bien 0 avec un graphe type cercle et 1 avec un graphe étoilé – mais elle donne des résultats erronés dans igraph

Dernier point qui concerne les graphes bipartis : statnet permet (difficilement) de les importer puis de calculer des indicateurs mais ceux-ci sont faux. Dans l’exemple du script, nous avons 7 liens, 3 événements et 4 acteurs. Si les 4 acteurs participent aux 3 événements, le graphe est complet, la densité de 1 et le nombre de liens égal à 12. Or, la commande gden donne comme résultat 0.33, soit 7*2 / 7*6 – dans un graphe simple non orienté de n sommets, la densité est égale à 2n / n(n-1). Et si vous tentez l’expérience avec igraph, le résultat est le même. Le seul – et très relatif – intérêt est de pouvoir passer ensuite à des matrices acteurs-acteurs ou événements-événements. Cette dernière opération est possible sans ces packages en multipliant la matrice m par sa transposée pour obtenir la matrice acteurs-acteurs

m %*% t(m)

ou en faisant l’opération inverse (matrice événements-événements)

t(m) %*% m

Les indices Kansky avec R : premières fonctions

K. Kansky est un géographe américain qui, dans sa thèse imprimée en 1963, proposa un certain nombre d’indicateurs permettant de caractériser des réseaux, indicateurs directement issus de la théorie des graphes. Si l’ouvrage est quasiment introuvable, le chapitre 2 Measures of network structure en 1989 paru dans la revue Flux est lui accessible sur Persée (http://www.persee.fr/web/revues/home/prescript/article/flux_1154-2721_1989_hos_5_1_913). Ces indices ont connu un succès foudroyant en géographie des transports et restent encore largement utilisés aujourd’hui. Bizarrement, aucun package R ne permet – à ma connaissance – de calculer ces indices. Le présent billet rappelle très brièvement la définition de ces indices, leur interprétation et fournit les fonctions R – basées sur le package igraph – permettant de les calculer.

Petite précision : il s’agit de ma toute première tentative d’écriture de fonction, merci  donc d’être indulgent-e, ma marge de progression dans ce domaine étant immense… Merci à Serge L’homme qui m’a donné la formule pour le diamètre kilométrique et corrigé mon indice de Shimbel, merci enfin à Timothée Giraud qui est toujours disponible pour répondre à mes nombreuses (et souvent naïves) questions. J’en profite pour signaler leurs sites respectifs : site pédagogique de Serge Lhomme (http://serge.lhomme.pagesperso-orange.fr/) ; blog R géomatique de Timothée Giraud (http://rgeomatic.hypotheses.org/).

Les indices de Kansky

Soit e le nombre de liens, v le nombre de sommets et p le nombre de composantes connexes d’un graphe g, Kansky propose notamment les indices suivants :

  • nombre cyclomatique mu = e – v + p. Plus il est élevé, plus le réseau est dense ;
  • alpha : ratio entre le nombre de circuits présents et le nombre de circuits possibles : il varie entre 0 (arbre) et 1 (graphe complet) ;
  • beta : e/v. Une valeur inférieure à 1 signale un arbre, beta ne peut dépasser 3 pour un graphe planaire ;
  • gamma : ratio entre le nombre de liens présents et le nombre de liens possibles. Une formule est proposée pour les graphes planaires et une autre pour les graphes non planaires. Le terme densité employé dans les autres disciplines est synonyme ;
  • degré de connectivité : v*(v-1) * 1/2(e). Ratio entre le nombre de sommets et le nombre de liens (son interprétation m’échappe – mais je suis loin d’être un spécialiste des graphes planaires et des infrastructures de transport…) ;
  • indice de Shimbel : somme des toutes les distances entre les sommets. Peut être mesuré de façon topologique (distance mesurée en nombre de liens) ou pondéré (par une distance, un flux, une intensité etc.)
  • eta : distance moyenne des liens, exprimé généralement en kilomètres ;
  • pi : flux total divisé par le long des plus courts chemins kilométriques ;
  • theta : flux total divisé par le nombre de liens (intensité moyenne)
  • iota : distance totale divisée par le flux total

Trois fonctions R

Les données doivent avoir la forme suivante : un data.frame (tableau de données) ayant au minimum 3 colonnes : origine, destination et distance. Des colonnes peuvent être ajoutées pour porter des poids divers (flux de marchandises, débit horaire etc.). Cette liste de liens doit être transformée en objet igraph : le package igraph doit donc être appelé avant leur exécution (library(igraph)).

La fonction kansky_un mesure les indicateurs ne prenant pas en compte les poids et les distances : diamètre (topologique), nombre cyclomatique, indice alpha (planaire et non planaire), beta, gamma (planaire et non planaire), le degré de connectivité et l’indice de Shimbel (topologique).

La fonction kansky_di mesure les indicateurs prenant en compte une distance entre les sommets : diamètre (kilométrique), eta, pi et indice de Shimbel (pondéré). Le tableau en entrée doit impérativement avoir trois colonnes (origine, destination, distance).

La fonction kansky_df mesure les indicateurs prenant en compte une distance et une intensité entre les sommets : indices theta et iota. Si le tableau de départ ne comprend pas 4 colonnes minimum, un message  d’erreur apparaît.

Exemple reprenant une des figures utilisées par Kansky pour détailler le nombre cyclomatique (mu)

library(igraph)
d <- rbind(c(1,2), c(2,3), c(2,5), c(2,7), c(3,4),
             c(4,5), c(5,6),  c(6,7)) 
g <- graph.data.frame(d, directed = FALSE)
kansky_un(g)

Fonctions téléchargeables ici.

Analyse de graphe et QGIS 2.0

Vous l’aurez sans doute déjà noté, la version 2.0 de QGIS est sortie en septembre 2013. Cette nouvelle version a suscité une certaine effervescence dans la sphère du SIG Open Source. Pour les développeurs du dimanche, cette évolution induisait quelques soirées de programmation afin d’apporter les modifications nécessaires au fonctionnement d’extensions développées dans des versions antérieures. Ainsi, je me suis attelé à mettre à jour les extensions présentées dans ce blog (billet 1, billet 2, billet 3). Il était selon moi dommage de laisser tomber en désuétude des extensions conçues il y a seulement un an. J’avais de surcroît plein de projets d’amélioration pour ces extensions, mais je me suis finalement laissé distraire (création d’un plugin d’analyse spatiale et d’un plugin permettant de calculer des modèles de localisation-allocation). Comble du ridicule, je n’ai pas pris le temps de publier cette mise à jour. QGIS est aujourd’hui à la version 2.2.

Je ne vais pas détailler dans ce billet les changements induits par cette nouvelle version en matière de programmation (vous pouvez trouver certains détails ici). Quoi qu’il en soit, je mets à votre disposition une nouvelle extension compatible avec QGIS 2.x (cliquez ici pour la télécharger). Une seule nouveauté est à noter : cette extension regroupe les trois extensions développées précédemment. De plus, cette extension incorpore directement NetworkX (bibliothèque sur laquelle ce plugin repose en partie). L’installation de cette extension est donc très simple.

qgisv2

Ce nouveau plugin regroupe le plugin d’analyse de réseau, celui d’analyse de vulnérabilité et permet de définir des communautés.

Pour installer cette extension, il faut tout d’abord dézipper le dossier téléchargé. Ensuite, il faut copier le dossier « network_analysis » dans « C:\Program Files\QGIS xxxx\apps\qgis\python\plugins\ ». Lancez QGIS et le tour est joué !

Une seule chose est à vérifier (car on apprend des commentaires passés). Lorsque vous dézippez un dossier sous Windows, il vous est proposé de l’envoyer vers un dossier du même nom. Ainsi, votre dossier « network_analysis » est lui-même dans un dossier « network_analysis ». Il faut éviter cette redondance. Ainsi, après avoir copié le dossier, vérifiez bien que votre arborescence est la suivante «C:\Program Files\QGIS xxxx\apps\qgis\python\plugins\network_analysis» et pas «C:\Program Files\QGIS xxxx\apps\qgis\python\plugins\network_analysis\network_analysis».

Arborescence correcte.

Arborescence correcte.

Pour effectuer les calculs, il est nécessaire de disposer de deux couches d’information géographique : une première couche contenant les nœuds et une deuxième couche contenant les arcs. La couche des nœuds doit impérativement contenir les identifiants des nœuds. Il est conseillé de choisir comme identifiant des valeurs numériques allant de 1 à n (vous pouvez télécharger un exemple ici). De même, la couche des arcs doit posséder un champ « identifiant ». De surcroît, cette couche doit impérativement comporter deux champs qui renvoient aux identifiants des nœuds formant les extrémités des arcs (vous pouvez télécharger un exemple ici).

A bientôt pour de nouvelles extensions…

Serge Lhomme

Maitre de conférences à l’université Paris Est Créteil – Docteur en géographie de l’université Paris VII spécialité « dynamique des milieux et risques »

More Posts - Website

Communautés, QGIS, R (igraph) et Python…

Toujours dans la continuité des billets précédents (billet 1) (billet 2) qui présentent des plugins d’analyse de réseaux spatiaux, ce billet présente un plugin permettant d’identifier des communautés au sein de ces réseaux. En effet, ce billet permet d’aborder d’une part un sujet important de l’analyse de réseaux, l’identification de communautés, et d’autre part les relations possibles entre R et QGIS via Python, puisque la bibliothèque choisie pour ces analyses est igraph.

igraph, qui est un package réputé sous R pour l’analyse de réseaux, dispose d’une extension Python. Dans un premier temps, il s’agit d’utiliser cette extension qui permet de travailler sous QGIS avec igraph et ce indépendamment de R. Pour utiliser igraph avec Python et par conséquent QGIS, c’est néanmoins un peu plus compliqué qu’avec NetworkX, car il faut cette fois-ci compiler la bibliothèque en question. Le plus simple consiste à télécharger l’installateur Windows pour Python 2.7 (http://igraph.sourceforge.net/download.html). Une fois l’installateur lancé, il suffit de bien préciser l’emplacement de Python dans QGIS « C:/Program Files/Quantum GIS Lisboa/apps/Python27 ». Normalement, dès le premier lancement de QGIS, tout fonctionne. Il est alors possible d’importer igraph via la console Python de QGIS.

Exemple de test du fonctionnement de IGraph sous QGIS via la console Python.

Exemple de test du fonctionnement de IGraph sous QGIS via la console Python.

Il convient ensuite de télécharger le nouveau plugin (cliquer ici pour l’obtenir). Il faut le dézipper, puis placer le dossier dans le répertoire suivant : « C:\Program Files\Quantum GIS Lisboa\apps\qgis\python\plugins ». Il suffit alors de lancer QGIS pour que le plugin apparaisse dans « Extension -> Analyse de Graphe… -> Detection de communautes ». S’il n’apparaît pas, il faut aller dans « Extension -> Gestionnaire d’extensions » et cocher l’extension « Detection de communautes (version 0.1) ». Après avoir lancé le plugin, une boîte de dialogue s’ouvre. Elle permet de sélectionner les nœuds et les arcs du graphe à étudier. Après avoir validé cette première boîte de dialogue, une deuxième boîte de dialogue apparaît. Elle permet de sélectionner les champs des extrémités des arcs, de préciser si le graphe est orienté ou non, de préciser si le graphe est pondéré ou non (et si oui, quel est le champ de cette pondération : distance, flux, importance…). Enfin, différents algorithmes sont proposés afin d’identifier les communautés au sein du réseau étudié.

Présentation de la deuxième boite de dialogue.

Présentation de la deuxième boite de dialogue.

Les différents algorithmes proposés par le plugin utilisant igraph sont les suivants :

  • Walktrap, algorithme développé par Latapy et Pons et fondé sur des marches aléatoires ;
  • Betweenness, algorithme fondé sur la centralité des arcs. Les arcs les plus centraux sont retirés progressivement, chaque nouvelle composante connexe définissant alors une nouvelle communauté ;
  • Multilevel, algorithme fonctionnant par agrégation de nœuds à une communauté afin d’améliorer la modularité ;
  • Eigenvector, algorithme séparant progressivement les communautés à partir des valeurs propres du graphe ;
  • Infomap, deuxième algorithme fondé sur des marches aléatoires.

Pour chaque algorithme, seule la solution optimale est retenue (c’est-à-dire celle offrant la modularité la plus forte). C’est pourquoi le nombre de communautés n’est pas demandé avant le calcul. La valeur finale de la modularité est indiquée afin de comparer les résultats obtenus par les différents algorithmes. Une fois le calcul terminé, les valeurs stockées correspondent aux communautés auxquelles appartiennent les nœuds. Il est alors possible d’afficher ces communautés :

Communautés issues de l’algorithme fondé sur la centralité des arcs.

Communautés issues de l’algorithme fondé sur la centralité des arcs.

Le fait d’utiliser l’extension Python de igraph interrogera sans doute quelques utilisateurs invétérés de R. Faut-il recoder en Python tous ses codes R et trouver une extension Python pour éventuellement les utiliser sous QGIS ? La réponse est non, puisqu’il est possible de faire interagir QGIS et R à l’aide notamment de la bibliothèque RPy. En effet, RPy est une interface Python pour R. Ainsi, cette interface permet en utilisant le langage Python d’interagir avec R et de bénéficier de l’ensemble des développements réalisés sous R. L’inverse est aussi possible. Pour cela, il faut repartir pour une opération de téléchargement (http://sourceforge.net/projects/rpy/files/OldFiles/) et de compilation. En partant du principe que le package igraph est bien installé sous R (il faut notamment vérifier sa présence dans « C:\Program Files\R\R-(x.xx.x)\library\ »), il est très simple d’exécuter des codes R sous QGIS et par conséquent de les exploiter. Par exemple, on va exécuter à partir de QGIS le code R suivant :

test <- function () {
library(igraph)
g <- barabasi.game(100,m=2)
eb <- edge.betweenness.community(g)
return(eb$bridges)
}

Ce code peut être sauvegardé comme étant le fichier test.r. Pour appeler ce fichier sous QGIS, il est possible d’utiliser la console Python de QGIS. Plus précisément, il suffit de taper les lignes de code suivantes :

>>> from rpy import * # Import de la bibliothèque Rpy
>>> r.source('C:/Users/serge.lhomme/Documents/R/test.r') # Chargement du fichier test.r
{'visible': False, 'value': <Robj object at 0x0CBBDA00>}
>>> resultat=r.test() # Exécution de la fonction test
Message d'avis :
package 'igraph' was built under R version 2.12.2
>>> resultat # Affichage du résultat
[198.0, 196.0, 194.0, 192.0, 190.0, 188.0, 186.0, 184.0, 182.0, 180.0, 178.0, 176.0, 174.0, 
172.0, 170.0, 169.0, 168.0, 167.0, 166.0, 165.0, 164.0, 163.0, 162.0, 161.0, 160.0, 159.0, 
158.0, 157.0, 156.0, 155.0, 154.0, 153.0, 151.0, 150.0, 149.0, 148.0, 147.0, 145.0, 144.0, 
143.0, 140.0, 138.0, 137.0, 136.0, 135.0, 134.0, 133.0, 132.0, 131.0, 130.0, 128.0, 126.0, 
124.0, 122.0, 121.0, 119.0, 118.0, 117.0, 115.0, 114.0, 112.0, 111.0, 110.0, 109.0, 108.0, 
106.0, 105.0, 104.0, 102.0, 101.0, 100.0, 98.0, 96.0, 94.0, 92.0, 90.0, 89.0, 87.0, 85.0, 
83.0, 81.0, 79.0, 77.0, 75.0, 73.0, 71.0, 69.0, 67.0, 65.0, 64.0, 62.0, 61.0, 59.0, 57.0, 
56.0, 52.0, 51.0, 47.0, 45.0]

A partir de ces simples commandes, il est possible d’appliquer ses codes R sous QGIS. Ce qui sera sans doute utile pour certains et ce qui augmente encore les possibilités offertes par QGIS.

Serge Lhomme

Maitre de conférences à l’université Paris Est Créteil – Docteur en géographie de l’université Paris VII spécialité « dynamique des milieux et risques »

More Posts - Website

Un plugin QGIS pour analyser la vulnérabilité structurelle de graphes spatiaux

Dans la continuité d’un précédent billet, qui présentait de manière succincte l’utilisation de Python afin de programmer un plugin QGIS d’analyse de réseaux spatiaux, ce billet présente un plugin permettant plus précisément d’analyser leur vulnérabilité structurelle. De plus, ce plugin permet d’introduire des recherches récentes en matière d’analyse de réseaux.

L’analyse de la vulnérabilité structurelle des réseaux connait un intérêt croissant. Les raisons de cet intérêt sont multiples. Ainsi, alors que le risque s’est imposé comme une problématique majeure au sein de nombreuses sociétés contemporaines, l’analyse structurelle constitue un moyen pertinent d’appréhender les problématiques de vulnérabilité. L’analyse structurelle se focalise sur l’étude de la forme, de la configuration des réseaux. La dynamique des flux au sein de ces réseaux n’est donc pas modélisée à l’aide d’équations physiques et les aspects purement matériels (fragilité intrinsèque des nœuds et des arcs) sont laissés de côté. L’approche structurelle permet alors une approche plus réflexive. De plus, elle peut constituer une alternative aux approches physiques modélisant les flux, ces approches physiques étant bien souvent trop complexes à mettre en œuvre sur de grands réseaux. Pour une bonne introduction, il est possible de consulter l’article de Jean-François Gleyze et Magali Reghezza (cliquer ici) ou cette synthèse du groupe fmr (cliquer ici).

Pour les personnes n’étant pas particulièrement intéressées par les problématiques de risques, il est important de préciser que l’utilisation des indicateurs et des algorithmes conçus à cet effet peut se révéler intéressante pour des analyses de réseaux plus classiques, afin notamment de mettre en évidence l’importance de certains composants (nœuds ou arcs).

Pour installer ce plugin, il convient tout d’abord de suivre les recommandations du précédent billet. Ainsi, la bibliothèque NetworkX (cliquer ici pour l’obtenir) doit être installée dans le répertoire des extensions Python de QGIS (par exemple : « C:\Program Files\Quantum GIS Lisboa\apps\Python27\Lib\site-packages »). En effet, pour assurer une certaine continuité avec le précédent plugin, j’ai choisi d’utiliser la même bibliothèque. Ainsi, si cette bibliothèque a déjà été installée pas besoin de recommencer cette installation. Ensuite, il convient de télécharger le nouveau plugin (cliquer ici). Il faut le dézipper, puis placer le dossier dans le répertoire suivant : « C:\Program Files\Quantum GIS Lisboa\apps\qgis\python\plugins ».

Il suffit de lancer QGIS pour que le plugin apparaisse dans « Extension -> Analyse de vulnérabilité ». S’il n’apparaît pas, il faut aller dans « Extension -> Gestionnaire d’extensions » et cocher l’extension « Analyse de la vulnérabilité d’un graphe (version 0.1) ». Après avoir lancé le plugin, une boîte de dialogue s’ouvre. Elle permet de sélectionner les nœuds et les arcs du graphe à étudier. Après avoir validé cette première boîte de dialogue, une deuxième boîte de dialogue apparaît. Elle permet de sélectionner les champs des extrémités des arcs, de préciser si le graphe est orienté ou non, de préciser si le graphe est pondéré ou non (et si oui, quel est le champ de cette pondération : distance, flux, importance…). Enfin, différents indicateurs sont proposés pour analyser la vulnérabilité structurelle du réseau étudié.

Boîte de dialogue et indicateurs

Concernant les indicateurs, l’allongement moyen calcule la moyenne des détours engendrés par la suppression d’un sommet. L’algorithme calcule alors cet indicateur pour tous les sommets. Ainsi, plus l’allongement moyen est important, plus le nœud doit théoriquement faire l’objet d’une attention particulière puisque sa défaillance engendre une reconfiguration importante du réseau. A noter que cet indicateur est corrélé avec l’indicateur de centralité intermédiaire (Gleyze, 2005 ; Lhomme, 2012). Cet indicateur présente le défaut de ne pas prendre en considération les éventuelles pertes de connexité liées à la suppression d’un nœud. En effet, certains détours seraient d’une valeur infinie. Pour pallier ce défaut, il convient alors d’analyser l’inverse de l’éloignement moyen (nommé par certains efficiency). Pour cet indicateur, le terme de centralité d’information a pu être utilisé (Crucitti et al., 2006). A partir de ces deux indicateurs, des cartes thématiques peuvent être réalisées.

Analyse de la centralité d'information

L’indicateur de vulnérabilité dynamique tente d’appréhender les problématiques liées à la reconfiguration du réseau après la défaillance d’un sommet. En effet, la reconfiguration imposée par la suppression d’un sommet peut engendrer une utilisation anormale (notamment une utilisation très supérieure par rapport à la normale) de certains sommets pouvant alors engendrer la défaillance de ceux-ci. Ainsi, l’objectif de cet indicateur est d’étudier les défaillances en chaîne pouvant être induites par la suppression initiale d’un sommet. Pour cela, il existe plusieurs solutions (Wang et Rong, 2011). Ici, l’indicateur calcule la tolérance minimale requise afin d’éviter que la moitié des composants du réseau soient perturbés par la suppression initiale d’un nœud.

Un dernier algorithme permet d’évaluer les conséquences induites par la suppression de plusieurs nœuds. Pour cela, des scénarios types sont étudiés : l’attaque de nœuds choisis aléatoirement, l’attaque des nœuds les plus centraux, l’attaque des nœuds de degrés les plus élevés, l’attaque des nœuds les plus centraux dynamiquement. L’ensemble de ces scénarios est présenté dans (Albert et al., 2004). En fonction du nombre de sommets perturbés, il convient alors d’évaluer les pertes de connectivité correspondantes. Comparativement aux autres indicateurs, une boîte de dialogue supplémentaire permet à l’utilisateur de préciser : le nombre d‘itérations à effectuer pour calculer le scénario aléatoire, le taux maximum de nœuds impliqués dans les scénarios (0 à 1), le répertoire où sera enregistrée l’image, le nom du fichier image, le type de format de l’image (.png ou .svg), la largeur et la hauteur de l’image, l’épaisseur des lignes et la taille des points.

Boîte de dialogue particulière à l'analyse de la vulnérabilité combinatoire

Ensuite, une boîte de dialogue affiche la progression du calcul. En effet, ce calcul, comme ceux des indicateurs d’allongement moyen et de vulnérabilité dynamique, peut être assez long. Enfin, le graphique final s’affiche.

Analyse de la vulnérabilité combinatoire

Pour finir, il est important de préciser que pour les deux plugins désormais implémentés, les nœuds des graphes doivent impérativement être numérotés en partant de 1, puis ne sauter aucune valeur.

Voici un jeu de données test (cliquer ici pour l’obtenir).

Références

Albert R., Albert I., Nakarado G., 2004, « Structural vulnerability of the North American power grid », Physical Review E, 69(2), pp. 1-4. http://arxiv.org/pdf/cond-mat/0401084v1.pdf

Crucitti P., Latora V., Porta S., 2006, « Centrality measures in spatial networks of urban streets », Physical Review E, 73(3), pp. 1-5. http://arxiv.org/pdf/physics/0504163.pdf

Gleyze J. F., 2005, La vulnérabilité structurelle des réseaux de transports dans un contexte de risque, Thèse de doctorat, Université Paris VII, Laboratoire COGIT-IGN, 826 p. http://tel.archives-ouvertes.fr/tel-00138991/

Lhomme S., 2012, Les réseaux techniques comme vecteur de propagation des risques en milieu urbain. Une contribution théorique et pratique à l’analyse de la résilience urbaine, Laganier R. & Serre D. (sous la dir.), Thèse de doctorat, Université Paris Diderot, 365 p. http://tel.archives-ouvertes.fr/docs/00/77/22/04/PDF/these_serge_lhomme.pdf

Wang J. W., Rong L. L., 2011, « Robustness of the western United States power grid under edge attack strategies due to cascading failures », Safety Science, 49(6), Elsevier, pp. 807-812. http://www.sciencedirect.com/science/article/pii/S0925753510002547

Serge Lhomme

Maitre de conférences à l’université Paris Est Créteil – Docteur en géographie de l’université Paris VII spécialité « dynamique des milieux et risques »

More Posts - Website

QGIS et l’analyse de graphe

Dans un précédent billet, il a été signalé qu’il existe peu d’outils intégrés au sein des SIG permettant d’analyser les réseaux. Néanmoins, les SIG libres peuvent être combinés à d’autres logiciels, programmes ou bibliothèques accroissant ainsi leur potentialité initiale. De plus, Grass (SIG libre qui, il est vrai, n’est pas des plus commodes à appréhender) permet d’effectuer différentes mesures sur des réseaux spatiaux.
Ce billet va présenter d’une part, comment réaliser un programme sur QGIS permettant d’analyser des réseaux spatiaux en ayant recours à des bibliothèques existantes d’analyse de réseaux (en l’occurrence, des bibliothèques développées en Python). D’autre part, ce billet va présenter un plugin réalisé à cet effet.

Utiliser une bibliothèque Python avec QGIS

Le SIG libre QGIS propose une console Python permettant d’interagir avec l’ensemble des outils de QGIS et par conséquent de développer de petits programmes automatisant plusieurs tâches. Dans les faits, il est plus souvent intéressant d’utiliser les nombreuses bibliothèques Python existantes afin d’étendre les capacités de QGIS. Par exemple, il existe de nombreuses bibliothèques permettant de réaliser de l’analyse de graphe (igraph, NetworkX…). Dans ce billet, nous allons utiliser NetworkX, bibliothèque développée originalement par Aric Hagberg, Dan Schult, et Pieter Swart (http://networkx.github.com/).

Pour utiliser une bibliothèque Python dans QGIS, il faut la télécharger (voici le lien pour la télécharger [1]), puis la dézipper et enfin la placer dans le répertoire des bibliothèques Python de QGIS. Ainsi, pour la version 1.8 (Lisboa) de QGIS, le répertoire se situe dans : « C:\Program Files\Quantum GIS Lisboa\apps\Python27\Lib\site-packages » (attention cette version de QGIS n’a pas été choisie par hasard, en effet elle implémente Python 2.7 qui est nécessaire pour la dernière version de NetworkX). Toutes les bibliothèques Python déjà présentes dans QGIS se trouvent dans ce répertoire (MapPlotLib, PyQt4…). La bonne complémentarité entre QGIS et Python est ici illustrée, car il n’y a plus rien à faire. Tout est installé et il est alors possible de réaliser de l’analyse de graphe avec QGIS, alors même qu’initialement ce n’était pas le cas.

Pour vérifier, il suffit d’ouvrir QGIS, d’aller dans « Extension -> Console Python », puis de taper les lignes de code suivantes :

>>> import networkx as nx  # Import de la bibliothèque NetworkX
>>> G=nx.Graph()  # Création d’un graphe G
>>> G.add_edges_from([(1,2),(1,3)])  # Création des arcs du graph G
>>> nx.clustering(G)  # Calcul les clustering coefficient des noeuds
{1: 0.0, 2: 0.0, 3: 0.0}
>>> nx.degree(G)  # Calcul les degrés des noeuds
{1: 2, 2: 1, 3: 1}

Analyser des graphes spatiaux dans QGIS avec Python

Cette étape est un peu plus complexe que la précédente. En effet, sans même avoir une grande connaissance de Python et de QGIS, il est possible d’utiliser la bibliothèque NetworkX dans QGIS en étudiant les tutoriels (voici le lien vers les tutoriels [2]). Faire interagir Python et QGIS ajoute un niveau de difficulté supplémentaire. Ainsi, partons du fait que deux couches ont été chargées dans QGIS. La première couche est composée de nœuds caractérisés par leur identifiant (vous pouvez télécharger un exemple ici). La deuxième couche est composée d’arcs caractérisés par leur identifiant, les deux identifiants des nœuds situés aux extrémités de l’arc et leur longueur (vous pouvez télécharger un exemple ici).

L’enjeu est alors de créer le graphe « NetworkX » à partir de ces deux couches. Pour cela, il faut récupérer les attributs des couches à l’aide de Python :

>>> canvas=qgis.utils.iface.mapCanvas() # Interaction avec la carte et les couches
>>> allLayers = canvas.layers() # Récupère les deux couches chargées dans QGIS
>>> couche_noeud = allLayers[0] # Récupère la couche supérieure
>>> couche_arc = allLayers[1]  # Récupère la couche inférieure
>>> provider = couche_arc.dataProvider() # Interaction avec la couche arc
>>> field=provider.fields() # Récupère les champs de la couche arc
>>> startind=provider.fieldNameIndex(‘depart’) # Index du champ « départ »
>>> endind=provider.fieldNameIndex(‘arrivee’)   # Index du champ « arrivée »
>>> distind=provider.fieldNameIndex(‘distance_K’) # Index du champ « distance_K »
>>> feat = QgsFeature()   # Permet de manipuler les entités des couches
>>> allAttrs = provider.attributeIndexes()   # Récupère les index des entités
>>> provider.select(allAttrs)   # Sélectionne toutes les entités 
>>> table=[]  # Création d’une table qui va récupérer les informations
>>> while provider.nextFeature(feat):  # Boucle parcourant l’ensemble des entités
…     attrs = feat.attributeMap()   # Récupère tous les attributs des entités
…     attrsstart=attrs.values()[startind].toInt()[0]  # Attribut nœud de départ
…     attrsend=attrs.values()[endind].toInt()[0]  # Attribut nœud d’arrivée
…     attrsdist=attrs.values()[distind].toInt()[0] # Attribut distance
…     table=table+[(attrsstart,attrsend,attrsdist)] # Ajoute les attributs récupérés à la table

Une fois les attributs des couches permettant de créer le graphe récupérés et après les avoir collectés dans un tableau, il est possible d’analyser ce graphe à l’aide la bibliothèque NetworkX :

>>> import networkx as nx  # Import de la bibliothèque NetworkX
>>> G=nx.Graph() # Création d’un graphe G
>>> G.add_weighted_edges_from(table) # Ajout des entités de la table au graphe G
>>> b=nx.betweenness_centrality(G,weight=’weight’) # Calcul de la centralité des noeuds du graphe G

Enfin, il est possible de récupérer les résultats obtenus pour les implémenter dans la couche nœuds :

>>> provider2 = couche_noeud.dataProvider() # interaction avec les données
>>> nbcol=int(provider2.fieldCount()) # Récupère le nombre de champs de la couche
>>> from PyQt4.QtCore import *   # Import de la bibliothèque PyQt4
>>> from PyQt4.QtGui import *   # Import de la bibliothèque PyQt4
>>> provider2.addAttributes([QgsField(« centralite », QVariant.Double)])
>>> taille=len(b)  # Récupère le nombre de nœuds
>>> couche_noeud.startEditing()  # Permet l’édition de la couche nœuds
>>> for k in range(taille):  # Boucle qui parcoure l’ensemble des nœuds
…      couche_noeud.changeAttributeValue(k,nbcol,b[k+1]) # Attribution des valeurs
>>> couche_noeud.commitChanges() # Enregistrement des changements

Toutes ces lignes de commande peuvent être récupérées et réorganisées pour créer un petit programme permettant de calculer les centralités des nœuds d’un graphe spatial (lien vers un fichier zip contenant un exemple ici) :

>>> execfile(« C:/Users/XXXXX/Desktop/centralite.py »)       # Appelle le programme à exécuter

Utiliser un plugin QGIS pour analyser des graphes spatiaux

Pour développer un plugin en Python dans QGIS, des connaissances supplémentaires en programmation sont nécessaires. Elles ne seront pas développées dans ce billet. Elles peuvent être trouvées sur ce site internet. Je vais plutôt introduire un petit plugin que j’ai développé. Il n’est pas très exhaustif, mais a pour objectif de montrer ce que l’on peut obtenir en combinant QGIS et une bibliothèque Python d’analyse de graphe.

Le dossier zip à télécharger se trouve ici. Une fois le dossier téléchargé et dézippé, il faut le placer dans le répertoire des plugins Python de QGIS. Pour la version Lisboa, le répertoire est le suivant : « C:\Program Files\Quantum GIS Lisboa\apps\qgis\python\plugins ». Le plugin est alors installé. Ensuite, il suffit de lancer QGIS pour que le plugin apparaisse dans « Extension -> Analyse de graphe ». S’il n’apparait pas, il faut aller dans « Extension -> Gestionnaire d’extensions » et cocher l’extension « Analyse de graphe version 0.1 ». Après avoir lancé le plugin, une boite de dialogue s’ouvre. Elle permet de sélectionner les nœuds et les arcs du graphe à étudier.

Après avoir validé cette première boite de dialogue, une deuxième boite de dialogue apparait. Elle permet de sélectionner les champs des extrémités des arcs, de préciser si le graphe est orienté ou non, de préciser si le graphe est pondéré ou non (et si oui, quel est le champ de cette pondération : distance, flux, importance…). Enfin, différents indicateurs sont proposés (ceux-ci peuvent être normalisés ou non).

Le tour est joué. Ainsi, un champ a été créé dans la couche des nœuds afin de stocker les résultats de la requête. Il est alors possible de réaliser une carte thématique finale.

Serge Lhomme

Maitre de conférences à l’université Paris Est Créteil – Docteur en géographie de l’université Paris VII spécialité « dynamique des milieux et risques »

More Posts - Website

Urban Network Analysis (UNA), la ToolBox qui enrichit l’analyse des réseaux sous ArcGis10x

Les outils du type de la ToolBox Urban Network Analysis (UNA) que nous présentons dans ce billet sont suffisamment rares pour être salués. Il n’existe pas, à notre connaissance, d’outil similaire intégré à un logiciel SIG – fusse t’il payant – permettant le calcul et la cartographie des mesures sur réseau.

Présentation :

UNA est donc un script fort utile car il augmente les possibilités de l’analyse des réseaux dans un SIG, en autorisant le calcul des mesures. Il a été développé par la Singapore University of Technology & Design en collaboration avec plusieurs structures telles que le City Form Lab du MIT, le Martin Center for Built Form and Land Use Studies, le Space Group de l’UCL et l’Human Space Lab, dans le cadre d’un programme de recherche présenté ici.

UNA s’appuie sur les résultats de différents travaux de recherche en analyse spatiale des réseaux (Martin & March 1971; Porta, Crucitti et al, 2005).

UNA est une Boite à outils (une Toolbox dans le jargon d’Esri) pour ArcGis10x – Service Pack 2 minimum – muni de l’extension Network Analysis.

Elle est téléchargeable ici ou , tout comme son fichier d’aide (UNA_help.pdf) sur lequel nous nous sommes appuyées pour rédiger le présent billet.

A noter que la Boite à outils est également fournie avec des fichiers exemples, ce qui facilite sa prise en main.

Développée à l’origine pour l’analyse de réseaux urbains non planaires, l’UNA peut être utilisée indifféremment pour l’analyse spatiale des réseaux planaires (de transport, par exemple) ou non planaire, mais seulement sous ArcGis10x. Les fichiers de départ nécessitent donc d’être aux format shapefiles (format natif de la gamme des produits ESRI dont l’extention est *.shp) et Network Analysis (*.nd). Si vous n’en disposez pas, de simples tableaux en (*.dbf par exemple) peuvent faire l’affaire, moyennant une transformation ultérieure (voir infra.)

Au sens de l’UNA, le fichier graphique du « réseau » est composé de trois éléments : des noeuds, des arêtes et des bâtiments. On ne présente plus les noeuds et les arètes (ou liens). Quant aux bâtiments, ils peuvent être assimilées à n’importe quel type d’unités spatiales surfaciques (polygones).

Intérêts de l’outil

Le premier intérêt de cet outil est de pouvoir combiner les trois types d’implantation spatiale (points, lignes, polygones) ainsi que leurs caractéristiques attributaires au sein d’une même analyse, ce qui n’est pas commun. La voie est ainsi ouverte vers une analyse de réseau dans un contexte multi niveaux, c’est-à-dire à différents échelons géographiques (l’ilôt, le quartier, la commune, etc.).

Le second intérêt porte sur le calcul de mesures pondérées puisque des caractéristiques attributaires des points et des polygones (des variables quantitatives de taille ou ordinales de rang, par exemple), et des lignes (longueur, etc.) peuvent être intégrées à l’analyse.

Le troisième intérêt porte sur les possibilités de visualisation cartographique des mesures.

L’UNA permet le calcul des cinq mesures classiques d’analyse des réseaux que sont : « R » (Reach) – « G » (Gravity) – « B » (Betweenness) – « C » (Closeness) et « S » (Straightness), en tenant compte de la position relative des nœuds les uns par rapport aux autres, selon une métrique discrète (topologique) ou continue (une distance kilométrique, par exemple). Elle restitue en retour différents fichiers des résultats.

Pré-requis :

1) Au préalable, il convient de vérifier que vous disposez bien de l’extension Network Analysis et que l’extension est bien « cochée ».

2) Ensuite, il faut dézipper l’UNA que vous avez téléchargé ( ici ou ) dans un répertoire de votre choix.

Deux fichiers et deux répertoires (qui contiennent des fichiers exemples) sont ensuite disponibles :

Urban_Network_Analysis_Tool.tbx qui est la boite à outils et son fichier d’aide. Tous deux s’installent automatiquement à la racine.

– les répertoires sont  <Code> et <TextFile>

<Code> contient les codes de sept scripts (rédigés en Python) qui permettant de générer les différents fichiers de résultats tels que, par exemple, le calcul de la matrice d’adjacence (tableau croisé de valeurs binaires ou la modalité ‘1’ correspond à l’existence d’un lien, sinon ‘0’) nommée : Adjacency_List_Computation ; ou encore les résultats des différentes mesures calculées. La centralité est stockée dans le fichier Centrality_List_Computation.

<TextFiles> contient deux sous-répertoires de fichiers d’exemples. Le premier contient les shapefiles d’unités spatiales surfaciques. Il est intitulé  Test_Buildings ; le second contient   les shapefiles des  ‘unités spatiales lineaires, nommé Test_Networks.

3) Enfin, il vous faut charger la Boite à Outils dans la fenêtre « Toolbox » de ArcGis10x. Pour ce faire, un click-droit sur la fenêtre des ToolBox / ajouter une ToolBox puis parcourir votre arborescence pour sélectionner l’UNA. Elle s’installe ensuite automatiquement à la fin de votre liste d’outils.

Préparation des données :

Deux fichiers minimum sont nécessaires : un fichier de points (ou un fichier de polygones) au format (*.shp) et un fichier de lignes au format (*.nd).
Il est important que tous vos fichiers soient stockés dans un même répertoire.

Si vous disposez de ces deux fichiers,vous pouvez aller directement au calcul des mesures, sinon, vous pouvez les générer en suivant la procédure suivante :

1) Sous ArcGis (ArcCatalog ou ArcMap), le fichier de points (*.shp) peut être généré à partir d’un tableau importé d’un tableur, puis converti par la procédure d’export en (*.dbf), avec un click-droit.

Ce tableau doit contenir au minimum 3 colonnes que sont : ID ; X_COORD ; Y_COORD. La colonne ‘ID’ correspond au numéro d’identification unique du point) ; la colonne ‘X_COORD’ à la valeur de la coordonnée en X du point dans un système de projection ; et ‘Y_COORD’ à la valeur de la coordonnée en Y du point dans le même système de projection, bien entendu.

En l’absence de fichier de points, si vous disposez d’une couche de polygone (*.shp), vous pouvez soit vous servir directement de celle-ci, soit calculer au préable les barycentres (ou centroïdes) des polygones (AddXYdata ou DataManagementTool/FeatureToPoint) puis en faire un fichier de points.

2) Le fichier de lignes (qui correspond aux arêtes du réseau) peut être également généré sous ArcGis (ArcCatalog ou ArcMap). Il importe beaucoup que celui-ci soit un fichier pour Network Analysis, c’est-à-dire muni de l’extention (*.nd).

Il est possible de générer ce fichier à partir d’un (*.shp) ou de fichiers externes importés (*.dwg/dxf) sous ArcCatalog, en les convertissant : click-droit sur le fichier de lignes (*.shp) / choisir « nouveau Dataset » et suivre les instructions.

Enfin, en l’absence de fichier de lignes, il est possible d’en créer un à partir du fichier de points précédent (voir supra, 1), en utilisant par exemple, l’outil  Point To Ligne pour obtenir un (*.shp), puis le convertir en (*.nd).

Calcul des mesures

Pour calculer les différentes mesures proposées par l’UNA, il est nécessaire de lancer l’outil (en double-cliquant dessus).

La fenêtre de menu suivante s’ouvre ensuite. Elle vous propose de

– sélectionner les deux fichiers préparés lors de l’étape précédente : celui des points (ou des surfaces) et celui des arètes (ou des lignes).

– sélectionner la mesure que vous souhaitez calculer.

Fenêtre de menus de l'UNA

Le descriptif des mesures proposé est le suivant :

1) la “Reach measure » telle que définie par Sevtsuk (2010) permet de sélectionner les nœuds (ou polygones) situés dans un certain rayon (r) de voisinage sur le réseau et de les dénombrer. Le voisinage est défini sur une base topologique (selon une métrique discrète). Il génère une zone tampon circulaire. Le voisinage fait référence ici au concept d’occasions interposées appliquées au réseau, fondé sur le principe du nombre de nœuds à franchir pour aller d’un point à un autre.

2) La « Gravity measure », ou « mesure gravitaire » telle qu’introduite par Hansen (1959), est une mesure importante dans l’analyse d’accessibilité dans un réseau de transport. Elle fait référence à l’attractivité des noeuds qui est inversement proportionnelle à la distance qui sépare chacun d’entre eux. Le modèle théorique sous-jacent est ainsi le modèle gravitaire, utilisé de manière pondérée ou non en fonction des caractéristiques des nœuds (taille) ou des liens (type de distance, par exemple). L’exposant  de la distance – qui prend en compte la rugosité de l’espace (ou friction) – peut être contrôlé en fonction du mode ou du temps de transport (à pied, en véhicule…). Les auteurs suggèrent, par exemple, de le porter à 0.1813 pour un déplacement à pied. Sans indication de l’utilisateur, cet indice est porté par défaut à 1.

3) La « Betweeness measure » ou « mesure de centralité d’intermédiarité » retient l’attention de nombreux chercheurs actuellement. Elle identifie les noeuds qui permettent une circulation rapide au sein du réseau, en limitant la distance à parcourir. Les supprimer conduit donc à réduire l’accessibilité de certains lieux. L’UNA utilise l’algorithme développé par Brandes ( 2001) qui renseigne sur « le potentiel »  des nœuds. On regrettera l’absence de cette mesure au niveau des liens.

4) La « Closeness measure » ou « mesure de centralité de proximité » d’un nœud correspond à l’inverse de la distance moyenne cumulée nécessaire pour atteindre ce nœud depuis n’importe lequel des autres nœuds du réseau si aucun critère de distance n’est défini. La centralité de proximité est donc l’inverse de l’indice de Shimbel. Il présente l’avantage d’être normalisé (de 0 à 1) avec 1 comme accessibilité maximale.

5) la « Straightness measure »  (Vragovic, Louis, et al. 2005 ; Porta, Crucitti et al.2005).

On notera que l’ensemble des résultats peut ensuite être normalisé.

La formulation des différentes mesures, fournie dans le fichier d’aide, est synthétisée dans le tableau ci-après :

Formule des mesures proposées par l'UNA

Pour en savoir plus sur les mesures locales de voisinage ou d’ensemble (globale) d’un réseau et leur utilisation, vous pouvez toujours lire et relire la synthèse du groupe f.mr. rédigée par César Ducruet (2010). Elle est disponible ici.

Cartographie des résultats

Dès lors que la mesure choisie a été calculée, la Betweenness par exemple, elle est stockée dans un fichier de résultat. La colonne ‘B’ correspondante vous permet ensuite de visualiser directement le résultat cartographique, puisque vous êtes dans un SIG…La figure suivante, extraite du fichier d’aide, présente un exemple de cartographie de la ‘B’.

Prise en compte des freins aux déplacements

La prise en compte de freins aux déplacements ou de contraintes à la circulation dans le réseau peut être réalisée. Elle nécessite toutefois une structuration particulière du fichier de points, par la définition « d’attributs d’impédance », sous Network Analysis.

Export des résultats

L’UNA fournit trois familles de fichiers de résultats forts utiles pour des analyses au sein de ArcGis10x ou avec d’autres outils.

1) Les fichiers des résultats du calcul des mesures, disponibles en (*.dbf), contiennent de nouvelles variables. Comme indiqué précédemment, elles sont intitulées R, G, B, C, S et correspondent respectivement aux mesures Reach, Gravity, Betweenness, Closeness et Straightness. C’est donc à partir de ce fichier que les cartes sont réalisées et « symbolisées » dans ArcMap.

2) Les matrices d’adjacence (Bldng_Adj*.dbf) qui contiennent les résultats des requêtes réalisées sur les noeuds obtenues suite au calcul de l’une ou de l’autre mesure.

3) Les graphes (*.shp) associés à la matrice d’adjacence qui permettent de spatialiser les relations de voisinage, de visualiser leurs différents attributs, etc.

A noter qu’une matrice d’adjacence et un graphe associé sont générés à chaque fois qu’une requête est réalisée sur les noeuds et/ou sur les lignes.

Pour finir, il s’agit d’un bel outil, simple d’utilisation mais qui nécessite toutefois une bonne compréhension de la signification des indicateurs associée à une bonne maîtrise d’ArcGis. Une compétence sur l’utilisation de Network Analysis est un vrai plus.

Références mentionnées

Brandes, U. , 2001,  A faster Algorithm for Betweenness Bentrality. Journal of Mathematical Sociology, Vol 25(2), 163-177.

Hansen, W. G. , 1959, How Accessibility Shapes Land Use. Journal of the American Planning Association, 25(2), 73-76.

Hillier & Hanson, 1984, The Social Logic of Space. Cambridge: Cambridge University Press.

Porta, Crucitti et al, 2005, The network analysis of urban streets: a primal approach. Environment and Planning B, 35(5), 705-725.

Sevtsuk, A. , 2010, Path and Place: A Study of Urban Geometry and Retail Activity in Cambridge and Somerville, MA. PhD Dissertation. Cambridge: MIT

Vragovic, I., Louis, E., & Diaz-Guilera, A., 2005, Efficiency of information transfer in regular and complex networks. Physics Review E., 71(026122).

Pour en savoir plus sur le projet:

http://cityform.mit.edu/projects/urban-network-analysis.html

 

Françoise Bahoken

Géographe et cartographe, je suis accueillie en détachement depuis décembre 2010 à l’IFSTTAR ; je suis par ailleurs membre-associée de l’UMR 8504 Géographie-Cités.

More Posts

La détection de communautés avec Pajek 3.6

La toute nouvelle version de Pajek réserve de bonnes surprises aux explorateurs de graphes complexes.

Les nouveautés dont il est question ici ont été ajoutées au logiciel entre Mai 2012 et Septembre 2012.

Deux algorithmes de détection de communautés ont été implémentés:

Network/Create Partition/Communities

-Louvain method

Le premier algorithme s’appuie sur la méthode de Louvain [1]. Cette méthode est bien adaptée aux graphes valués, y compris s’ils sont signés. Comme pour beaucoup d’algorithmes de clustering (ou détection de communautés), la qualité du partitionnement est optimisée par maximisation de la fonction de modularité [2].  L’indice de modularité mesure la différence entre la densité d’un graphe (ou sous-graphe) donné et la densité d’un graphe aléatoire possédant les même caractéristiques (nombre et poids des liens).

Principe

Maximiser la fonction de modularité pour partitionner un graphe revient à s’assurer que le nombre et le poids des liens est plus important à l’intérieur des partitions qu’entre les partitions. Pour l’exprimer autrement, il faut que la densité intra-communautaire dépasse la densité inter-communautaire.

La méthode de Louvain est « bottom-up » et multi-niveau. Au lancement de l’algorithme, tous les sommets appartiennent à une partition différente. Ils sont regroupés, par itération, dans des partitions de modularité optimale. Arrivé à une première situation d’optimum, le processus se poursuit au niveau supérieur : chaque partition est traitée comme un sommet et ainsi de suite. L’opération continue jusqu’à ce qu’il n’y ait plus aucun gain de modularité possible.

Avantage

Cette approche permet d’échapper au défaut de la modularité connu sous le nom de « limite de résolution » [3]. Si le fait d’associer une partition du premier niveau avec une autre devait engendrer une perte de modularité, cette association ne serait pas réalisée par la méthode de Louvain. De ce fait, les petites structures dotées d’une bonne modularité ne risquent pas d’être noyées dans des structures plus grandes et moins significatives. C’est ce qu’il risquerait d’arriver avec un algorithme « top-down » appliqué à un trop grand graphe.

Inconvénient

Le défaut connu de la méthode provient de la sensibilité du résultat à l’ordre de traitement des sommets. L’ordre dans lequel les sommets sont considérés a une influence sur le partitionnement. De ce fait, il peut être utile de procéder à des permutations pour s’assurer de la robustesse du résultat [4]. L’autre piste est d’imposer un ordre qui tienne compte des caractéristiques du réseau [5].

En plus de l’algorithme classique de Louvain,  Pajek propose une variante (avec rafinement multi-niveau) qui tend à donner de meilleurs résultats pour les grands graphes [6]. La méthode de Louvain est également disponible dans les logiciels Gephi et NetworkX ainsi que dans la librairie igraph du logiciel R, l’algorithme étant plus rapide sur Pajek que sur Gephi.

-VOS Clustering

La méthode de détection de communauté VOS optimise une autre fonction de qualité que la fonction de modularité: la fonction dite de VOS clustering [7]. Il s’agit d’une variante pondérée de la modularité qui a été pensée pour le traitement des données bibliométriques. Cette technique donne moins de poids qu’une méthode classique aux sommets à très haut degré. Contrairement à la méthode de Louvain, elle dépend d’un algorithme « « top-down« . Pour cette raison, le seul moyen d’échapper au problème de limite de la résolution consiste à faire varier le paramètre de résolution.

Cette méthode convient en fait aux graphes pour lesquels la valeur des liens rend compte de l’importance de la similarité entre les sommets. Ses concepteurs ont souhaité que soit associé à l’algorithme de clustering un algorithme de visualisation faisant appel au même principe mathématique. Cette approche unifiée permet d’éviter qu’ils y ait des incohérences entre le résultat du processus de clustering et la représentation graphique de ce résultat.

Tout naturellement, Pajek propose aussi l’algorithme de visualisation VOS. Il fonctionne sur le principe de l’optimisation de la fonction de qualité dite de VOS mapping. Étant donné la proximité des deux méthodes de détection, il n’est surement pas impensable d’utiliser le VOS mapping pour représenter un graphe partitionné suivant la méthode de Louvain.

Aller plus loin

Une comparaison des deux méthodes de détection est accessible sur: http://mrvar.fdv.uni-lj.si/pajek/community/LouvainVOS.htm

Des précisions sur l’utilisation de ces méthodes avec Pajek se trouvent aussi sur le site: paramétrage, évaluation de la qualité du partitionnement, représentation graphique du résultat.

Depuis l’interface graphique de Pajek, un raccourci vers le logiciel VOSviewer offre la possibilité d’obtenir rapidement une représentation  élégante du partitionnement.

Ci-dessous, un exemple à partir des données de collaborations scientifiques néerlandaises (Source: Web of Science 2006-2008):

Références
[1] V. Blondel, J.-L. Guillaume, R. Lambiotte, et E. Lefebvre, « Fast unfolding of communities in large networks », Journal of Statistical Mechanics, vol. 2008, no 10, sept. 2008.
[2] Il existe de nombreuses autres méthodes de détection de communautés. Pour avoir un aperçu de l’existant: S. Fortunato, « Community detections in graphs », Physics and Society, vol. 2010, no 486, p. 75-174.
[3] S. Fortunato et M. Barthélemy, « Resolution limit in community detection », Proceedings of the National Academy of Sciences of the United States of America, vol. 104, no 1, p. 36-41, janv. 2007.
[4] V. Blondel, G. Krings, et I. Thomas, « Régions et frontières de téléphonie mobile en Belgique et dans l’aire métropolitaine bruxelloise », Brussels studies, no 42, oct. 2010.
[5] P. De Meo, E. Ferrara, G. Fiumara, et A. Provetti, « Generalized Louvain method for community detection in large networks », under review, 2011. URL: http://www.emilio.ferrara.name/wp-content/uploads/2011/07/isda2011-k-path.pdf
[6] R. Rotta et A. Noack, « Multilevel local search algorithms for modularity clustering », Journal of Experimental Algorithmics, vol. 16, no 2.3, 2011.
[7] L. Waltman, N. J. van Eck, et E. C. M. Noyons, « A unified approach to mapping and clustering of bibliometric networks », Journal of Infometrics, vol. 4, no 4, p. 629-635, oct. 2010.