lucidiot's cybrecluster

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

Lucidiot Informatique 2021-11-14
Réimplémentons un bouton de l'interface web d'OpenStreetMap !


Nous revoilà pour un nouvel opus de l'incroyable saga qui n'intéresse presque personne. Nous avons précédemment appris à compter dans des cases, et maintenant on va pouvoir commencer à essayer de récupérer un fond de carte par dessus lequel nous pourrons afficher nos cases.

Si vous utilisez l'interface Web d'OpenStreetMap, vous observerez peut-être le bouton de partage qui permet de partager de plusieurs façons l'endroit où vous êtes sur la carte et le point affiché : en copiant l'URL, en incluant un objet HTML au milieu d'une page Web, en partageant une URI geo: compatible avec beaucoup de logiciels de cartographie, ou en générant un rendu de la carte en tant qu'une seule image. Ce dernier est celui qui nous intéresse.

Déconstruction d'une URL d'export OpenStreetMap

OpenStreetMap fonctionne avec un système de tiles ou carreaux qui représentent des carrés de la carte, donc j'avais cru en faisant mes premières recherches que j'allais devoir bricoler quelque chose qui récupère des tiles et les assemble en faisant des calculs pour positionner correctement la grille sur l'image obtenue. Mais cet export change tout ; il a une URL qui ressemble à ça :

https://render.openstreetmap.org/cgi-bin/export
  ?bbox=
    -11.843261718750002,
    39.774769485295465,
    10.78857421875,
    52.38901106223458
  &scale=6145195
  &format=png

J'ai rajouté des sauts de ligne pour rendre l'URL plus lisible. Le lien n'a pas l'air très compliqué ; on demande un rectangle englobant la zone à exporter, un facteur d'échelle, et un format d'image. Mais comment obtient-on ce facteur d'échelle ?

Calcul de l'échelle

J'ai été fouiller dans le code source du site web d'OpenStreetMap pour trouver le code qui génère cette URL d'export et notamment le facteur d'échelle. Je suis tombé sur cette fonction, qui fait des calculs assez étranges avec des nombres magiques non définis.

function getScale() {
  var bounds = map.getBounds(),
    centerLat = bounds.getCenter().lat,
    halfWorldMeters = 6378137 * Math.PI * Math.cos(centerLat * Math.PI / 180),
    meters = halfWorldMeters * (bounds.getEast() - bounds.getWest()) / 180,
    pixelsPerMeter = map.getSize().x / meters,
    metersPerPixel = 1 / (92 * 39.3701);
  return Math.round(1 / (pixelsPerMeter * metersPerPixel));
}

C'est alors que j'ai sorti un éditeur LaTeX pour me représenter un peu mieux tous ces calculs en une seule équation qui détermine le facteur d'échelle :

Équation de calcul d'échelle originale
fig. 1: Équation de calcul d'échelle originale

Équation LaTeX

\frac{1}{
  \frac{pixel\ width}{
    6378137 \pi
    \times \cos(\frac{lat_{center} \times  \pi}{180})
    \times \frac{lat_{east} - lat_{west}}{180}
  }
  \times
  \frac{1}{92 \times 39.3701}
}

J'ai seulement reproduit les calculs faits dans toutes ces variables JavaScript et n'ai aucunement touché à la partie mathématique. Mais on peut maintenant commencer à la simplifier :

  1. Les deux grandes fractions en bas sont multipliées entre elles, donc on peut les fusionner facilement en une seule fraction dont le seul numérateur sera pixel width et le dénominateur est la multiplication des deux dénominateurs actuels.

  2. On se retrouve avec une fraction en 1 / fraction, ce qui veut dire l'inverse d'une fraction. On peut donc retirer le 1 / et inverser les numérateur et dénominateur de la « sous-fraction ».

  3. Le premier π sur la gauche peut aussi s'écrire π / 1, ce qui rend plus clairement visible qu'on peut l'inclure dans la fraction de droite vu qu'il y a une multiplication de fractions là aussi.

  4. On a donc deux fractions qui utilisent toutes les deux une multiplication par π et une division par 180. Cette construction est utilisée quand on peut convertir des degrés en radians. PostgreSQL nous fournit déjà une fonction RADIANS() pour faire ça, donc je vais l'utiliser dans le reste de la notation pour raccourcir. J'ai tout à fait le droit de faire ça puisqu'il n'existe aucun véritable standard dans la notation mathématique ; chacun fait ce qu'il veut, ce qui rend souvent beaucoup de travaux mathématiques encore moins compréhensibles.

Début de simplification
fig. 2: Début de simplification

Équation LaTeX

\frac{
  6378137
  \times \cos(radians(lat_{center}))
  \times radians(lat_{east} - lat_{west})
  \times 92 \times 39.3701
}{pixel\ width}

Ça commence à être un peu plus buvable. On peut cependant encore finioler un peu et réduire au minimum le nombre de ces constantes magiques : on commence par les multiplier toutes entre elles, puis on va utiliser le principe inverse de celui qu'on a utilisé pour fusionner des fractions ; on va sortir les constantes du numérateur pour les mettre à côté de la fraction.

Équation finale
fig. 3: Équation finale

Équation LaTeX

\frac{23101926018.3404}{pixel\ width}
\cos(radians(lat_{center}))
\times radians(lat_{east} - lat_{west})

On peut maintenant commencer à traduire tout ça en SQL. La latitude correspond à l'axe X dans PostGIS. On peut obtenir le barycentre d'une géométrie avec ST_Centroid, et la soustraction des deux latitudes des bords de la géométrie peuvent se faire en utilisant de nouveaux ST_XMax et ST_XMin. On a donc une fonction qui prend en paramètre un rectangle englobant et une taille d'image désirée en pixels et qui renvoie un facteur d'échelle acceptable pour OpenStreetMap :

CREATE OR REPLACE FUNCTION osm_scale (bbox BOX2D, width INTEGER DEFAULT 1000)
RETURNS DOUBLE PRECISION AS
'SELECT 23101926018.3404 * COS(RADIANS(ST_X(ST_Centroid(bbox)))) * RADIANS(ST_XMax(bbox) - ST_XMin(bbox)) / width'
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT
PARALLEL SAFE;

Sachez qu'il aurait bien sûr été tout à fait possible de représenter cette fonction avec l'équation originale :

CREATE OR REPLACE FUNCTION osm_scale (bbox BOX2D, width INTEGER DEFAULT 1000)
RETURNS DOUBLE PRECISION AS
'SELECT 1. / ((width / (6378137. * PI() * COS(RADIANS(ST_X(ST_Centroid(bbox)))) * (ST_XMax(bbox) - ST_XMin(bbox)) / 180.)) * 1. / (92. * 39.3701))'
LANGUAGE SQL
IMMUTABLE
RETURNS NULL ON NULL INPUT
PARALLEL SAFE;

Mais pourquoi faire ça quand je peux faire peur au lecteur avec un petit aparté sur des équations ?

Reconstruction d'une URL OpenStreetMap

Il ne reste plus qu'à créer une nouvelle fonction pour générer une URL OpenStreetMap. Pas besoin de LaTeX cette fois, c'est une question de formatage de chaîne de caractères :

CREATE TYPE osm_render_format AS ENUM ('png', 'jpg', 'pdf', 'svg');

CREATE OR REPLACE FUNCTION osm_render_url(
  bbox box2d,
  width integer DEFAULT 1000,
  format osm_render_format DEFAULT 'png'::osm_render_format
) RETURNS text AS
'SELECT format(
  ''https://render.openstreetmap.org/cgi-bin/export?bbox=%s,%s,%s,%s&scale=%s&format=%s'',
  ST_XMin(bbox),
  ST_YMin(bbox),
  ST_XMax(bbox),
  ST_YMax(bbox),
  osm_scale(bbox, width),
  format
)'
LANGUAGE SQL
IMMUTABLE
PARALLEL SAFE;

Le paramètre bbox est le rectangle englobant, qui demande les 4 coordonnées séparées par des virgules. On utilise notre fonction de facteur d'échelle pour... le facteur d'échelle, et on peut choisir parmi 4 formats. On va se contenter de PNG pour cette carte, c'est le format le plus simple.

Maintenant qu'il nous est possible de construire une URL pour télécharger l'image du fond de carte, il va falloir télécharger cette image, construire une image représentant notre nombre de réseaux, et les superposer. Il aurait peut-être été possible de faire quelque chose avec un raster GDAL, mais je n'avais pas envie de passer énormément de temps à apprendre l'utilisation des rasters dans PostGIS. À la place, j'ai fait quelque chose de bien étrange.

Télécharger l'image

Je peux maintenant commettre quelques crimes supplémentaires en essayant d'utiliser maintenant cette génération d'URL pour exporter une image :

$ psql -c "COPY (SELECT 'curl ''' || osm_render_url(ST_EstimatedExtent('geometries', 'geo')) || '''') TO STDOUT" | sh -x
+ curl 'https://render.openstreetmap.org/cgi-bin/export?bbox=-11.843261718750002,39.774769485295465,10.78857421875,52.38901106223458&scale=6145195&format=png'
<html>
<head>
<title>Error</title>
</head>
<body>
<h1>Error</h1>
<p>Missing or invalid token</p>
</body>
</html>

Cette construction étrange avec psql me permet d'exécuter une requête SQL qui renvoie une chaîne de caractères et d'en récupérer le résultat directement. Cette technique ne fonctionne pas toujours puisqu'il y a un format spécial utilisé par PostgreSQL, notamment pour gérer les valeurs nulles ou de multiples colonnes, mais dans ce cas précis ça fonctionne. Je génère une commande curl que j'envoie directement au shell, et j'utilise sh -x pour qu'on puisse voir quelle commande est exécutée, pour rendre les choses un peu plus explicites.

L'URL générée est valide mais pourtant OpenStreetMap me renvoie une erreur. Après avoir fait un petit tour dans la console de développeur en effectuant un export normal sur leur site, j'ai pu découvrir qu'un jeton TOTP (la même chose que ce qui fait fonctionner l'authentification à deux facteurs sur la plupart des sites) est fourni sous la forme d'un cookie. Et quelques essais avec curl plus tard, j'ai pu confirmer que c'est tout ce qu'il me manquait.

Générer une ligne de commande

Pour récupérer ce token dans un script shell, on peut faire une requête vers la page d'accueil d'OpenStreetMap, demander à curl de nous donner seulement les en-têtes HTTP, et récupérer juste la valeur de _osm_totp_token :

$ curl 'https://...' \
  -H "Cookie: _osm_totp_token=$(
    curl -s --fail -I https://www.openstreetmap.org/ |
    sed -En '/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p'
  )" \
  -o map.png

J'utilise une expression sed qui va rechercher _osm_totp_token et supprime tous les caractères de la ligne sauf 6 chiffres situés après un =. Avec l'argument -n, sed n'affiche rien par défaut, donc je peux avoir un comportement similaire à grep juste en ajoutant p après mon remplacement ; seule la ligne qui contiendra _osm_totp_token sera affichée après la substitution. Je repasse ensuite ces caractères dans un header de cookie pour la véritable requête. Enfin, j'ajoute au curl principal un -o map.png parce que le PNG, ça ne s'affiche pas très bien dans un terminal.

Convertissons à nouveau tout ça en une requête COPY imbriquée dans un script shell, accompagnée de sa quantité absurde d'échappements de caractères :

$ psql -c "COPY (SELECT
  'curl '''
  || osm_render_url(ST_EstimatedExtent('geometries', 'geo'))
  || ''' -H \"Cookie: _osm_totp_token=\$(
    curl -s --fail -I https://www.openstreetmap.org/ |
    sed -En ''/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p''
  )\" -o map.png'
) TO STDOUT WITH (FORMAT CSV, QUOTE E'\\t')" | sh -x
+ curl -s --fail -I https://www.openstreetmap.org/
+ sed -En '/_osm_totp_token/s/^.*=([0-9]+);.*$/\1/p'
+ curl 'https://render.openstreetmap.org/cgi-bin/export?bbox=-11.843261718750002,39.774769485295465,10.78857421875,52.38901106223458&scale=6145195&format=png' -H 'Cookie: _osm_totp_token=043230' -o map.png

On obtient cette fois-ci une carte parfaitement valide !

L'instruction COPY a des arguments un peu étranges maintenant : on ajoute WITH (FORMAT CSV, QUOTE E'\\t'). Les deux controbliques à la suite sont là parce que le tout est imbriqué dans une chaîne de caractères de shell POSIX avec des apostrophes doubles, mais et le reste ?

Par défaut avec COPY, le format texte spécial de PostgreSQL est utilisé, et il utilise la controblique comme caractère d'échappement. Pour avoir une chaîne de caractères que sh pourra comprendre, il faut donc s'en débarrasser. Je pourrais suivre le psql d'un sed juste pour faire le nettoyage, mais pourquoi faire ça quand je peux activer le mode FORMAT CSV ?

Je me retrouve cependant encore coincé : cette fois, les apostrophes doubles sont utilisées comme caractères d'échappement. CSV est très mal standardisé, et bien que son nom soit Comma-Separated Values, certains utilisent des points-virgule, des tabulations, des espaces, ou je ne sais quoi encore. Excel peut par exemple utiliser le point-virgule si votre système est en français et la virgule en anglais, sans que ce soit pour autant clairement documenté. Mais PostgreSQL suit les règles courantes du CSV sous Unix, où """a"", b", c signifie deux colonnes : une qui contient "a", b et une qui contient c.

Et comme CSV n'est pas du tout standardisé et tout le monde fait n'importe quoi, PostgreSQL vous autorise à tordre sa définition de CSV pour qu'elle corresponde à ce que vous préférez, donc on peut utiliser QUOTE pour changer le caractère qui délimite des chaînes de caractères, ou ESCAPE pour changer le caractère d'échappement.

J'ai défini QUOTE à la tabulation, ce qui est possible avec PostgreSQL en utilisant le préfixe E sur la chaîne de caractères SQL pour que PostgreSQL se mette à interpréter \t comme une tabulation. J'avais déjà utilisé cette technique pour importer des données TSV venues de nos scans Wi-Fi dans un article précédent.

Ainsi, on se retrouve avec une commande curl qui commence et se termine par des tabulations, et rien de plus, sans que ça ne perturbe le reste de notre ligne de commande. Les tabulations sont ignorées par sh, donc nous sommes tranquilles.

Conclusion

À cette étape de notre génération de heatmaps, il est vrai que toute cette partie sur la gestion des échappements de caractères est totalement inutile ; j'aurais juste pu faire un SELECT osm_render_url(...), copier le résultat moi-même dans un script ou même juste le mettre dans un navigateur qui s'occupera des cookies à ma place. Mais mon objectif finalement est non seulement de générer une heatmap mais aussi de pouvoir le faire aussi simplement que je fais les tâches courantes dans ce projet de warwalking : lancer un script shell avec quelques arguments, et rien de plus. Plus tard dans cette incroyable saga, le script ainsi généré se complexifiera largement plus.

Après avoir utilisé du SQL, du shell POSIX, du LaTeX et même un peu de Python pour faire des heatmaps, nous allons dans le prochain épisode découvrir un format d'image ésotérique...


Commentaires

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