lucidiot's cybrecluster

Heatmaps de réseaux Wi-Fi, épisode 4

Lucidiot — Informatique — 2021-11-30
Le sprint final, ou plutôt le marathon final, vers la tant désirée heatmap.

---

Nous revoici pour la suite et la fin de l'incroyable épopée de génération de heatmaps en SQL. Si vous avez suivi tous les articles jusque ici sans carboniser l'intégralité vos neurones, bravo. Je vais aujourd'hui m'assurer que le peu de neurones qu'il vous reste soient exterminés.

Précédemment, nous avons découpé la carte en pixels, calculé une densité de réseaux, téléchargé un fond de carte depuis OpenStreetMap et généré une image représentant cette densité. Nous allons maintenant tout assembler pour former une seule image, notre fameuse heatmap.

Construire un PNG transparent

Je comptais à la base envoyer directement le fichier PGM à ImageMagick pour pouvoir en faire le traitement. Cependant, ImageMagick semble avoir du mal avec ce format et je n'ai jamais réussi à faire une image cohérente avec. À la place, on va faire appel à un utilitaire de la librarie netpbm, dont le seul but est de traiter ce format de fichier. pnmtopng a un nom relativement explicite et peut convertir n'importe quel fichier de la famille PNM vers du PNG.

Pour rappel, notre image PGM ressemble à ceci :

Il serait assez étrange d'utiliser ça directement pour construire l'image ; dans une heatmap, il y a le fond de carte, et par dessus une couleur qui devient plus ou moins transparente ou change de teinte en fonction de l'intensité de la valeur à représenter. Il existe cela dit un argument à pnmtopng qui permet d'exploiter les données de cette image comme de la transparence : -alpha.

Avec -alpha, on utilise notre image comme étant un masque par dessus une autre image. C'est comme un pochoir : plus le pixel dans le masque se rapproche du blanc, moins le pixel de l'image qui est convertie sera transparent. Un pixel blanc du masque donnera donc un pixel entièrement opaque, et un pixel noir donnera un pixel entièrement transparent.

Pour pouvoir utiliser cette fonctionnalité, il va nous falloir du coup avoir une autre image qui passera à travers le masque. On peut pour cela utiliser un autre utilitaire qui va créer une image pour nous : ppmmake. Cet outil génère des fichiers Portable Pixmap, c'est-à-dire des images en couleur, donc on va pouvoir en profiter pour quitter le domaine du monochrome et utiliser une couleur.

J'ai choisi le bleu, qui semblait être assez visible sur le fond de carte d'OpenStreetMap. ppmmake demande la couleur, puis la largeur et la hauteur de l'image, et il renverra une image entièrement bleue. On utiliserait donc une ligne de shell comme celle-ci :

ppmmake blue $width $height | pnmtopng -alpha $alphamask > image.png

Afin de nous épargner pour l'instant toute l'intégration en SQL, je vais utiliser des variables pour certains endroits où la requête SQL se chargera pour nous de remplir les champs.

On obtient maintenant ce genre d'image :

Assembler des calques

Nous disposons donc maintenant de deux images : le fond de carte, et la densité des réseaux. On va donc les superposer, comme on superposerait des calques dans un logiciel de traitement d'image classique. Pour pouvoir le faire de façon automatique dans un script shell, on va utiliser ImageMagick.

Si vous ne connaissez pas ImageMagick, c'est un outil incroyablement puissant qui est utilisable dans beaucoup de langages mais fournit surtout des outils en ligne de commande comme convert, mogrify ou identify qui permettent un grand nombre d'opérations sur des images. Il est possible de s'en servir pour des GIF animés ou des effets 3D, ou pour des opérations toutes simples comme des redimensionnements. Mes deux outils préférés pour l'édition d'image sous Linux sont ImageMagick et Inkscape. Oui, j'utilise un éditeur d'images vectorielles pour modifier des images matricielles histoire de bien faire peur aux graphistes. Mais revenons à nos moutons.

ImageMagick a une documentation très complète, bien qu'il soit souvent difficile de naviguer dedans en raison de la complexité de certains des sujets abordés et surtout de la quantité de fonctionnalités disponibles. En général, la documentation la plus fréquemment consultée est plutôt StackOverflow.

On va utiliser convert, qui est conçu pour prendre zéro, une ou plusieurs images en entrée et donner une ou plusieurs images en sortie. Ici, il aura deux images en entrées et une seule en sortie.

La superposition classique comme celle qu'on veut peut se faire assez facilement en utilisant -composite. Les outils en ligne de commande d'ImageMagick ne travaillent pas du tout comme des commandes classiques POSIX et les arguments sont traités dans l'ordre, comme des pipes d'un shell. Ici, -composite va s'attendre à ce que tous les arguments précédents soient des images, et les superposera ensemble pour donner une seule image.

En considérant que nous avons au préalable enregistré le fond de carte dans un fichier PNG comme nous l'avons fait pour le fichier PGM, nous pouvons rajouter convert en bout de chaîne de notre ligne précédente :

ppmmake blue $width $height |
  pnmtopng -alpha $alphamask |
  convert $map - -composite - > image.png

On dit d'abord à convert de prendre notre carte ainsi que l'entrée standard, c'est-à-dire nos données de réseaux Wi-Fi, de superposer les images, et de renvoyer le résultat sur la sortie standard.

On obtient quelque chose comme ceci :

C'est bon, on a une heatmap, c'est gagné, on remballe tout... Quoi, comment ça cette carte fait n'importe quoi ?

Redimensionner le PGM

Pour pouvoir redimensionner le calque à la carte OpenStreetMap, on peut utiliser -resize. C'est un autre opérateur d'ImageMagick qui permet de redimensionner une image selon des dimensions données. Cependant, -composite s'attend à avoir vraiment juste des images, pas des opérations précédant la composition.

ppmmake blue $width $height |
  pnmtopng -alpha $alphamask |
  convert - -resize 990x1131 layer.png
convert $map layer.png -composite image.png

On utilise cette fois deux commandes convert. La première va seulement faire le redimensionnement puis stocker une image PNG qui contiendra seulement notre calque redimensionné. La seconde prend les deux images et les superpose. Cela pourrait suffire, sauf que les outils d'ImageMagick ont vraiment été conçus pour permettre de faire énormément de traitements avec une seule image. Il serait dommage de ne pas en profiter.

ppmmake blue $width $height |
  pnmtopng -alpha $alphamask |
  convert $map \( \
    - \
    -resize 990x1131 \
  \) -composite image.png

On a ici placé la première commande entre parenthèses dans la deuxième. On utilise des \( pour garantir que le shell ne va pas nous jouer des tours et qu'on envoie bien les parenthèses à ImageMagick directement. Cette commande unique évite d'utiliser un fichier temporaire ; on récupère la carte, on construit l'image redimensionnée depuis l'entrée standard, et on assemble le tout. Avec ça, on donne l'illusion à -composite qu'il y a bien juste deux images.

Détecter les dimensions de la carte

Le rendu d'OpenStreetMap est non seulement bien plus grand que le PGM, mais en plus sa taille peut varier d'un export à l'autre. Un nouveau nettoyage de la base de données changeant ses statistiques d'étendue des données géographiques, ou une différence dans la façon dont OpenStreetMap arrondit les coordonnées peuvent nous conduire à observer des différences de quelques pixels sur plusieurs exécutions. Il n'est donc pas possible comme je viens de le faire d'utiliser directement des dimensions dans le script.

Puisqu'on ne pourra pas prédire même dans la requête SQL la taille finale de la carte, on peut à la place utiliser un autre outil d'ImageMagick pour obtenir les dimensions de l'image automatiquement. identify est l'outil qu'il nous faut.

Avec identify, on passe en paramètre des noms d'image et on récupère en sortie toutes sortes d'informations sur l'image. On peut utiliser -verbose pour que l'outil vomisse une quantité incroyable d'informations diverses et variées sur l'image, mais ce qui nous intéresse plutôt ici serait de récupérer seulement la largeur et la hauteur.

L'argument -format va nous permettre de construire une valeur qui ressemble à un argument valide pour l'opérateur -resize de convert :

ppmmake blue $width $height |
  pnmtopng -alpha $alphamask |
  convert $map \( \
    - \
    -resize $(identify -format '%wx%h' $map) \
  \) -composite image.png

-resize s'attend à des dimensions du type HxL, ou hauteur fois largeur. L'argument -format permet d'utiliser une chaîne de caractères particulière où on peut placer des symboles particuliers qui seront traduits en valeurs venus de l'image : ici %w sera remplacé par la largeur et %h par la hauteur. On récupère donc une valeur semblable à 990x1131, et convert l'acceptera sans problème.

On fait du progrès ! Mais ce n'est pas encore ça. Si vous regardez attentivement, vous remarquerez que tout est encore décalé : il y a des tracés liés à des rues qui se retrouvent au milieu d'un cours d'eau et les cercles ont l'air écrasés. Et un autre petit détail, pas très difficile à corriger, est qu'on ne voit pas bien la résolution de la carte : le redimensionnement ajoute un flou qui pourrait être esthétiquement plaisant mais qui ne m'aide pas vraiment à me représenter la résolution de ma heatmap.

Éviter le flou

Pour éviter le flou, qui est surtout peu pratique quand on réduit la résolution de la carte (donc qu'on a de très grands carrés) pour la générer plus vite, il suffit d'utiliser le paramètre -filter, qui permet de spécifier un filtre de rééchantillonage. Les filtres de rééchantillonage servent à déterminer la couleur finale de chaque pixel : des opérations comme le redimensionnement peuvent faire que des pixels originaux se retrouvent à cheval entre plusieurs pixels de destination, ou inversement qu'un seul pixel de destination soit maintenant composé de plusieurs pixels source. On appelle donc filtre de rééchantillonage un algorithme particulier qui permet de calculer la couleur d'un pixel de destination.

Il existe moult filtres et aucun n'est parfait, tout dépend des opérations effectuées, du style de résultat souhaité, de la vitesse d'exécution souhaitée, etc. La documentation sur les filtres de rééchantillonage d'ImageMagick est très longue et ne donne pas vraiment de conclusion sur quel est le meilleur filtre, juste les avantages et les inconvénients de chaque. Dans notre cas, on va cependant pouvoir utiliser le filtre Point, qui est probablement le plus simple : le pixel source le plus proche est utilisé pour le pixel de destination. Cela va directement agrandir le pixel.

convert $map \( \
  - \
  -filter Point \
  -resize $(identify -format '%wx%h' $map) \
\) -composite image.png

Il faut bien faire attention à configurer le filtre avant de déclencher le redimensionnement, sinon le paramètre serait ignoré.

Déformer en un rectangle

J'avais rapidement mentionné au tout début de cette épopée que la grille utilisée pour grouper les réseaux Wi-Fi par cases n'était pas forcément une grille de carrés. En fait, les cases de la grille sont mesurés en degrés, donc leur forme dépend de la projection de carte utilisée et des coordonnées de la case. Par exemple, au niveau des pôles, les carrés seraient en fait des triangles. Au niveau de l'équateur, les cases seront bien carrées, et elles se déformeront progressivement en rectangles, trapèzes puis triangles en se rapprochant des pôles.

Je vais considérer ici que les trapèzes et les triangles n'existent pas puisque la grille en Europe de l'Ouest n'a que des rectangles. Techniquement ce sont déjà des trapèzes en raison de la courbure de la Terre, mais l'angle est tellement faible qu'on ne le constate pas encore. Travailler uniquement avec des angles droits va simplifier beaucoup le travail.

L'opérateur de redimensionnement -resize n'est capable que de faire des redimensionnements simples où l'image n'est pas déformée : le facteur de redimensionnement est le même pour la hauteur que pour la largeur. Si on demande un redimensionnement avec des dimensions farfelues, par exemple passer d'une image de 32×32 à 420×69, l'image obtenue sera incluse dans le rectangle mais pas déformée, donc on aura ici une image de 69×69 pixels. Pour pouvoir faire une déformation, il va falloir utiliser l'opérateur de… déformation : -distort.

L'opérateur de déformation ou distorsion est un opérateur extrêmement générique qui permet de tordre une image dans tous les sens et avec beaucoup d'options. On peut par exemple créer une représentation en carte polaire à partir d'un planisphère, s'adapter aux paramètres particuliers d'une lentille d'appareil photo qui déformerait les bords de l'image, ou manipuler des images sur trois dimensions. Oui, on peut faire du rendu 3D, avec exactement les mêmes techniques qui rendent possible de créer des jeux en trois dimensions uniquement en CSS sans avoir recours à des technologies comme WebGL.

Nous allons utiliser la méthode de distorsion SRT, ou Scale, Rotate, Translate. C'est la distorsion la plus simple disponible, et elle simule les mêmes opérations que -resize, -rotate, -scale, etc. avec plus de flexibilité. C'est en fait un alias qui va générer une matrice de transformation affine. Ce genre de matrice est assez complexe à maîtriser quand vous n'y connaissez pas grand chose, et pour des opérations relativement simples de déformations, l'alias SRT est bien plus simple. Cette distorsion regroupe en fait trois distorsions séparées, appliquées l'une après l'autre : on met à l'échelle (autrement dit, on étire ou rétrécit) l'image, on la tourne, et on la déplace vers une position donnée.

Les paramètres de la déformation SRT sont assez inhabituels vu que beaucoup d'entre eux sont optionnels, par exemple pour permettre de seulement tourner l'image, seulement la déplacer, etc. Nous souhaitons étirer l'image différemment sur ses deux axes pour qu'elle rentre bien dans le cadre, donc nous avons besoin d'une syntaxe qui permet de donner deux facteurs d'échelle. La documentation de SRT nous donne ceci :

-distort SRT "X,Y ScaleX,ScaleY Angle"

Il n'est pas possible d'avoir la translation (X,Y) et la mise à l'échelle selon l'axe (ScaleX,ScaleY) toutes deux optionnelles ici puisqu'on ne saurait plus si 1,1 correspond à la mise à l'échelle ou la translation, donc on va se retrouver à devoir juste indiquer 0,0 en translation. On mettra aussi 0 en angle, puisque l'angle est toujours obligatoire. Il nous reste donc à déterminer ScaleX et ScaleY. Pour cela, on va plonger encore un peu plus dans la syntaxe étrange des formats de ImageMagick.

Calculer les facteurs d'échelle

ImageMagick s'attend à un facteur d'échelle, c'est à dire un multiplicateur pour les dimensions de l'image : 1 correspond à aucune déformation, 2 deux fois plus grand, 0.5 deux fois plus petit.

Nous avons deux images dont les dimensions sont inconnues. On veut redimensionner une des deux pour qu'elle corresponde à l'autre. En considérant que les dimensions de l'image d'origine sont x₁ et y₁ et celles de l'image redimensionnée x₂ et y₂, on peut obtenir les facteurs d'échelle à appliquer en calculant x₁ / x₂ et y₁ / y₂.

On a vu précédemment qu'on pouvait utiliser %w et %h dans les options de format de identify pour obtenir la largeur et la hauteur de l'image. Les échappements de caractères particuliers de ImageMagick ne sont pas utilisables seulement dans les options de format de sortie mais en fait dans quasiment tous les arguments d'ImageMagick. On peut donc directement s'en servir dans notre distorsion SRT. Mais comment faire un calcul sur les dimensions de deux images ?

On va aller pour cela un peu plus loin dans cet échappement de caractères. Les paramètres d'une seule lettre comme %w s'épuisent rapidement, donc des paramètres aux noms plus longs sont utilisables avec des crochets, comme par exemple %[Orientation] pour une orientation déterminée par des métadonnées EXIF. Il existe aussi des préfixes, donc on peut aussi utiliser %[exif:Orientation]. Un préfixe particulier va nous être utile ici : %[fx:…], qui permet d'utiliser des expressions FX.

Le langage FX a été initialement conçu pour l'opérateur -fx (Special Effects). Le langage FX dispose de conditions, de boucles, de variables, etc., et est Turing-complet, c'est-à-dire qu'il est théoriquement possible de programmer n'importe quoi avec. La puissance de ce langage fait qu'il a été étendu pour être mis à disposition dans un grand nombre d'opérateurs d'ImageMagick, et donc on peut s'en servir ici dans les chaînes de caractères formatées. On va s'en servir ici pour faire une simple division : quelque chose comme %[fx:1/2] par exemple renvoie 0.5.

Nous allons maintenant revenir à notre bon vieux identify qu'on a utilisé pour le redimensionnement :

identify -format '%wx%h' $map

On peut transformer nos deux paramètres en des expressions FX :

identify -format '%[fx:w]x%[fx:h]' $map

On n'utilise pas de % ici car w et h nous sont déjà mis à disposition comme des variables. Cette syntaxe forcera ImageMagick à interpréter chaque expression FX, trouver la variable, et simplement en sortir la valeur. On pourrait aussi écrire ça en imbriquant du formatage dans du formatage :

identify -format '%[fx:%w]x%[fx:%h]' $map

Dans ce cas-ci, ImageMagick effectuera deux passes sur la chaîne de caractères : d'abord pour transformer %w et %h en leurs nombres respectifs, ce qui donnerait quelque chose comme '%[fx:990]x%[fx:1131]', puis pour évaluer les expressions FX, qui sont cette fois-ci des constantes. On obtient alors 990x1131.

Plaçons cet identify de sorte à l'utiliser directement comme argument de notre distortion SRT :

-distort SRT "$(identify -format '0,0 %[fx:%w],%[fx:%h] 0' $map)"

Nous allons générer ici 0,0 990,1131 0, une distortion valide mais qui redimensionnerait beaucoup trop fort l'image vu qu'on l'agrandirait 990 fois en largeur et 1131 fois en hauteur. Pour pouvoir faire notre ratio, nous allons profiter du fait que FX soit disponible à la fois dans -format et dans -distort : On va générer, et pas évaluer, une expression FX dans -format, et elle sera évaluée ensuite dans -distort. Puisque notre identify travaille sur la carte, il pourra fournir les dimensions de la carte, et puisque -distort travaille sur le calque transparent, on disposera des dimensions de ce calque.

-distort SRT "$(identify -format '0,0 %%[fx:%w/w],%%[fx:%h/h] 0' $map)"

On utilise ici des %%, qui sont la syntaxe pour mettre littéralement un % sans qu'il soit interprété comme une expression de formatage. Ainsi, seuls %w et %h seront véritablement interprétés, et on va cette fois générer une expression de ce type :

-distort SRT "0,0 %[fx:990/w],%[fx:1131/h] 0"

Cette expression sera alors réévaluée, et w et h seront les dimensions du calque. Nous avons donc pu effectuer notre division et avons donc un facteur d'échelle, en imbriquant une chaîne formatée d'ImageMagick dans une chaîne formatée d'ImageMagick… On n'a plus qu'à intégrer le tout à l'expression convert originale :

convert $map \( \
  - \
  -filter Point \
  -distort SRT \
    "$(identify -format '0,0 %%[fx:%w/w],%%[fx:%h/h] 0' $map)" \
\) -composite -

On obtient maintenant une image redimensionnée correctement… mais on a créé deux nouveaux problèmes.

Éviter la perte de données

On a de nouveau une image quasiment vide ! On peut observer quelques données en haut à droite, mais la plupart des pixels bleus sont absents. Une bonne lecture de la documentation de l'opérateur de distorsion indique que la distorsion va essayer de rogner automatiquement l'image déformée pour correspondre à la taille de l'image d'origine, ce qui peut être une action utile dans certains cas mais ne nous aide pas du tout ici. La mise à l'échelle va du coup faire complètement sortir l'image de son cadre.

Certains opérateurs d'ImageMagick disposent de variantes, qui soit travaillent juste différement soit de façon totalement inverse. Les variantes sont disponibles en remplaçant le - par un +. L'opérateur +distort fonctionne exactement comme une distorsion normale, mais n'applique pas ce rognage automatique. À la place, il utilise le canevas virtuel, un paramètre particulier qu'on retrouve dans certains formats d'image.

Le canevas virtuel permet de rogner une image pour que seule une portion des pixels stockés dans le fichier soient véritablement affichés : une autre technique étrange si vous voulez faire du rognage très rapide sans recalculer tous les pixels. On peut aussi utiliser cette technique pour travailler avec des coordonnées négatives. À la différence d'un rognage classique, le rognage n'affecte pas du tout les pixels réellement stockés, donc on peut juste annuler ce rognage complètement pour le retirer et retrouver notre image entière.

On peut utiliser l'opérateur -repage pour redéfinir le canevas virtuel, et sa variante +repage pour supprimer ce canevas, c'est donc ce qu'on va faire après la distorsion :

convert $map \( \
  - \
  -filter Point \
  +distort SRT \
    "$(identify -format '0,0 %%[fx:%w/w],%%[fx:%h/h] 0' $map)" \
  +repage \
\) -composite -

Le canevas virtuel n'est pas seulement qu'une hauteur et une largeur mais aussi des coordonnées X et Y, puisqu'un rognage se fait avec une paire de coordonnées. L'opérateur +distort prend cette notion en compte, et si par exemple nous avions fait une translation sur notre image pour la déplacer vers le haut et la gauche, des pixels se seraient retrouvés en coordonnées négatives, et nous aurions pu les récupérer avec ce système de canevas ou au moins positionner l'image correctement par rapport à la carte, sans que le point d'origine n'aie changé. Mais ici, la mise à l'échelle ne s'est effectuée que sur des coordonnées positives, donc on peut juste détruire le canevas virtuel sans réfléchir.

Éviter à nouveau le flou

On a retrouvé nos données, mais cette fois-ci le flou est de retour. C'est assez beau, plus doux que des rectangles colorés, mais le flou est trop fort : les petits rectangles sur certaines rues avec des valeurs peu élevées sont quasiment invisibles. Pourtant, on avait défini un filtre de rééchantillonage pour retirer le flou !

Il va nous falloir maintenant mentionner l'interpolation. C'est un peu le même principe que le filtre de rééchantillonage mentionné plus haut. L'interpolation est une version plus avancée du rééchantillonage qui s'applique pour des cas particuliers de la distorsion, où on peut se retrouver avec des pixels dont les coordonnées se retrouvent à la frontière entre deux pixels source. Avec notre distortion faisant une mise à l'échelle, beaucoup de points vont se retrouver dans cette situation, et l'interpolation entre donc en jeu : par défaut, c'est une interpolation bilinéaire qui est utilisée. C'est une méthode assez classique qui est un bon compromis entre performances et qualité du flou.

On peut définir nous-mêmes la méthode d'interpolation à utiliser pour dire à ImageMagick comment inférer la couleur d'un pixel à partir des pixels source situés à proximité. La plupart des méthodes vont conduire à des flous divers et variés, mais la méthode nearest neighbor (voisin le plus proche) va juste chercher le ou les points les plus proches et prendre la couleur de là, sans faire de calculs qui conduiraient à un flou. C'est donc exactement le même comportement que notre filtre de rééchantillonage Point ! On peut le définir en ajoutant -interpolate Nearest tout comme on a mis -filter Point.

ppmmake blue $width $height |
  pnmtopng -alpha $alphamask |
  convert $map \( \
    - \
    -filter Point \
    -interpolate Nearest \
    +distort SRT \
      "$(identify -format '0,0 %%[fx:%w/w],%%[fx:%h/h] 0' $map)" \
    +repage \
  \) -composite -

C'est bon, le mal de tête des manipulations d'image est terminé, et nous obtenons maintenant la heatmap tant attendue :

Un script shell générable en SQL

Il reste encore quelques efforts : il faut maintenant faire le lien entre la requête SQL qui génère un fichier PNM et notre génération d'image, pour finalement arriver à un script que je peux lancer en une seule ligne pour générer une heatmap avec toutes les données disponibles sur la base de données. Commençons par convertir le script shell en quelque chose qui pourra être utilisé au sein même de la requête SQL :

#!/bin/sh -e
alphafile="$(mktemp)"
mapfile="$(mktemp)"
trap "rm \"$alphafile\" \"$mapfile\"" EXIT
curl -s --fail \
  '[URL de rendu OSM ici]' \
  -H "Cookie: _osm_totp_token=$(
    curl -s --fail -I https://www.openstreetmap.org/ |
    sed -En '/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p'
  )" \
  -o "$mapfile"
printf "P2\n [WIDTH] [HEIGHT]\n65535\n[PIXELS]\n" > "$alphafile"
ppmmake blue [WIDTH] [HEIGHT] \
  | pnmtopng -alpha "$alphafile" \
  | convert "$mapfile" \( \
      - \
      -filter Point \
      -interpolate Nearest \
      +distort SRT \
        "$(identify -format '0,0 %%[fx:%w/w],%%[fx:%h/h] 0' $mapfile)" \
      +repage \
  \) -composite -

1. On définit deux fichiers temporaires : un pour la carte OpenStreetMap (mapfile) et un pour les données PNM (alphafile).

2. On s'assure que ces fichiers soient supprimés automatiquement quand le script s'arrête, même s'il plante à cause d'une erreur, avec la commande trap.

3. On utilise la ligne de commande curl qu'on avait construit dans l'épisode 2 pour télécharger le fond de carte et le stocker dans $mapfile.

4. On stocke les données PNM dans le fichier $alphafile. On utilise printf puisque cela permettra à la requête SQL de remplir les trous avec les vraies données.

5. On exécute tout notre bazar de construction d'image !

Requête SQL finale

On peut maintenant intégrer ce script shell au sein de la requête SQL. Au lieu de renvoyer un fichier PNM comme la dernière fois, la requête va maintenant renvoyer le script directement, rempli avec toutes les bonnes valeurs :

WITH output (i, j, geom, value) AS (
  SELECT i, j, geom, COUNT(geo)
  FROM ST_SquareGrid(0.0001, ST_SetSRID(ST_EstimatedExtent('geometries', 'geo'), 4326)) AS grid
  LEFT JOIN geometries ON ST_Intersects(grid.geom, geo)
  GROUP BY i, j, geom
), stats (width, height, ratio, bbox) AS (
  SELECT
    MAX(i) - MIN(i) + 1,
    MAX(j) - MIN(j) + 1,
    MAX(value) / 65535.0,
    ST_SetSRID(ST_Extent(geom), 4326)
  FROM output
), lines (line) AS (
  SELECT string_agg(CEIL(value / ratio)::text, ' ' ORDER BY i)
  FROM output, stats
  GROUP BY j
  ORDER BY j DESC
)
SELECT '#!/bin/sh -e
alphafile="$(mktemp)"
mapfile="$(mktemp)"
trap "rm \"$alphafile\" \"$mapfile\"" EXIT
curl -s --fail \
  ''' || osm_render_url(bbox) || ''' \
  -H "Cookie: _osm_totp_token=$(
    curl -s --fail -I https://www.openstreetmap.org/ |
    sed -En ''/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p''
  )" \
  -o "$mapfile"
printf "P2\n'
|| width || ' ' || height
|| '\n65535\n' || string_agg(line, '\n' ORDER BY j DESC)
|| '\n" > "$alphafile"
ppmmake blue ' || width || ' ' || height || ' \
  | pnmtopng -alpha "$alphafile" \
  | convert "$mapfile" \( \
      - \
      -filter Point \
      -interpolate Nearest \
      +distort SRT \
        "$(identify -format ''0,0 %%[fx:%w/w],%%[fx:%h/h] 0'' $mapfile)" \
      +repage \
  \) -composite -'
FROM lines, stats
GROUP BY width, height, bbox

Script shell final

On peut enfin appliquer une belle dose d'échappement de caractères, et envelopper notre requête SQL au sein d'une commande COPY dans une ligne psql. En toute fin de ce script, qui est techniquement un one-liner, on exécute le script shell renvoyé par la requête SQL et on stocke le résultat dans un fichier PNG.

#!/bin/sh -e
psql -c 'COPY (
WITH output (i, j, geom, value) AS (
  SELECT i, j, geom, COUNT(geo)
  FROM ST_SquareGrid(0.0001, ST_SetSRID(ST_EstimatedExtent('"'geometries'"', '"'geo'"'), 4326)) AS grid
  LEFT JOIN geometries ON ST_Intersects(grid.geom, geo)
  GROUP BY i, j, geom
), stats (width, height, ratio, bbox) AS (
  SELECT
    MAX(i) - MIN(i) + 1,
    MAX(j) - MIN(j) + 1,
    MAX(value) / 65535.0,
    ST_SetSRID(ST_Extent(geom), 4326)
  FROM output
), lines (line) AS (
  SELECT string_agg(CEIL(value / ratio)::text, '"' '"' ORDER BY i)
  FROM output, stats
  GROUP BY j
  ORDER BY j DESC
)
SELECT '"'"'#!/bin/sh -e
alphafile="$(mktemp)"
mapfile="$(mktemp)"
trap "rm \"$alphafile\" \"$mapfile\"" EXIT
curl -s --fail \
  '"'''"' || osm_render_url(bbox) || '"'''"' \
  -H "Cookie: _osm_totp_token=$(
    curl -s --fail -I https://www.openstreetmap.org/ |
    sed -En '"''"'/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p'"''"'
  )" \
  -o "$mapfile"
printf "P2\n'"'"'
|| width || '"' '"' || height
|| '"'"'\n65535\n'"'"' || string_agg(line, '"'"'\n'"'"' ORDER BY j DESC)
|| '"'"'\n" > "$alphafile"
ppmmake blue '"'"' || width || '"' '"' || height || '"'"' \
  | pnmtopng -alpha "$alphafile" \
  | convert "$mapfile" \( \
      - \
      -filter Point \
      -interpolate Nearest \
      +distort SRT \
        "$(identify -format '"''"'0,0 %%[fx:%w/w],%%[fx:%h/h] 0'"''"' $mapfile)" \
      +repage \
  \) -composite -'"'"'
FROM lines, stats
GROUP BY width, height, bbox
) TO STDOUT WITH (FORMAT CSV, QUOTE E'"'"'\t'"'"')' "$@" |
  sh -e > image.png

Il n'y a plus qu'à lancer le script avec les mêmes paramètres de connexion que pour un client psql classique :

./heatmap.sh -h hostname -p 1234 -d database -U user

Et la magie opère, générant une heatmap exactement comme celle montrée au-dessus.

Conclusion

On a donc imbriqué une chaîne de caractères formatée de ImageMagick dans une chaîne de caractères formatée de ImageMagick dans une commande convert dans un script shell POSIX dans une requête SELECT dans une requête COPY dans une commande psql dans un script shell POSIX. De vraies poupées russes informatiques, même si les poupées russes sont quand même bien plus agréables à regarder que ce… tas de texte.

Je ne pensais pas que j'allais me retrouver à écrire plus de 11 000 mots sur cette seule requête SQL, et c'est très probablement mon record de mots par ligne pour une requête dans tout ce projet de warwalking. Le développement de cette requête a pris plusieurs soirées, sans compter celles passées à déchiffrer toute la syntaxe pour expliquer progressivement la requête dans des articles : quand j'écris la requête SQL, je vais un peu tester des fonctionnalités sans trop réfléchir à comment elles fonctionnent, mais dans des articles, il faut bien que je dise de quoi il s'agit.

J'ai encore quelques pistes que je pourrais explorer pour améliorer cette heatmap ou en générer d'autres, et sur quelques autres sujets liés au warwalking, mais on verra ça à d'autres occasions. Pour l'instant, on peut considérer que la saga s'achève ici !

---

Commentaires

---

Soumettre