Journal le dessous des cartes

Posté par  . Licence CC By‑SA.
Étiquettes :
37
14
août
2015

Sommaire

Cher journal,

J'aime bien lire sur ce site les histoires de création de carte comme celle-ci ou celle-là, alors je vais te conter la mienne.

quelle est la question ?

Partant d'une question pertinente : "Est-il facile de se ravitailler en GPL en France", j'en suis venu à une question idiote : "Est ce qu'il y a une grande différence des prix du carburant en fonction d'où on habite".

J'imagine une carte de France des points de vente, catégorisés en "bien en dessous de la moyenne" en vert, "dessous de la moyenne" en vert clair, "dans la moyenne" en orange/jaune, "au dessus de la moyenne" en rouge clair, "bien au dessus de la moyenne" en rouge.

Une carte par type de carburant, en utilisant les derniers prix constatés.

obtenir les données

source

Le site prix des carburants du ministère de l'économie, de l’industrie et du numérique (Oui oui, tout ça), permet de connaître les prix en points de vente (par la suite PdV). Ceux-ci ont obligation de déclarer à l'état les changements de prix.

C'est bien évidement parce que je connais l'existence de ce site que ces questions plus ou moins tordues me sont venues.

La section open data permet de télécharger l'historique de l'ensemble des déclarations des points de ventes.

traitement

Les données sont au format XML et plutôt volumineuses. Sur les sept premiers mois de l'année, un peu plus de 2 millions d'enregistrements pour plus de 120MB. J'aurai dû sortir le parseur adéquate genre xml pull parser. Mais le fait qu'il soit bien indenté m'a permis de faire un traitement un peu plus brutal en ligne à ligne avec awk. Sans plus d'optimisation, les données exploitables sont produites en quelques secondes mais avec un léger embonpoint mémoire d'environ 100MB. Le code, plutôt quick and dirty, est déposé ici et invoqué avec cela.

Les données exploitables sont formatées dans des fichiers csv, un par type de carburant, dont les colonnes sont l'identifiant du PdV, les coordonnées latitude/longitude du PdV, le dernier prix déclaré pour ce type de carburant ; une ligne par PdV.
Un fichier de statistiques contient la moyenne et l'écart type pour chaque type de carburant.

statistiques

Comme je ne voyais pas de raison pour qu'il en soit autrement, j'ai supposé que les prix suivaient une distribution en Gaussienne et donc que je pouvais procéder à un découpage en intervalles comme proposé par la théorie (μ est la moyenne et σ l'écart type):
intervals

Ce qui me fait mes cinq intervalles :

  1. x < μ−2σ,
  2. μ−2σ ≤ x < μ−σ,
  3. μ−σ ≤ x ≤ μ+σ,
  4. μ+σ > x ≥ μ+2σ,
  5. x > μ+2σ

Dans les faits, l'hypothèse est à peu près vérifiée.

construire la carte

fond de carte

Wikipedia fournit une grande variété de carte de France. Je suis parti de celui-ci, un peu raboté et converti en PNG:
régions et départements

placer les points

Les données de base permettent de travailler avec des coordonnées GPS (latitude/longitude) ou avec les communes (code postal). J'ai choisi de travailler avec des coordonnées GPS. Il va donc falloir trouver les emplacements des points d'intérêts sur une carte en pixels.

projection

La terre étant ronde, pour la représenter sur une carte plate, il faut faire une projection. Et cette projection est toujours imparfaite : si une partie de la carte est correcte, le reste est forcément déformé. A ce titre, le fond de carte présenté sur le Wikimedia Commons est intéressant car on voit que l'échelle varie selon la latitude.
échelle

En effet, pour la France, on utilise la projection Lambert-93 (avec comme parallèles "standards", CC44 et CC49). Elle a le mérite de ne pas créer de déformation horizontales (le long d'une latitude, d'est en ouest) mais crée des déformations verticales (de long d'une longitude, du nord au sud).

J'ai donc glané sur le web une fonction effectuant la conversion de WGS84 (système GPS) vers Lambert-93. Le code était en php, je l'ai réécrit en awk pour mon usage. Il est déposé ici pour les curieux. Attention, ça pique un peu. Je conseille de repartir d'une implémentation propre si vous avez ce besoin. L'article wikipedia donne les pointeurs chez l'IGN pour avoir les paramètres et les formules officiels.

zoom et centrage

Les coordonnées Lambert-93 représentent une position x,y dans un plan dont l'origine (centrage) et l'échelle (zoom) sont différentes de celles mon fond de carte.
Pour déterminer la transposition que je dois faire, il me faut connaître les coordonnées de trois points dans chaque plan.
Sur le fond de carte, il est ardu de positionner précisément des villes à l'œil nu. J'ai donc choisi de faire correspondre les points les plus au bord.
En utilisant le site suivant, j'obtiens:

             |        lambert       | fond de carte
pointe ouest | x= 124754, y=6831715 |  X=  0, Y=293
pointe nord  | x= 668091, y=7110469 |  X=557, Y=  0
pointe est   | x=1082912, y=6886783 |  X=988, Y=233
pointe sud   | x= 714399, y=6148363 |  X=608, y=951

J'utilise les trois premiers pour paramétrer la transposition et le quatrième pour vérifier. Le code en python utilise numpy pour la résolution du système d'équations. Il est disponible ici. La vérification n'est pas très concluante d'ailleurs mais poursuivons.

dessiner

processus

Je vais dérouler la procédure suivante:

  • constantes : type de carburant, moyenne, écart type, couleurs, couleur de fond
  • charger les pixels du fond de carte depuis PNG
  • pour chaque PdV:
    • transposer les coordonnes GPS en pixel dans le plan de la carte
    • si le pixel est couleur fond (je ne veux pas abîmer les bordures)
      • calculer la couleur, changer la couleur
  • dumper les pixels vers PNG

outil

J'ai choisi d'utiliser python et son module PIL.
Le code, à nouveau quick and dirty, est déposé ici pour les curieux.

résultat

L'image produite par tout ceci est celle-ci pour le SP95:
SP95

Première remarque: les points sont un peu trop petits pour que la visualisation se fasse bien. Il y a aussi beaucoup (trop ?) de points jaunes, c'est à dire dans la moyenne.
Ensuite, la projection est clairement fausse ! Les points rouges qui devraient se trouver sur Paris sont légèrement décalés. Plus on descend vers le sud, plus le décalage est important pour finir avec une multitude de points dans l'eau :(
Pour me rassurer, je pense que la projection du fond de carte n'est pas en lambert93, contrairement à ma projection des points d'intérêt; d'où le décalage.

Mais quel était l'intérêt de la manœuvre ?
Partant d'une question plutôt anodine, j'ai du pas mal me triturer le cerveau pour arriver à produire quelque chose. J'imagine que les personnes dont c'est le métier sont mieux outillés et plus expérimenté et peuvent réaliser cette tâche de manière relativement triviale. Cela m'a mis en appétit pour m'intéresser plus longuement à OSM et aux logiciels de GIS. Sûrement pour le prochain épisode !

  • # Le rabot, cet outil si magique

    Posté par  . Évalué à 3.

    un peu raboté et converti en PNG

    Tu as raboté la Corse, tu vas avoir des problèmes :/ Je ne m’attarderais pas sur le fait que cette carte ne fait pas référence à la Réunion ou à la Guadeloupe mais celles-ci devaient déjà être rabotées avant que tu n’interviennes…

    • [^] # Re: Le rabot, cet outil si magique

      Posté par  . Évalué à 3.

      Oui, c'était le clin d'œil derrière cette expression, en effet.
      Mais bon pour la majorité des gens il faut prendre le ferry, ça fausse pas mal les statistiques :)

  • # Noooooooooooooooooooon

    Posté par  (site web personnel) . Évalué à 10.

    Je suis en vacances, j'essaye de me sortir la tête de mes cartes et voila qu'on arrète pas d'en parler sur LinuxFr …

    Tu as fait ce qui me semble presque la partie la plus dure : Transformer le XML en quelque chose d'exploitable avec en plus, magie, des coordonnées.

    Tu te lance ensuite, plus ou moins a la main, dans une transformation de projection, qui est casse gueule, et ou tu fait un raté.

    C'est le boulot d'un SIG (GIS), et QGIS en est un très bon, libre et très actif, qui ne fera qu'une bouchée de pain d'un si «simple» problème.

    Trouver des données

    Tu as déjà les données sur les stations de carburants, avec des coordonnée, tu as donc la moitié du travail.

    Il faut trouver un fond de carte géo-référencé. La gamme «GEOFLA» de l'IGN, fournie dans une licence libre dans le câdre de l'OpenData, fournira un très bon point de départ : http://professionnels.ign.fr/geofla
    QGIS (ou un autre SIG) sera rapidement capable d'en tirer une carte au moins équivalente a celle que tu as pris ici, mais avec un géo-référencement direct.

    Suivant la finesse de ce que tu veux faire, d'autres données de l'IGN peuvent être utilisé (voir leur Catalogue et surtout ce qui sont des données Libre), que l'on pourra compléter par des donnée OpenStreetMap si l'on veux aller dans le détail (Dont Géofabrik fournis des exports intéressants)

    Mixer les données

    Difficile d'aller plus loin sans faire un court complet d'utilisation de QGIS (ce qui risque d'être un peu long), mais globalement, tu y intègre tes jeux de donnée et c'est dans l'interface même du SIG que tu paramètre tes filtrage et tes mises en forme.

    QGIS (et probablement tout SIG qui se respecte) prendra de lui même en charge les modèles statistiques courants.

    Vers l'infini et au delà

    Avec peu être un peu plus d'outils (libre) il est possible de faire beaucoup beaucoup plus, et de monter un véritable projet.

    L'association PostgreSQL, PostGis, QGIS permet de faire beaucoup de choses très intéressantes.
    Il faudrait que je prenne le temps de transcrire un peu mon aventure (qui dure depuis plus de deux ans maintenant) a ce niveau.
    Mais vu que je suis en vacances, et que le but est de me sortir un peu l'esprit de tout cela, ça sera pour plus tard.

    • [^] # Re: Noooooooooooooooooooon

      Posté par  (site web personnel, Mastodon) . Évalué à 10.

      Pour les curieux, avec QGIS 2.10 que je découvre pour l'occasion:

      • télécharger les limites des communes métropolitaines ici
      • extraire le fichier, et avec QGIS, ouvrir LIMITE_COMMUNE.SHP
      • télécharger les prix du carburant qui vous intéresse, et en créer un fichier contenant les 3 colonnes: latitude longitude prix
      • ajouter un layer à partir du fichier en question: Layer > Add Layer > Add delimited text layer…
      • choisir le bon délimiteur, faire correspondre latitude et longitude avec vos colonnes, confirmer
      • choisir la projection EPSG:2154, RGF93 / Lambert-93
      • pour le color coding: clic droit sur le layer > Style > Graduated, puis choisir la 3ème colonne (celle du prix, y'en a peut-être un ou deux qui suivent?), l'échelle de couleur et la façon de générer ces dernières à partir de la valeur des prix

      Bon, vous aurez reconnu dans mon explication une implémentation de la méthode la R.A.C.H.E., c'est juste pour mettre le pied à l'étrier de ceux qui hésite, pendant que ceux qui savent vraiment faire sont en vacances :)

    • [^] # Re: Noooooooooooooooooooon

      Posté par  . Évalué à 3.

      Le décalage des points vient du fait que le fond de carte n'est pas projeté en Lambert 93, mais en projection cylindrique équidistante. A part ça, le cheminement est bon, à mon avis.
      Mais avec QGis, ce sera plus simple, c'est sûr !

      • [^] # Re: Noooooooooooooooooooon

        Posté par  . Évalué à 3.

        Merci, c'est une super info. Comment as tu deviné ?
        Je vais me mettre au gis mais en attendant, sais tu où je peux trouver l'algo de projection ?

  • # R et maptools

    Posté par  . Évalué à 7.

    L'utilisation de R et de la librairie maptools me semble tout à fait adaptée dans ce cas.
    Il y a des fonds de cartes inclus par défaut. Tu peux changer de projection en appelant une seule fonction ou simplement demander dans quel projection est telle carte ou tel jeu de données.
    En utilisant R, tu as également toute la gamme d'outils statistiques qui peuvent de permettre de faire une partition un peu plus intéressante de tes données. Ou alors de la faire plus rapidement si tu veux garder le même système.

    Par contre, j'ai voulu vérifier d'où venaient les fonds de carte et j'ai pas réussi à retrouver. La librairie est sous licence GPL>=2 mais pas moyen de trouver la source exacte.

    • [^] # R et ggplot2

      Posté par  . Évalué à 1.

      Gazole

      library(XML)
library(ggplot2)
library(maps)
library(scales)
library(RColorBrewer)

day <- "20150827"
baseurl <- "http://donnees.roulez-eco.fr/opendata/jour/"

# Récupération des données
temp <- tempfile()
download.file(paste0(baseurl, day), temp)
xmlsource <- unzip(temp, paste0("PrixCarburants_quotidien_", day, ".xml"))
unlink(temp)
parsed <- xmlParse(xmlsource)

gazole <- data.frame(
    id = xpathSApply(parsed, "//prix[@nom='Gazole']/../@id"),
    cp = xpathSApply(parsed, "//prix[@nom='Gazole']/../@cp"),
    latitude = as.numeric(xpathSApply(parsed, "//prix[@nom='Gazole']/../@latitude")) / 100000,
    longitude = as.numeric(xpathSApply(parsed, "//prix[@nom='Gazole']/../@longitude")) / 100000,
    prix = as.numeric(xpathSApply(parsed, "//prix[@nom='Gazole']/@valeur")) / 1000,
    date = as.Date(xpathSApply(parsed, "//prix[@nom='Gazole']/@maj"))
)

# Ajout départements. Cas particuliers : Corse et COM
gazole$dep <- ifelse(substr(gazole$cp, 1, 2) == "97" | substr(gazole$cp, 1, 2) == "20",
                     substr(gazole$cp, 1, 3),
                     substr(gazole$cp, 1, 2)
                     )
gazole$dep[which(gazole$dep == "200" | gazole$dep == "201")] <- "2A"
gazole$dep[which(gazole$dep == "202" | gazole$dep == "206")] <- "2B"

gazole$dep <- as.factor(gazole$dep)

# suppression Réunion (et mauvaises loc ?)
gazole.metro <- subset(gazole, longitude < 10 & latitude > 40)

# Moyenne prix et intervalle de confiance
mean(gazole.metro$prix)
confint(lm(gazole.metro$prix ~ 1))

# GLM : différences entre départements
res <- glm(data = gazole, prix~dep)
summary(res)

# Graphiques boxplot par département et par prix
#+ fig.width=15
ggplot(gazole, aes(dep, prix)) + 
    geom_boxplot() +
    stat_summary(fun.y = mean,  geom = "point") +
    ggtitle(paste("Gazole", day)) +
    xlab("Département") +
    ylab("Prix (€)") +
    theme(axis.text.x=element_text(angle = 90, hjust = 0))

#+ fig.width=15
ggplot(gazole, aes(x=reorder(dep, prix, FUN=mean), prix)) + 
    geom_boxplot() +
    stat_summary(fun.y = mean,  geom = "point") +
    ggtitle(paste("Gazole", day)) +
    xlab("Département") +
    ylab("Prix (€)") +
    theme(axis.text.x=element_text(angle = 90, hjust = 0))

# Histogramme
ggplot(gazole, aes(x = prix)) + 
    geom_histogram() +
    ggtitle(paste("Gazole", day)) +
    xlab("") +
    ylab("") 

# Carte
#+ fig.height=10
france.map <- map_data("france")
ggplot(france.map, aes(long,lat)) +
    geom_polygon(aes(group = group), fill = NA, colour = "grey60") +
    geom_point(data = gazole.metro, aes(longitude, latitude, color = cut(prix, quantile(prix, seq(0,1,0.2))))) +
    scale_colour_manual("Prix (€)", values=rev(brewer.pal(5, "RdBu")), guide=guide_legend(reverse = TRUE)) +
    coord_fixed(ratio = 1.45) +
    ggtitle(paste("Gazole", day)) +
    xlab("") +
    ylab("")

  • # Remarques en vrac

    Posté par  (site web personnel) . Évalué à 8.

    Dans les faits, l'hypothèse est à peu près vérifiée.

    Alors ça … il y a 12.000 distributions qui ressemblent à des gaussiennes quand tu regardes un nuage de points. Par exemple tu ne peux souvent pas trop faire la différence avec une variable loi log-normale (variable dont le logarithme suit une loi normale), et au lieu du logarithme tu peux prendre n'importe quel “homémorphisme concave ou convexe de R sur un intervalle de R et dont la dérivée à la valeur moyenne vaut 1” pour transformer un nuage de points gaussien: on voir toujours un nuage de points gaussien.

    Ensuite l'affirmation “c'est proche d'une distribution gaussienne” est en soi assez vide, ce qui compte c'est de savoir si c'est plus proche de telle ou telle autre distribution de référence.

    Pour y voir plus clair, tu peux comparer les distributions cumulées et les fonctions quantiles: pour la gaussienne tu peux la calculer avec ta bibliothèque de calcul scientifique préférée, pour la distribution c'est l'intégrale de l'histogramme ou, pour le dire autrement, les probabilités des sous-niveaux. Tu peux simplement superposer leurs graphes. La comparaison des fonctions quantiles est en général très parlante on réalise en général un Q-Q-plot, c'est à dire la fonction association les quantiles d'une distribution à celle de la seconde. Si je note PHI et PSI les deux distributions cumulées, c'est le graphe de PSI ROND (INVERSE PHI) et la fonction (INVERSE PHI) est dans toutes les bibliothèques de calcul scientifique. (C'est chiant qu'on ne puisse pas poster du UNICODE …)

    Cf. https://en.wikipedia.org/wiki/Q–Q_plot

    D'autres calculs importants à faire sont les moments d'ordre supérieur (skewness et kyrtosis pour 3 et 4) qui vont t'aider à distinguer ta distribution de la distribution gaussienne, là encore ta bibliothèque de calcul scientifique devrait tout savoir faire.

    Actuellement les prix sont plutôt modélisés par des variables log-normales (qui ne prennent que des valeurs positives, propriété importante des prix!) mais certaines personnes pensent que cette distribution devrait avoir une queue épaisse, ce qui n'est pas le cas des lois log normales.

    Elle a le mérite de ne pas créer de déformation horizontales (le long d'une latitude, d'est en ouest) mais crée des déformations verticales (de long d'une longitude, du nord au sud).

    Je ne sais pas trop ce que tu appelles une “déformation” mais les cartes qu'on utilise en géographie sont en général des projections conformes c'est à dire que les angles mesurés sur la carte sont ceux observés sur la terre. C'est une propriété cruciale pour la navigation, par exemple! Les distances elle n'ont rien à voir. Avec Lambert, si je mesure une distance le long d'une longitude (direction SN) et reporte la même distance carte aux même longitudes mais une autre longitude, alors la distance réelle ne change pas. Si je fais la même opération en inversant les rôle des longitudes et latitudes alors les distances réelles changent!

    • [^] # Re: Remarques en vrac

      Posté par  . Évalué à 0.

      Salut,

      Je vais faire plus court.

      Dans les faits, l'hypothèse est à peu près vérifiée.

      Pas en temps de trajet aoûtiens !

      Sur l'année je ne sais pas, mais je pense qu'en été, les zones d'approvisionnement comme les aires d'autoroutes changent leurs tarifs. Soyons bien d'accord, il s'agit de tarif, pas de passages.

      Une info intéressante serait d'afficher en "gros" les outliers. Une règle au doigt mouillé est 3-sigmas, ou 2.5 sigmas corrigés (on normalise sur la MAD et non la moyenne).

      En tout cas, jolies cartes :)

      Matricule 23415

  • # Marque

    Posté par  (site web personnel) . Évalué à -6.

    le dessous des cartes

    toi, tu cherches les ennuis

    • [^] # Re: Marque

      Posté par  . Évalué à 4.

      C'est un hommage.
      J'aime beaucoup cette émission.

    • [^] # Re: Marque

      Posté par  (site web personnel) . Évalué à 7.

      C'est une très vieille expression, déjà présente dans le dictionnaire de l'Académie française, 4e édition (1762) : « On appelle "Le dessous des cartes", La carte ou les cartes qui sont au-dessous du jeu de cartes après qu'on a coupé. » (source). J'imagine que le sens a dérivé de « ce qui est caché sous les cartes d'un jeu » à « ce qui est caché sous les cartes géographiques ».

      • [^] # Re: Marque

        Posté par  . Évalué à 1.

        t'es pas très cartomancien ;)

        Matricule 23415

      • [^] # Re: Marque

        Posté par  (site web personnel) . Évalué à 1.

        bah… Tu sais, l'histoire et le droit des marques, ça ne marche pas toujours bien ensemble.
        (exemple avec Visa).

  • # Utilise QGIS, c'est fait pour ça !

    Posté par  (site web personnel) . Évalué à 10.

    Hello,

    J'imagine que les personnes dont c'est le métier sont mieux outillés et plus expérimenté et peuvent réaliser cette tâche de manière relativement triviale. Cela m'a mis en appétit pour m'intéresser plus longuement à OSM et aux logiciels de GIS. Sûrement pour le prochain épisode !

    Quand je vois l'effort que tu as fait pour produire cette carte et tes conclusions:, je me dis qu'il faut que je fasse un prochain journal démo sur QGIS !
    Toute la partie de "rendu" aurait pu être créée en quelques clics…

    Dans tous les cas, bravo, tu as bien travaillé, notamment sur la phase acquisition de données qui est clairement le point noir pour débuter: elle implique de récupérer des données réparties chez plusieurs fournisseurs, d'y appliquer des traitements pour obtenir ce qu'on veut et enfin, de convertir tout ça dans un format lisible.

    En attendant, je te conseille de lire deux documents de référence:

    Ces documents, qui sont très majoritairement traduits en français, te permettront de te familiariser avec les concepts des SIG que tu as déjà abordés en partie tout seul, ainsi que de prendre en main QGIS, le logiciel de SIG libre de référence (enfin, c'est mon point de vue assumé).

  • # Juste merci

    Posté par  . Évalué à 10.

    meme si le resultat est perfectible (station service dans la mer)
    l'idée et le cheminement est là.

    et ca fait du bien de voir des gens encore capable de decomposer une problematique, pour en sortir un solution.

    Bravo.

  • # moins geek mais fonctionne correctement

    Posté par  . Évalué à 3.

    bonjour

    comme je suis flemmard, j'ai besoin aussi d'une représentation géographique avec des données style big data, avec des couleurs, des icones personnalisé etc …

    bon c'est pas libre, mais ca fait sont boulot, ca prend en entré un fichier csv ou autre et le représente. Sans l'elegance d'une conversion gps/lambert-93 je le conçois

    http://www.gpsvisualizer.com/

    c'est un peu dans le style logiciel libre, documentation incomplète, obligé de reflechir un peu pour une utilisation avec coloration automatique du bleu au rouge, et pour la création de menu contextuel. Habituer au page man incomplete de 1993 ? c'est pour vous.

    sinon pour ta carte, regarde du coté de pain au chocolat/chocolatine qui donne de bon exemple de jolie carte pour plein de données.

    http://blog.adrienvh.fr/2012/10/16/cartographie-des-resultats-de-chocolatine-ou-pain-au-chocolat/

  • # 37

    Posté par  . Évalué à 1.

    Quoi qu’il en soit dans le 37 c’est deux, pas quatre…

  • # En Python, Pandas et Matplotlib

    Posté par  . Évalué à 6.

    J'ai joué un peu avec pandas et matplotlib et voilà ce que ça pourrait donner:

    Gazole

    Le reste est ici : https://gist.github.com/saimn/70edc8916c205dd3fbd8
    avec un notebook pour générer tout ça (s'il y a des amateurs pour améliorer le rendu).
    Pour les couleurs j'ai utilisé une barre de couleur centrée sur la moyenne limité à 5 sigma, mais on voit sur l'histogramme que c'est pas vraiment gaussien.

  • # Un GIS à la rescousse

    Posté par  . Évalué à 4.

    Suivant vos recommandations, je me suis intéressé à QGIS.
    Mon expérience ne mérite pas encore un journal mais au moins un petit commentaire pour vous montrer ce qu'un néophyte peut obtenir en quelques minutes.
    Le fond de carte est obtenu avec les données de l'IGN (merci SlowBrain pour le lien).
    Les données CSV que j'avais préparées s'importent directement.
    Ensuite il faut jouer un peu avec les styles.

    projet QGIS

    Vraiment bien foutu ce logiciel.

Suivre le flux des commentaires

Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.