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

Lucidiot — Informatique — 2021-11-07
La première étape de fabrication d'une énième requête SQL incompréhensible.

Non, je ne suis toujours pas à court de choses à écrire concernant mon warwalking. Je voulais commencer à analyser de façon plus générale mes données, et une façon de faire ça sans voir une montagne de points incompréhensible est d'utiliser une heatmap. C'est une carte où on va représenter non pas les points de données mais leur densité, et ça peut donner quelque chose qui ressemble à une carte de radar météorologique ou de caméra thermique.

Comme vous l'avez probablement remarqué si vous avez lu les articles précédents, j'aime faire les choses en utilisant des technologies pas du tout adaptées à ça, donc on va essayer de faire le maximum en SQL. Et cette tâche a été suffisamment complexe pour mériter plusieurs articles !

Il existe des librairies notamment en JS capable de générer une heatmap à partir d'un jeu de données et de l'afficher une carte comme celle d'OpenStreetMap. Ces librairies vont souvent utiliser des représentations vectorielles; elles dessinent des contours ou des courbes de niveau qui représentent les couleurs à afficher, et le navigateur se débrouillera pour générer des dégradés de couleurs qui suivent ces courbes. Représenter ce genre de choses en SQL me semble vraiment complexe et j'ai préféré prendre une approche un peu plus simple et plus compatible avec une requête en base de données : découper la carte en pixels, et faire des statistiques sur chacun de ces pixels. C'est ce qu'on va commencer à faire aujourd'hui.

Construire une grille

On va essayer de générer une grille composée de beaucoup de petits carrés, et mesurer le nombre de points dans chaque carré. Pour cela, il nous faudra deux choses : d'abord, pouvoir connaître l'étendue de toutes les données sur la carte, c'est-à-dire récupérer un rectangle qui englobe l'intégralité des données, et ensuite pouvoir découper ce rectangle en carrés.

Pour obtenir une boîte englobante (ou bounding box), on peut utiliser des fonctions d'agrégation de PostGIS, notamment ST_Extent et ST_EstimatedExtent. Cette dernière va utiliser non pas les données elles-mêmes mais les statistiques fournies par PostgreSQL pour la colonne entière, statistiques qui peuvent s'avérer imprécises et ne vont souvent pas englober 100% des données mais souvent un petit peu moins. Ces statistiques peuvent être mises à jour avec VACUUM ANALYZE.

SELECT ST_Extent(geo::geometry) FROM all_map;

-- Ne fonctionne que si `geo` était déjà de type `geometry`
SELECT ST_EstimatedExtent('all_map', 'geo');

J'ai commencé à travailler en utilisant ST_Extent directement car la colonne geo de la vue all_map est de type geography (mieux adapté pour l'affichage automatique sur une carte que fournit pgAdmin). Utiliser ST_EstimatedExtent renvoie alors une erreur. On reviendra plus tard sur l'optimisation de toute façon.

On obtient avec ces deux méthodes un type box2d, qui n'a aucune notion d'un SRID et qu'on ferait mieux de convertir maintenant pour éviter d'avoir le problème des œufs de Pâques qu'on a eu lors de la dernière tentative de cartographie. Je vous laisse l'honneur de lire les articles précédents si je vous ai déjà perdu avec les SRID ou les œufs…

SELECT ST_SetSRID(ST_Extent(geo::geometry)::geometry, 4326)
FROM all_map;

Maintenant qu'on a le rectangle englobant toutes nos données, on va pouvoir découper en plein de petits carrés.

Réinventons la roue

S'il me fallait implémenter une génération de grille à partir de deux points et d'une dimension donnée pour les carrés en Python, je m'y prendrais probablement d'une manière similaire à celle-ci :

def makegrid(x1, y1, x2, y2, size):
   for x in range(x1, x2, size):
      for y in range(y1, y2, size):
         yield (x, y, x + size, y + size)

Pas très compliqué ; on utilise deux fois range pour pouvoir boucler sur toutes les coordonnées des droites qui vont trancher le rectangle, et on renvoie deux paires de coordonnées représentant le carré. C'est l'approche que j'ai pris quand j'ai commencé à essayer d'implémenter une grille en SQL :

CREATE OR REPLACE FUNCTION makegrid(geometry, numeric)
RETURNS setof geometry AS
'SELECT ST_MakeEnvelope(x-$2, y-$2, x, y, ST_SRID($1)) AS square FROM
generate_series(ST_XMin($1)::numeric, ST_XMax($1)::numeric, $2) AS x,
generate_series(ST_YMin($1)::numeric, ST_YMax($1)::numeric, $2) AS y'
LANGUAGE sql;

J'utilise ST_XMin, ST_XMax, ST_YMin et ST_YMax pour accéder aux quatre coordonnées de la boîte. generate_series est l'équivalent pour PostgreSQL de range; c'est une fonction qui renvoie une table d'une seule colonne où chaque ligne est un chiffre. En faisant une jointure sans aucune condition (un FULL OUTER JOIN), je calcule toutes les combinaisons entre les deux séries x et y, donc je fais une boucle imbriquée. J'utilise ensuite ST_MakeEnvelope qui accepte cinq paramètres : deux paires de coordonnées et un SRID que je récupère auprès de la géométrie originale. La fonction renvoie une géométrie qui n'est rien d'autre qu'un rectangle ayant ces deux coordonnées, donc je me retrouve avec une fonction qui renvoie une table d'une seule colonne comprenant des géométries, autrement dit un setof geometry.

Utilisons une roue existante

C'est à ce moment-là que j'ai... laissé un peu tomber l'idée de la heatmap. J'avais commencé à réfléchir à cette heatmap il y a un long moment, alors que le projet de warwalking n'avait pas encore vraiment maturé, et je n'ai pas vraiment été plus loin que cette fonction. Quand je suis revenu dessus un mois ou deux plus tard, j'ai trouvé un peu par accident dans la documentation de PostGIS une fonction qui génère toute seule une grille exactement comme je le souhaite. Je ne sais pas vraiment pourquoi je ne l'avais pas vue avant, du coup j'ai réinventé une roue sans vraiment le faire exprès.

WITH extent (bbox) AS (SELECT ST_SetSRID(ST_Extent(geo::geometry)::geometry, 4326) FROM all_map)
SELECT grid.*
FROM extent, ST_SquareGrid(0.0001, bbox) AS grid;

Cette fonction ST_SquareGrid prend les mêmes paramètres que la mienne, mais dans l'ordre inverse ; la taille des carrés avant le rectangle englobant. Elle renvoie non pas un setof geometry mais un setof record : une façon de dire qu'elle renvoie une table avec plusieurs colonnes. Elle renvoie en fait trois colonnes : geom, i, et j. geom est la géométrie d'une seule case, et i et j identifient la case de façon unique. On peut trier par i et j pour garantir un certain ordre des cases, par exemple ligne par ligne. Ça peut s'avérer assez utile pour la génération d'images, plutôt que d'essayer d'ordonner selon un polygone.

Notons que i et j sont des index où le 0 correspond au point d'origine du système de référence spatiale. En d'autres termes, pour nos données en coordonnées GPS classiques (donc EPSG 4326), le point d'origine se trouvera à 0° Nord et 0° Ouest. On y trouvera une case dont le bord inférieur droit touche l'équateur et le méridien de Greenwich, et dont la hauteur et la largeur sont de 0.0001 degrés. Vu que le découpage se fait en degrés de coordonnées WGS84, et pas par exemple en mètres, les carrés pourront avoir des superficies différentes et ne pas toujours être carrés au fur et à mesure qu'on se rapproche des pôles. Cela ne pose problème que lorsque la projection de carte utilisée n'est pas la même partout ; on se retrouverait alors à devoir convertir des coordonnées.

Compter par case

On peut maintenant s'atteler au comptage du nombre de réseaux Wi-Fi présents dans chaque case de la grille. J'utilise le résultat de notre vue all_map, qui utilise ST_Buffer pour générer des cercles lorsqu'une information de précision est disponible, donc je dois m'attendre à avoir soit des points soit des polygones (vu que les cercles sont approximés en tant que polygones). J'ai l'embarras du choix pour comparer les polygones ; par intersection, contenu, contenu exclusif (sans toucher les bords), etc.

La fonction qui est la plus rapide et la plus compatible avec mes besoins est ST_Intersects. Il est possible qu'avec un cercle suffisamment grand, il se retrouve à être sur plusieurs cases en même temps, mais ça n'est pas très gênant pour une heatmap. Avec une autre fonction comme ST_Contains, où la géométrie devrait être intégralement à l'intérieur de la case de la grille, je risquerais de perdre un bon nombre de réseaux qui se trouveraient accidentellement en travers des frontières de la grille.

La rapidité de ST_Intersects vient du fait que les polygones ont quelques octets supplémentaires de mémoire qui contiennent les coordonnées d'un rectangle englobant, destiné exactement pour accélérer les comparaisons d'intersection de ce genre. Si le carré autour du cercle ne croise pas du tout celui de ma case, on peut abandonner immédiatement avant de commencer à comparer point par point. De même, si on compare point par point et qu'un seul se trouve dans la case, c'est gagné et on n'a pas besoin de continuer à comparer d'autres points. Pour ST_Contains, il faut vérifier qu'aucun point ne déborde, et bien qu'on puisse imaginer qu'une optimisation avec ce rectangle englobant incorporé serait possible, elle ne semble pas vraiment être faite pour des polygones.

WITH extent (bbox) AS (SELECT ST_SetSRID(ST_Extent(geo::geometry)::geometry, 4326) FROM all_map)
SELECT grid.geom, COUNT(geo)
FROM extent, ST_SquareGrid(0.0001, bbox) AS grid
LEFT JOIN all_map ON ST_Intersects(grid.geom, geo)
GROUP BY grid.geom;

Je collecte la grille dans une table et en fais une nouvelle jointure avec notre table de positions de réseaux. J'utilise un LEFT JOIN ici car il est tout à fait possible pour des cases de grilles de ne contenir aucun réseau, et je veux garder les cases quand même avec un décompte de zéro. La jointure fera que chaque case de la grille sera répétée une fois pour chaque réseau, ou renvoyée une seule fois avec aucun réseau. On peut donc utiliser un COUNT() tout de suite pour récupérer le résultat qu'on veut. On peut utiliser ST_Contains directement comme condition de jointure, puisqu'on cherche les réseaux qui sont contenus dans la case pour chaque case. Des conditions de jointure simples comme celles-ci peuvent aussi souvent être optimisées par PostgreSQL.

Optimisation

Techniquement, cette requête fonctionne. Et si vous avez de l'expérience avec les gens qui utilisent le mot « techniquement, » vous savez que si quelqu'un utilise ce mot, vous pouvez en toute sécurité ignorer le reste de la phrase. J'ai revérifié toutes les requêtes de cet article tout en écrivant, puisque j'écris quasiment toujours mes articles post-mortem et je me retrouve à revenir sur mes pas pour comprendre tout ce que j'ai fait jusque là. J'ai commencé ce paragraphe alors que la requête tournait encore, depuis plus d'une demi-heure ; j'ai finalement abandonné et j'ai regardé le résultat de EXPLAIN pour comprendre un peu mieux la situation.

La vue all_map, et je le savais déjà, n'a pas beaucoup d'opportunités pour utiliser efficacement des index, donc toutes les jointures sont soient basées sur des hashs, soit sur des boucles imbriquées — le type de jointure le plus basique de PostgreSQL, qui est garanti de fonctionner à tous les coups mais qui peut être vite très lent. Quand une requête utilise des boucles imbriquées, soit le planificateur pense qu'il y a extrêmement peu de lignes et donc que les autres jointures ne valent pas la peine, soit il a juste abandonné devant la bizarrerie de la situation. Vu qu'il affiche un nombre de lignes très proche de ce qui est vraiment en base, et vu ce que je lui demande de toute façon, je penche pour la deuxième option.

Dans cette requête, on ne se contente pas d'utiliser all_map ; on le fait deux fois, et on fait une jointure de l'une sur l'autre, ainsi qu'avec la grille qui fait environ 100 000 pixels. Sachant que le résultat final, une agrégation, nécessitera que PostgreSQL récupère l'intégralité des lignes pour faire le compte, ça devient vite assez sanglant. On a cependant quelques actions simples qui peuvent améliorer drastiquement les performances et calmer les hurlements du serveur.

J'ai parlé de ST_EstimatedExtent tout à l'heure ; vu que la fonction accède aux statistiques de la table, ça veut dire qu'on ne ferait plus du tout appel aux lignes renvoyées par all_map mais seulement à l'estimation que PostgreSQL fait tout seul des données contenues dans la table. Si on peut trouver un moyen pour que notre vue soit compatible avec, on gagne déjà beaucoup. Vu qu'on a eu une erreur pour le type geography, on pourrait se dire que juste créer une sorte d'alias suffirait :

wifi=# CREATE VIEW lol AS (SELECT geo::geometry FROM all_map);
wifi=# SELECT ST_EstimatedExtent('lol', 'geo');
WARNING:  stats for "lol.geo" do not exist
st_estimatedextent
------------------

(1 rows)

La ligne vide ici correspond en fait à un NULL, et le message est assez explicite. La vue vient à peine d'apparaître, donc PostgreSQL n'a aucune statistique. On peut essayer de les générer de force avec VACUUM ANALYZE:

wifi=# vacuum analyze lol;
WARNING:  skipping "lol" --- cannot vacuum non-tables or special system tables
VACUUM

Dommage ! Il n'y a pas de statistiques disponibles pour les vues. Mais on peut faire plus intelligent que ça avec une vue matérialisée : c'est juste une table qui garde un lien vers la vue d'origine, ce qui permet plus tard si on le souhaite de faire un REFRESH MATERIALIZED VIEW pour la régénérer.

CREATE MATERIALIZED VIEW geometries AS SELECT geo::geometry FROM all_map;
VACUUM ANALYZE geometries;

Cette fois, pas de plaintes de PostgreSQL. La vue matérialisée s'est créée en 35 secondes et on peut continuer en réécrivant notre requête pour qu'elle utilise ce nouveau jouet :

SELECT grid.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 grid.geom;

Le retour de EXPLAIN est cette fois largement plus digeste. On garde encore une boucle imbriquée, mais il n'y en a plus qu'une seule ; on joint la vue à la grille et on fait une agrégation. J'ai aussi vérifié à la main et le retour de ST_EstimatedExtent est tout à fait cohérent ; vu que PostgreSQL vient tout juste de préparer la vue matérialisée, il a des statistiques toutes fraîches et très précises.

Le coût indiqué par le planificateur, qui est son unité de mesure pour choisir quoi faire pour exécuter une requête, est passé de 26 milliards à 630 millions. Ça reste assez conséquent et je n'ai toujours pas pu avoir la moindre idée du temps d'exécution de la requête tant c'était long.

On peut aller un peu plus loin avec la technique magique de l'index GIST qui marche presque à tous les coups lorsqu'on parle de PostGIS. Il n'y a pas de questions d'extension btree_gist comme dans mes précédentes requêtes; ici on ne parle que de géométries et GIST peut parfaitement travailler avec ça. Mon expérience personnelle sur ce projet et sur d'autres me fait dire que dès qu'on utilise la moindre fonction qui compare des géométries, ou si on groupe ou agrège par une géométrie, un index GIST qui ne coûte pas grand chose sera très utile.

CREATE INDEX geomgist ON geometries USING GIST (geo);

La requête n'est pas changée, et EXPLAIN me dit qu'un nœud Materialize, qui signifie que PostgreSQL doit concrétiser les lignes à cet endroit et les avoir quelque part en mémoire ou sur le disque pour pouvoir continuer, a disparu. On n'accède en fait même plus vraiment à la vue geometries et on utilise à la place directement geomgist dans la jointure. On atteint un coût proche de 82000, et la requête s'exécute enfin en un temps raisonnable : 9.975 secondes.

On peut encore appliquer une dernière modification : on groupe actuellement par chaque case de la grille, mais ST_SquareGrid ne nous renvoie pas que des polygones comme montré plus haut. Tout ce qu'il nous faut, c'est avoir quelque chose qui identifie de façon unique chaque case de la grille, et on peut grouper par ça, trier par ça pour récupérer nos pixels. On peut donc grouper par i et j au lieu de grouper par géométrie, ce qui évite des comparaisons coûteuses entre des polygones qui n'ont pas d'index GIST :

SELECT i, j, 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;

Le coût estimé pour cette requête n'a pas changé d'un pouce, mais le temps d'exécution est quand même passé à 3.902 secondes. Voilà qui est bien mieux et plus dans la norme des requêtes qu'on fait sur ce projet.

Conclusion

On a passé beaucoup trop de temps à construire juste une requête bien plus petite que celles auxquelles on a l'habitude dans ce projet, et ce n'est que le début. Dans le prochain épisode, on récupérera le fond de carte sur lequel on pourra afficher nos pixels, et pour ça, il va me falloir écrire des équations en LaTeX.

Commentaires

Il n'y a pour l'instant aucun commentaire. Soyez le premier !