lucidiot's cybrecluster

Importer du GPX dans SpatiaLite

Lucidiot — Informatique — 2021-06-17
Utiliser autant de SQL que possible, ce qui n'est pas toujours la meilleure idée…

Cet article a été bien plus long que je ne l'avais prévu, et assez technique. Et histoire de vous « rassurer », c'est le premier d'une petite série…

Dans le cadre d'un projet de cartographie de réseaux Wi-Fi dont je parlerai plus en détail prochainement, j'ai eu besoin d'importer des fichiers GPX dans une base de données SpatiaLite. Les fichiers GPX peuvent décrire des waypoints (des points uniques nommés, avec des métadonnées), des routes (des chemins allant de waypoint à waypoint; des itinéraires) et des tracks (des enregistrements de la position GPS réelle à intervalles réguliers). Dans mon cas, je n'ai que des traces, et je ne m'intéresse même pas à isoler chaque trace vu qu'il peut y en avoir plusieurs dans le même fichier ; je veux juste récupérer les coordonnées et la précision de tous les points, avec leur horodatage.

Solutions existantes

Il y a moult façon de procéder à cela. La première option de DuckDuckGo est d'utiliser gpx2spatialite, un package Python qui importe tout dans une structure de tables assez compliquée mais d'où je pourrais extraire les informations dont j'ai besoin. Ce package utilise cependant Python 2, et je n'ai pas envie d'avoir à réinstaller Python 2 partout sur mes machines pour pouvoir m'en servir, surtout qu'idéalement j'aimerais pouvoir exécuter cet import sur Termux, où les paquets disponibles sont assez restreints et installer des dépendances peut vite être long et complexe.

L'option suivante serait d'utiliser du GeoJSON comme format intermédiaire : passer de GPX à GeoJSON à l'aide d'un script, puis écrire un one-liner avec jq pour extraire les bonnes valeurs et générer un CSV ou des commandes INSERT qu'on peut envoyer à SQLite. Les seuls convertisseurs que j'ai pu trouver sont soit des services en ligne ne m'inspirant pas trop confiance, soit des librairies NodeJS. JavaScript n'est pas tellement le problème dans mon projet vu que je n'essaie pas de faire du code très solide, mais c'est l'écosystème NPM qui me pose problème, donc je préfère éviter. Je ne veux pas non plus avoir à utiliser un navigateur pour faire cette étape intermédiaire, ce qui demanderait beaucoup de travail manuel.

Une autre option encore serait d'écrire ma propre version de gpx2spatialite, et je pourrais très bien faire une XSLT pour transformer le GPX en instructions SQL, ou utiliser mon script Python existant de conversion XML vers JSON puis utiliser un script jq. Mais tout ça implique encore et toujours des dépendances. N'y a-t-il pas un moyen de tout faire avec SpatiaLite ?

Le package Debian spatialite-bin fournit moult outils, y compris spatialite_xml_load. Je ne sais pas trop pourquoi ils ont fait ça, mais cet outil peut créer autant de tables que nécessaire dans une base de données pour représenter l'intégralité d'un fichier XML. En important mon GPX, je me retrouve avec des dizaines de tables telles que gpx_trk_trkseg_trkpt_time et les jointures deviennent plutôt complexes. Ne parlons même pas du nettoyage ; je ne peux pas utiliser une base en mémoire comme espace temporaire pour préparer un CSV que je transmettrais ensuite à ma vraie base, donc je me retrouve à soit devoir gérer des fichiers temporaires ou supprimer moult tables de la base. Il doit forcément y avoir quelque chose de plus simple.

VirtualXPath

Vu que c'est ma première fois avec cette extension, j'essaie juste de jouer un peu avec SpatiaLite pour comprendre un peu plus comment ça marche, tandis que je continue à réfléchir à ce problème. Je me demande pourquoi l'interpréteur de SpatiaLite m'affiche à chaque fois une vingtaine de lignes au lieu de juste trois lignes comme celui standard de SQLite :

SpatiaLite version ..: 4.3.0a   Supported Extensions:
    - 'VirtualShape'    [direct Shapefile access]
    - 'VirtualDbf'      [direct DBF access]
    - 'VirtualXL'       [direct XLS access]
    - 'VirtualText'     [direct CSV/TXT access]
    - 'VirtualNetwork'  [Dijkstra shortest path]
    - 'RTree'       [Spatial Index - R*Tree]
    - 'MbrCache'        [Spatial Index - MBR cache]
    - 'VirtualSpatialIndex' [R*Tree metahandler]
    - 'VirtualElementary'   [ElemGeoms metahandler]
    - 'VirtualXPath'    [XML Path Language - XPath]
    - 'VirtualFDO'      [FDO-OGR interoperability]
    - 'VirtualGPKG' [OGC GeoPackage interoperability]
    - 'VirtualBBox'     [BoundingBox tables]
    - 'SpatiaLite'      [Spatial SQL - OGC]
PROJ.4 version ......: Rel. 6.3.1, February 10th, 2020
GEOS version ........: 3.8.0-CAPI-1.13.1 
TARGET CPU ..........: x86_64-linux-gnu
the SPATIAL_REF_SYS table already contains some row(s)
SQLite version ......: 3.31.1
Enter ".help" for instructions
SQLite version 3.31.1 2020-01-27 19:55:54
Enter ".help" for instructions
Enter SQL statements terminated with a ";"
spatialite>

Je lis un peu plus attentivement ces lignes et je m'arrête sur VirtualXPath. Je sais qu'il existe des extensions de SQLite pour ouvrir des fichiers, d'autres bases de données, ou même des sites web comme des tables ; peut-être que cette extension qui parle de XML peut m'aider ?

SpatiaLite rajoute le XmlBlob, qui est un type spécial permettant de représenter des documents XML dans un format binaire. Ce n'est pas sans me rappeler le format CXML utilisé par Sony pour certaines métadonnées dans les systèmes PSP et PS3, et qu'il me faudra finir de documenter un jour.

XmlBlob permet de gérer des métadonnées standardisées des systèmes d'information géographique, ou des informations de style pour afficher des données géométriques avec des couleurs différentes. L'extension VirtualXPath rajoute en plus de ça la possibilité d'utiliser une table un peu particulière pour filtrer des documents XML en utilisant la notation XPath. Ça tombe bien, je maîtrise XPath vu que j'ai dû m'en servir pour ITSB.

Pour importer du GPX, on pourrait donc juste envoyer le GPX directement à la base de données et le traiter en SQL, comme je pourrais le faire dans PostgreSQL ? Ça semble totalement maudit, donc ça me convient parfaitement !

Réinventons la roue

J'ai mis quelques heures à développer un script shell compatible POSIX qui génère des instructions pour l'interpréteur SpatiaLite pour cet import. J'ai notamment perdu pas mal de temps parce qu'une colonne affichée comme vide peut être un NULL, une chaîne de caractères vide, ou toute valeur qui ne peut pas être représentée sous forme de texte, telle qu'un XmlBlob. On va lire pas à pas toutes ces requêtes.

Charger le XML

D'abord, il nous faut charger le XML. Pour avoir un XmlBlob, la seule façon est d'utiliser la fonction XB_Create, qui demande un BLOB. L'encodage d'un fichier XML est parfois décrit dans son header <?xml encoding="us-ascii"?>, donc il est important d'utiliser un format binaire pour laisser à l'interpréteur XML la capacité d'utiliser n'importe quel encodage. Cela dit, je m'en fiche un peu puisque tout ce que j'utilise sera en UTF-8, donc je vais utiliser du texte, et le transformer en BLOB pour faire plaisir à SpatiaLite :

CREATE TEMPORARY TABLE gpx_import_data AS SELECT XB_Create(CAST('<gpx>…</gpx>' AS BLOB)) AS document;

Vu que SQLite n'a pas de notion de variables et plus généralement pas de notion de procédures stockées, on utilise une table temporaire d'une ligne et d'une colonne pour tout stocker.

Activer XPath

On doit ensuite créer une table virtuelle pour utiliser XPath. Une table virtuelle pourrait être appelée une vue dans tous les autres moteurs de base de données, mais SQLite y ajoute de la syntaxe particulière qui permet en fait de choisir un "moteur" qui fera fonctionner cette table. Ici, on va utiliser VirtualXPath et lui donner le nom de notre table et la colonne où se trouve le fichier GPX :

CREATE VIRTUAL TABLE gpx_import_xpath USING VirtualXPath('gpx_import_data', 'document');

Utiliser XPath

Nous voulons récupérer deux choses différentes dans chaque point des traces GPX : le contenu de deux balises, et la valeur de deux attributs. Voici à quoi ressemble un unique point dans mon fichier GPX :

<trkpt lat="32.8656424" lon="8.754631">
  <ele>122</ele>
  <time>2021-05-31T17:21:28Z</time>
  <hdop>10</hdop>
  <extensions>
    <osmand:speed>0.7</osmand:speed>
  </extensions>
</trkpt>

On a les deux coordonnées en tant qu'attributs, et l'élévation au dessus du sol en mètres, l'horodatage, la précision horizontale et la vitesse comme balises.

Étant donné que je me fiche complètement d'à quelle trace appartient chaque point, je peux ignorer l'intégralité de la hiérarchie et aller chercher tout de suite tous les points, puis récupérer ces quatre attributs. Cela se traduit en une expression XPath assez concise : //trkpt/(@lat | @lon | hdop | time). Le // indique qu'on fait une recherche récursive. On s'arrête à la balise trkpt, et on y cherche deux attributs et deux balises à l'aide d'unions (|).

Il est possible de faire des requêtes en lecture sur gpx_import_xpath, à condition de toujours indiquer une expression XPath à utiliser ; ne pas mettre d'expression du tout ne renvoit rien du tout. Je ne m'attendais pas vraiment à ce comportement qui n'est pas documenté, donc ça m'a un peu confus au début. On obtient donc :

SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//trkpt/(@lat | @lon | hdop | time)';

Gérer les espaces de noms

Si on essaie d'exécuter la requête ci-dessus, on se rend compte que rien du tout n'est renvoyé. En fait, même juste //trkpt ou /gpx ne renvoient rien. Pourtant, faire la même expression dans une XSLT ou dans l'interpréteur XPath intégré aux navigateurs fonctionne parfaitement. Mais pourquoi ?

Mon fichier GPX utilise moult espaces de noms. En XML, les espaces de noms permettent de mélanger plusieurs schémas ; disons par exemple que vous vouliez inclure dans votre fichier GPX des informations sur des personnes à l'aide de FOAF : c'est possible en ajoutant un attribut xmlns:foaf pointant vers le schéma XML de FOAF, puis en préfixant toutes les balises avec foaf:, par exemple <foaf:Person>. Si vous voulez obtenir un élément appartenant à un espace de noms dans une expression XPath, il faudra rajouter son préfixe aussi : //foaf:Person. Mais si l'élément n'a pas d'espace de noms, ou est dans l'espace de noms par défaut (celui défini par xmlns="…", alors vous ne devriez pas avoir à vous préoccuper de ça.

Dans mon GPX, tous les éléments standards des GPX utilisent l'espace de noms par défaut, et seule l'extension ajoutée OsmAnd pour indiquer la vitesse utilise un espace de noms. Mais l'extension VirtualXPath a décidé que s'il y a des espaces de noms, alors tout le monde devrait en avoir un ! Toutes les autres balises se retrouvent alors préfixées par dflt, abréviation de default. Je suis donc obligé de modifier l'expression XPath :

SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/(@lat | @lon | dflt:hdop | dflt:time)';

Nous n'avons pas les mêmes valeurs

Toujours aucune valeur n'est renvoyée pour les balises hdop et time, donc je n'arrive pas à récupérer l'horodatage et la précision ! En essayant un peu au hasard et en lisant attentivement les quelques exemples dans la documentation, je me rends compte que demander explicitement la valeur texte des éléments renvoie bien la valeur : dflt:time/text() au lieu de seulement dflt:time. La raison pour cela est assez logique, mais pas forcément visible dans tous les contextes où on utilise XPath.

Dans XPath, il y a quatre types : les nœuds ou ensembles de nœuds, les chaînes de caractères, les nombres et les booléens. Mais les nœuds ont eux-mêmes différents types : ce sont les balises, les attributs, les commentaires, etc. d'un fichier XML. Alors que la notion de valeur pour les chaînes, nombres et booléens est assez explicite, celle pour les nœuds dépend du type de nœud. Un attribut a une valeur, c'est même quasiment la seule chose qu'il a. Sachant qu'un commentaire n'est qu'une chaîne de caractères délimitée, sa valeur est assez logique à voir. Mais pour une balise XML, ça n'est pas si facile.

La valeur de <a>12</a> pourrait être interprétée comme étant "12", mais qu'en est-il d'une structure plus complexe comme <a><b>12</b>34</a> ? La valeur serait "1234" ? "34" ? "<b>12</b>34" ? Ou peut-être même l'élément tout entier ? Et quid de <a ></a> ?

Dans une XSLT, où on transforme du XML en XML, on ne se pose pas trop la question : si on utilise <xsl:value-of="." ></xsl:value-of="."> par exemple, on mettra <b>12</b>34 directement. La valeur est un ensemble de nœuds, ce que le générateur peut tout à fait accepter et produire en sortie.

Je m'attendais donc à la même chose dans mes requêtes SQL, mais je me suis heurté à un problème mentionné en passant tout à l'heure : ce qui n'a pas de représentation directe en chaîne de caractères pour SQLite est affiché comme vide. Je ne pouvais juste pas exploiter la valeur du tout.

On allonge donc encore un peu la requête :

SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/(@lat | @lon | dflt:hdop/text() | dflt:time/text())';

La sortie de l'Union

La requête ne renvoie toujours rien ! Mais qu'est-ce que j'ai encore raté ?

Je n'ai pas creusé beaucoup, mais il semblerait que VirtualXPath ne prenne pas du tout en charge les unions. Même avec une forme plus explicite telle que //dflt:trkpt/@lat | //dflt:trkpt/@lon ne fonctionne pas. Tant pis, on va faire plusieurs requêtes séparées pour récupérer chaque type de valeur :

SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:time/text()';
SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lat';
SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lon';
SELECT * FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:hdop/text()';

Mais on se retrouve alors avec 4 requêtes complètement indépendantes, et on ne sait plus à quel point appartient chaque donnée ! Il nous faut un moyen de les relier à nouveau.

Quand on veut ajouter des lignes de deux requêtes les unes après les autres, c'est une UNION. Quand on veut récupérer les lignes communes à deux requêtes, c'est un INTERSECT. Quand on veut retirer les lignes d'une requête d'une autre requête, c'est un EXCEPT. Mais lorsqu'il s'agit de manipuler des colonnes, peu importe ce qu'on veut faire, c'est un JOIN.

Jusqu'ici, j'ai utilisé SELECT *, mais regardons un peu plus en détail ce que ces requêtes XPath nous renvoient :

pkid
La clé primaire de la ligne qui contient le XmlDocument dans la table originale gpx_import_data.
sub
Un index supplémentaire, parce qu'on peut avoir plusieurs correspondances d'une expression XPath dans le même document.
parent
Nom du nœud parent, s'il existe.
node
Nom du nœud sur lequel on a eu une correspondance.
attribute
Nom de l'attribut sur lequel on a eu une correspondance, si l'expression XPath devait trouver des attributs.
value
Valeur du nœud ou de l'attribut trouvé.
xpath_expr
L'expression qu'on a exécuté.

Dans mon contexte, je n'ai qu'une seule ligne et donc un seul document XML dans ma table temporaire, donc je peux ignorer pkid et seul sub compte un identifiant unique. Je n'ai pas besoin de xpath_expr vu que c'est moi qui la tape, ni de node ou attribute parce que je sais ce que je cherche. Et je n'ai pas besoin de parent non plus car la recherche des résultats avec XPath traversera le document toujours dans le même ordre, donc si je fais deux requêtes séparées, le sub correspondra à la même position dans les résultats. Je n'ai donc besoin que de sub et value ; value aura le résultat de mes requêtes, et sub me permettra de faire toutes les jointures entre les requêtes.

Pour une raison que je n'ai pas vraiment pu déterminer, faire une jointure classique (INNER JOIN) n'a pas fonctionné, donc j'ai dû utiliser une LEFT JOIN. Vu que j'ai toujours le même nombre de résultats dans toutes les requêtes et que le sub correspond toujours, je n'ai pourtant pas besoin de gérer la moindre valeur nulle. Cette affirmation ne serait pas vraie si je n'avais pas un fichier GPX bien structuré, où des attributs seraient manquants…

On se retrouve donc avec cette nouvelle requête :

SELECT * FROM
    (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:time/text()') time
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lat') lat USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lon') lon USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:hdop/text()') hdop USING (sub);

Insérer dans la table

Dernière ligne droite ! Il faut maintenant récupérer les valeurs que je veux et en faire des lignes, une pour chaque point, dans les bons types. Pour l'instant, je n'ai que des chaînes de caractères.

SQLite a une gestion assez particulière des dates, en utilisant la notion de jour julien : le nombre de jours à partir du 1er janvier -4712 à 12h00 UTC. Il est cependant très flexible dans les données qu'il peut recevoir et transformer en ce format, donc je peux lui donner directement ma date au format ISO 8601 : On utilisera donc JULIANDAY(time.value) pour obtenir un REAL qui plaira à SQLite.

En ce qui concerne la position, je veux garder un point tel qu'utilisé par SpatiaLite, et pas seulement deux coordonnées séparées. Il me faudra donc d'abord les transformer en nombre à virgule double précision avec CAST(lat.value AS DOUBLE), et je peux utiliser la fonction MakePoint(<latitude>, <longitude>) pour construire un point.

Cela dit, si je crée des points ainsi, j'ai utilisé des coordonnées cartésiennes : la Terre est alors considérée comme complètement plate ! Quoi que vous puissiez croire, ça ne plaît pas à ce système et les cartes générées ainsi pourraient avoir des points placés incorrectement avec une déviation allant jusqu'à des milliers de kilomètres. Un cercle tracé sur la France par exemple aurait la forme d'un œuf.

La solution est d'utiliser un référentiel. Les systèmes d'information géographiques acceptent souvent en entrée un très grand nombre de coordonnées, qui bien qu'elles soient toutes des chiffres peuvent correspondre à des positions très différentes. Certains systèmes de positionnement par satellite ont leur propre référentiel, parfois conçu pour volontairement rendre inutilisable le positionnement (celui de la Chine par exemple). Il existe plusieurs bases de données qui listent ces référentiels, et celle qui nous intéresse le plus est l'EPSG Geodetic Parameter Dataset, du nom du groupe européen d'exploration pétrolière. L'EPSG n'existe plus et la base est aujourd'hui maintenue par le comité géomatique de l'association internationale des producteurs de pétrole et de gaz (IOGP). L'industrie pétrolière simplifie le boulot des géographes !

Le référentiel que tout le monde utilise quand on parle de « coordonnées GPS », c'est WGS 84, qui a aussi le doux nom de EPSG 4326. J'ai fini par apprendre ce dernier numéro par cœur parce qu'on finit par l'utiliser partout quand on joue dans des bases de données géographiques. Je vais donc indiquer aussi à SpatiaLite que mes coordonnées utilisent ce référentiel avec MakePoint(<latitude>, <longitude>, 4326), et mes points seront ainsi bien placés.

Enfin, pour la précision, c'est un nombre entier positif, donc je peux juste le convertir en entier : CAST(hdop.value AS INTEGER). Tout ça n'est pas vraiment obligatoire vu que SQLite est très loin d'être strict dans sa gestion des types dans les tables, mais je préfère toujours être bien explicite pour éviter d'avoir des données peu consistantes et des surprises plus tard.

Je peux donc sélectionner ces trois colonnes et tout insérer en une seule fois dans ma table finale avec une construction INSERT SELECT :

INSERT INTO trip (time, position, precision)
SELECT JULIANDAY(time.value), MakePoint(CAST(lat.value AS DOUBLE), CAST(lon.value AS DOUBLE), 4326), CAST(hdop.value AS INTEGER)
FROM
    (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:time/text()') time
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lat') lat USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lon') lon USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:hdop/text()') hdop USING (sub);

Script final

On peut maintenant prendre tout ce bazar et le mettre dans un petit script shell que je pourrai utiliser pour insérer dans ma table :

#!/bin/sh
# Generates SQL statements to import a GPX file from stdin to SpatiaLite
printf "BEGIN;
CREATE TEMPORARY TABLE gpx_import_data AS SELECT XB_Create(CAST('"
sed "s/'/''/g"
echo "' AS BLOB)) AS document;
DROP TABLE IF EXISTS gpx_import_xpath;
CREATE VIRTUAL TABLE gpx_import_xpath USING VirtualXPath('gpx_import_data', 'document');
INSERT INTO trip (time, position, precision)
SELECT JULIANDAY(time.value), MakePoint(CAST(lat.value AS DOUBLE), CAST(lon.value AS DOUBLE), 4326), CAST(hdop.value AS INTEGER)
FROM
    (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:time/text()') time
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lat') lat USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/@lon') lon USING (sub)
    LEFT JOIN (SELECT sub, value FROM gpx_import_xpath WHERE xpath_expr = '//dflt:trkpt/dflt:hdop/text()') hdop USING (sub);
DROP TABLE gpx_import_xpath;
COMMIT;"

Avec un fichier GPX nommé data.gpx et une base de données déjà existante nommée wifi.sqlite, on peut l'utiliser ainsi :

$ ./gpx_import.sh < data.gpx | spatialite wifi.sqlite

Et c'est ainsi que je me retrouve avec des tables remplies de points en n'utilisant que du SQL, avec un peu de shell pour économiser quelques appuis de touches.

Et spoiler : J'ai fait tout ça pour rien ! On verra ça la prochaine fois…

Commentaires

fluffy, 2021-06-18

La dernière ligne de cet article est la meilleur.