Cet article sous sa forme complète (codes et liens) sera un des résultats du travail résultant du financement collaboratif que nous avons mis en place: https://www.sigterritoires.fr/index.php/financement-collaboratif-pour-lintegration-de-donnees-marines-dans-qgis/. Dès que le montant visé sera atteint, l’article sera mis à jours dans sa version complète.
Dans le cas où nous n’atteindrions pas le montant visé, nous enverrons la version complète à tous les contributeurs.
Le contenu de l’article indique, au jour d’aujourd’hui, le contenu du projet pour les fichiers Geopackage, à l’issue des multiples échanges avec nos contributeurs.
…………………………………………………………………………………………………..
Le format S57 est complexe et ne respecte pas totalement les standards des formats SIG. Sa conversion avec GDAL est donc soit complexe, soit incomplète. Dans cet article nous allons rentrer dans le détail de la transformation complète du S57 en tables de fichier Geopackage.
Les types de géométrie
Les « couches » S57 sont des classes d’objets. Par exemple, les différentes zones terrestres sont codées dans la classe d’objets LNDARE.
La définition de cette classe d’objets est:
Geo object: Land area (LNDARE) (P,L,A)
Attributes: CONDTN OBJNAM NOBJNM STATUS INFORM NINFOM
Si tous les objets LNDARE possèdent les mêmes attributs, il n’en est pas de même du type de géométrie. L’information (P,L,A) indique que, dans cette classe d’objets on peut trouver des points, des lignes et des polygones. Contrairement aux standards SIG, les trois types de géométrie coexistent au sein de la même classe d’objets.
Quel qu’il soit le format choisi pour l’intégrer dans QGis, nous aurons à créer une couche par type de géométrie. Si on ne le fait pas, GDAL va créer le type de couche en fonction de la première entité trouvée lors de la conversion. Si elle est de type point, la couche créée sera de type point et les entités lignes et polygones seront ignorées. De même, si la première est de type ligne, ce seront les points et les polygones qui seront ignorés.
Il faut aussi noter que le format S57 n’a aucune contrainte sur les classes d’objets: si dans la zone couverte par le fichier S57 il n’y a pas une classe d’objet, il n’y aura rien la concernant dans le fichier (pas de couche vide). De même, s’il n’y a qu’un type de géométrie présente, alors que les trois types sont possibles, il n’y aura pas de trace des autres types de géométrie.
Il faudra alors décomposer le traitement d’un fichier S57 avec ogr2ogr en trois étapes, une pour chaque type de géométrie. Les options suivantes permettent de traiter chaque classe d’objet S57 en ne sélectionnant qu’un type de géométrie :
-where « OGR_GEOMETRY=’POINT’ or OGR_GEOMETRY=’MULTIPOINT’ »
-where « OGR_GEOMETRY=’LINESTRING’ or OGR_GEOMETRY=’MULTILINESTRING’ »
-where « OGR_GEOMETRY=’POLYGON’ or OGR_GEOMETRY=’MULTIPOLYGON’ »
GDAL permet pour certains types de drivers de créer des préfixes dans les tables en sortie. Dans ce cas on pourrait créer toutes les tables (points, lignes, polygones) dans un seul geopackage, en les préfixant par exemple avec pt_, li_ et pl_. Le problème est que le driver S57 de GDAL ne permet pas cette option. On doit alors créer trois geopackages distincts, dans lesquels les tables seront créées selon le type de géométrie. Dans chaque geopackage il y aura alors une table de même nom mais de géométrie différente. Si cette solution est la plus simple, elle n’est pas la plus élégante et pratique. Nous allons utiliser ici trois geopackages d’import pour les commandes ogr2ogr dans lesquels nous allons importer les tables d’un fichier S57. Nous les appellerons pointsENC, LinesENC et PolysENC. D’autre part nous allons créer un geopackage unique appelé ENC pour la base de données complète. Pour chaque carte ENC nous allons importer son contenu dans les trois geopackages d’import avec ogr2ogr puis nous allons exécuter des scripts Python pour mettre à jour la base de données du geopackage ENC avec les nouvelles tables importées. Nous sommes obligés d’utiliser les scripts Python au lieu de requêtes et fonctions SQL car la version SQL de SQLite, utilisée pour les geopackages, est très limitée.
Pour le premier fichier S57 à importer le processus est le suivant:
- Création du geopackage ENC, vide.
- Chargement des classes d’objets Point dans le geopackage pointsENC avec ogr2ogr
- Chargement des classes d’objets Lignes dans le geopackage LinesENC avec ogr2ogr
- Chargement des classes d’objets Polygone dans le geopackage PolysENC avec ogr2ogr
- Suppression des tables vides des trois geopackages.
- Clonage des tables des trois geopackages d’import dans le geopackage ENC, en préfixant les tables de pointsENC par pt_, de LinesENC par li_ et celles de PolysENC par pl_.
Pour les fichiers S57 suivants, le processus est le suivant:
- Chargement des classes d’objets Point dans le geopackage pointsENC avec ogr2ogr
- Chargement des classes d’objets Lignes dans le geopackage LinesENC avec ogr2ogr
- Chargement des classes d’objets Polygone dans le geopackage PolysENC avec ogr2ogr
- Suppression des tables vides des trois geopackages.
- Mise à jour des tables existantes dans le geopackage ENC à partir des trois geopackages d’import ET
- Clonage des tables des trois geopackages d’import qui n’étaient pas présentes dans le geopackage ENC,
- Détection et effacement des doublons éventuels résultant de la mise à jour des tables du geopackage ENC
Le traitement des sondes bathymétriques
Les sondes bathymétriques posent plusieurs problèmes lors de la conversion de format. Le premier est que la valeur de profondeur n’est pas une contenue dans un attribut, mais en tant que Z dans la géométrie (XYZ). Le deuxième c’est que les sondes ne sont pas du type Point mais Multipoint. Pour avoir directement la valeur des sondes dans un champ attributaire, il faut donc ajouter deux paramètres à la ligne de commande ogr2ogr:
Récupération des identifiants « parents »
Certains objets peuvent être regroupés par un objet non géographique parent. Par exemple, les secteurs de feu sont codés avec un secteur par enregistrement de la table « LIGHTS ». Les différents secteurs d’un même feu n’ont aucun attribut particulier qui puisse indiquer qu’ils correspondant à la même entité. Cette information se trouve dans un enregistrement « parent ». Pour récupérer cet identifiant en tant qu’attribut de la table LIGTHS, il faut ajouter une option à la ligne de commande:
Cette option crée un attribut name_rcid qui est commun à tous les secteurs d’un même feu, mais créé aussi une série de champs (name_rcnm,ORNT,USAG,MASK).
Le traitement des type « Liste »
Le format S57 utilise, en plus des types de champ classiques (entier, réel, texte) un type particulier: les listes, de chaînes de caractères ou d’entiers. Si on retrouve ces types de champ dans PostgreSQL, il n’en va pas de même dans les shapefiles et le geopackage.
En absence des options spécifiques, il y aura des messages d’avertissement indiquant que le format n’est pas supporté et le contenu de ces champs sera ignoré. Pour ne pas perdre le contenu, il faut ajouter dans la ligne de commande ogr2ogr l’option :
Cette option prend une liste {..,…} et génère une chaîne de texte avec un premier chiffre indiquant le nombre d’éléments de la liste, deux points, puis les éléments de la liste séparés par des virgules:
{3} -> ‘(1:3)’
{3,1} -> ‘(2:3,1)’
Forcément, ceci complique un peu le traitement de ces champs mais l’information n’est pas perdue. Par contre, le driver GDAL présente un problème pour la gestion du format String résultant.
Tout d’abord, pour ce qui est des attributs des objets S57 il n’y a aucun problème. Le champ de la table geopackage est de type ‘Text’, tout court, ce qui permet la création d’une chaîne de caractères de n’importe quelle longueur.
Par contre, pour les « identifiants parents », la définition utilisée par le driver est ‘Text()’, c’est à dire une chaîne de caractères avec une longueur maximum. Cette longueur est, a priori, déterminée par le contenu du premier enregistrement traité. C’est dire que ce n’est que du hasard si vous avez un champ avec la longueur nécessaire pour intégrer toutes les valeurs du fichier.
Vous verrez alors des messages du type:
WARNING : Value of field ‘NAME_RCID’ has 371 characters, whereas maximum allowed is 30
La solution n’est pas compliquée mais implique pas mal de temps de travail. Partez de la base que pour l’affichage classique des cartes ENC, ceci ne représente pas un problème et vous pouvez tout simplement ignorer ces messages d’avertissement.
Par contre si vous avez besoin de traitements plus poussés incluant ces informations, voici le seul moyen que j’ai trouvé, pour le moment:
Commandes ogr2ogr pour créer les geopackages d’import
Fort de tous ces éléments, voici donc les trois lignes de commandes ogr2ogr pour créer les tables dans les geopackages d’import:
Commandes ogr2ogr
Pour exécuter ces lignes de commande il suffit d’ouvrir la fenêtre de shell d’OSGeo4W:
La fenêtre de Shell s’ouvre
Rentrez alors les lignes de commande, en faisant attention à bien rentrer le chemin et nom de fichier avec l’extension .000 .
On aura alors comme résultat une série de tables créées dans chaque geopackage d’import. Par contre il y aura des tables complètement vides. En effet, si le type de géométrie est possible pour une classe d’objet, la table sera créée mais s’il n’y a aucune occurrence dans le fichier S57 traité, la table n’aura aucun enregistrement.
Charger plusieurs fichiers S57 en même temps
Les lignes de commandes du paragraphe précédent assument que vous avez un seul fichier S57 à charger. Mais vous pouvez charger plusieurs fichiers simultanément. L’utilisation des options -append -update font que vous pouvez traiter en boucle tous les fichiers .000 d’un répertoire. Vous aurez ainsi tous les fichiers chargés dans les geopackages d’import.
Les fichiers .bat pour réaliser cette tâche sont les suivants:
Fichier .bat pour charger les géométries point
Fichier .bat pour charger les géométries ligne
Fichier .bat pour charger les géométries polygone
Faites attention à modifier les chemins pour qu’ils correspondent à vos geopackages d’import et au répertoire à charger. On pourrait faire un seul fichier.bat pour les trois traitements, mais il convient de bien vérifier le résultat de chacun des types de géométrie.
Création du geopackage ENC
Nous allons créer un geopackage vide qui recevra les tables importées, puis les mises à jour successives lors des chargements d’autres fichiers S57.
Dans le panneau Explorateur de QGis, faites un clic droit sur Geopackages et sélectionnez Créer une base de données
Dans la fenêtre qui s’ouvre, renseignez le nom du geopackage, laissez l’option par défaut pour le nom de la table et sélectionnez Pas de géométrie dans Type de géométrie.
Cliquez sur OK pour avoir le nouveau geopackage vide ajouté aux connexions geopackage de QGis. Il est conseillé d’effacer la table vide ENC qui est créée lors de la création du geopackage. Pour cela ouvrez le gestionnaire de base de données de QGis à partir du menu principal, positionnez-vous sur la table ENC du geopackage ENC.gpkg et avec un clic droit sélectionnez Effacer…
Utilisation des scripts Python pour la gestion de la base de données
Pour toutes les opérations de gestion de la base de données nous allons utiliser la console Python de QGis.
Pour l’ouvrir, allez dans le menu Extensions -> Console Python. Les panneaux de la console Python s’ouvrent:
Pour exécuter un script Python, chargez le code dans le panneau Éditeur. Vous pouvez faire un copier-coller ou bien utiliser l’icône Répertoire pour pointer sur le fichier .py. L’indentation étant cruciale pour les scripts Python, privilégiez cette deuxième option.
Suppression des tables vides et ajout du nom de la carte et de l’échelle
Comme nous l’avons vu plus haute, lors de l’exécution des commandes ogr2ogr, des tables vides sont créées pour les couches où il n’y a pas la géométrie demandée dans la ligne de commande. Nous allons donc supprimer ces tables vides et profiter du script Python pour ajouter deux attributs à toutes les tables: le nom de la carte ENC et l’échelle de celle-ci.
Pour ce qui est du nom, il peut être utile dans la maintenance ultérieure de la base de données de pouvoir sélectionner les entités correspondantes à un fichier source S57. Cette information n’est contenue dans aucun attribut du fichier S57. Nous devons donc la renseigner manuellement.
Il en va de même pour l’échelle de la carte. Si un attribut est bien prévu dans le S57 pour rentrer l’échelle de la carte « papier », l’attribut CSCALE dans la table (M_CSCL, il est bien rare que cette table et son attribut soient présents dans les fichiers S57.
Si le nom de la carte est d’une utilité relative, il n’en va pas de même de l’échelle, car c’est un critère indispensable dans la recherche de doublons. En effet, lors du mélange des entités en provenance de plusieurs fichiers S57 il y aura une coexistence plus ou moins heureuse de données issues de différentes échelles.
Le code Python suivant permet de supprimer les tables vides d’un schéma d’import et d’ajouter deux attributs (enc_chart et scale) à toutes les tables:
Attention à respecter l’indentation du code Python. Pour simplifier la tâche, téléchargez directement le fichier .py avec ce lien: delete_empty_tables_in_geopackages.py.
Pour l’exécution, veillez à modifier les paramètres en entrée du script:
Les geopackage_paths vous n’aurez qu’à les rentrer que la première fois que vous exécuterez le script. Par contre le nom de la carte ENC et l’échelle seront à modifier lors de l’import de chaque nouveau fichier S57.
A l’issue de cette exécution vous aurez dans chaque table les nouveaux attributs:
Clonage et/ou ajout des tables importées dans le geopackage ENC
Dans le cas du premier chargement S57, le script Python suivant créé toutes les tables présentes dans les trois geopackages d’import. Pour les chargements ultérieurs, si la table existe déjà dans le geopackage ENC, les enregistrements de la table du geopackage d’import sont ajoutés aux enregistrements déjà présents dans ENC. Si, par contre, la table n’existe pas encore, elle est créée et alimentée par les enregistrements du geopackage d’import. Dans le fichier geopackage ENC les tables sont préfixées par pt_, li_ et pl_ pour indiquer leur type de géométrie.
Le clonage est réalisé avec le code suivant:
Vous pouvez télécharger le fichier .py avec ce lien: clone_or_append_tables_with_prefix.py.
Ajout de la symbologie au Geopackage
La symbologie par défaut des tables d’un geopackage est contenue dans une table nommée layer_styles qui n’est créée que si vous utilisez l’option enregistrer par défaut du bouton style.
Après le clonage des tables dans le fichier ENC, la table n’existe pas. Nous allons donc importer la table contenant toute la symbologie des cartes ENC. Pour cela; téléchargez le geopackage layer_styles en cliquant ici.
Ce geopackage ne contient qu’une table, nommée layer_styles, avec les symbologies de toutes les tables S57.
Cette symbologie utilise des symboles svg que vous devez télécharger sur votre machine. En cliquant ici vous téléchargerez un répertoire ‘nautical’ contenant tous les symboles svg nécessaires.
Par défaut, les références aux symboles svg sont faites vers un répertoire ‘C:/nautical’. Vous avez plusieurs possibilités:
1- Décompresser le fichier nautical dans un répertoire c:/nautical de votre machine
Dans ce cas vous n’avez plus rien à faire. Les layer_styles trouveront les symboles svg sans problème.
2- Vous souhaitez enregistrer le répertoire nautical dans un répertoire de votre choix
Dans ce cas, vous devez modifier la table layer_styles pour qu’elle trouve les symboles.
Ouvrez le gestionnaire de base de données de QGis.
- Sélectionnez la table layer_styles de votre geopackage
- Ouvrez une fenêtre SQL
- Entrez la requête:
UPDATE layer_styles
SET styleQML = REPLACE(styleQML, ‘C:/nautical’, ‘Votre chemin’); - Cliquez sur le bouton Exécuter
3- Vous souhaitez enregistrer le répertoire nautical dans un le répertoire de votre geopackage ENC
Dans ce cas, les références aux symboles svg seront relatives au Geopackage et vous pourrez l’utiliser à partir de plusieurs machines. On fera alors la même démarche que précédemment mais la requête SQL sera:
UPDATE layer_styles
SET styleQML = REPLACE(styleQML, ‘C:/nautical’, ‘nautical‘);
Bien sûr, il faut que le projet soit configuré pour stocker les chemins relatifs (Propriétés du projet -> Général->Enregistrer les chemins->Relatif)
Gestion des doublons
Quand deux cartes ENC se superposent, les entités présentes dans la zone de superposition vont être insérées autant de fois que de superpositions il y en a. Il faut donc, pour éviter des doublons et plus, nettoyer la base de données après chaque import de fichier S57.
Mais il y a deux types de doublons différents. D’une part, les zones de recouvrement de deux cartes d’échelle similaire vont produire de « vrais » doublons qu’il faudra nettoyer. Mais la superposition de cartes d’échelle très différente est àç traiter différemment. La carte moins détaillé (petite échelle) aura une sélection des données de la carte plus détaillé (grande échelle). Dans ce deuxième cas il n’est pas conseillé de supprimer cette information.
Le code Python ci-après va supprimer les « vrais » doublons. Pour ce qui est des autres doublons, il est conseillé de travailler avec des filtres QGis. Selon vos besoins, vous pourrez définir, grâce à l’attribut « scale » que nous avons ajouté lors de la suppression des tables vides, quelles sont les données affichées ou pas sur votre carte.
Dans cet exemple, seules les entités correspondantes à une échelle comprise entre 12500 et 45000 seront affichées sur votre carte.
Le code suivant crée, pour chaque table du geopackage ENC, une table temporaire (temp_table). Puis il remplit cette table avec les enregistrements distincts basés sur name_rcid. Ensuite, il supprime la table d’origine et renomme la table temporaire avec le nom de la table d’origine.
Ce traitement s’applique à toutes les tables S57, sauf les tables des feux (LIGHTS) et des sondes bathymétriques (SOUNDG). En effet, chaque secteur de feu a le même name_rcid mais correspond à une couleur est un secteur géographique différents. De même, le name_rcid des sondes bathymétriques est le même pour tous les points de sondes de l’entité multipoint d’origine. Dans le premier cas nous devons chercher les couleurs et secteurs en double pour un même name_rcid, et pour les sondes nous devons trouver les sondes ayant les mêmes XY dans leur géométrie.
Vous pouvez télécharger le fichier .py à partir de ce lien : supprime_doublons.py. N’oubliez pas de modifier les chemins et noms de fichiers pour qu’ils correspondent aux vôtres. Remarquez aussi la ligne :
if table_name not in [« pt_LIGHTS », « pt_SOUNDG », »layer_styles »]:
En principe, vous aurez que les tables préfixées et la table layer_styles. Cette ligne permet de traiter différemment les tables LIGTHS et SOUNDG, mais surtout, elle empêche de vider la table layer_styles. Si vous êtes amené à créer d’autres tables que celles issues du fichier s57 dans le geopackage, il faut prendre en compte les noms de ces tables dans ce script, en les ajoutant aux exceptions.
