La première partie du projet Financement Collaboratif pour l’Intégration de Données Marines dans QGIS est arrivée à son terme grâce aux contributions d’ORANGE Marine, Geoconceptos Uruguay, Janez Avzec et Gaetan Mas. Nous les remercions vivement. Nous publions donc le contenu final de cette partie du travail.
Etant donné le volume du résultat, nous le publions en deux parties: la première concerne la création et la gestion de la base de données, la deuxième la symbologie S57 sous QGis.
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 et ne permet pas la création de procédures ou fonctions.
Si vous avez un seul 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 et ajout des attributs carte_enc et echelle..
- 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_.
Si vous avez un lot de fichiers S57 à charger simultanément, le processus est le suivant:
- Chargement des classes d’objets Point dans le geopackage pointsENC avec ogr2ogr avec suppression des tables vides et ajout du nom du fichier enc et de l’échelle dans les attributs de toutes les tables.
- Chargement des classes d’objets Lignes dans le geopackage LinesENC avec ogr2ogr avec suppression des tables vides et ajout du nom du fichier enc et de l’échelle dans les attributs de toutes les tables.
- Chargement des classes d’objets Polygone dans le geopackage PolysENC avec ogr2ogr avec suppression des tables vides et ajout du nom du fichier enc et de l’échelle dans les attributs de toutes les tables.
- 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 (optionnel)
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 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:
-oo SPLIT_MULTIPOINT=ON -oo ADD_SOUNDG_DEPTH=ON
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:
-oo RETURN_LINKAGES=O
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 :
-mapFieldType StringList=String,IntegerList=String
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:
- Exécutez la commande ogr2ogr pour créer le geopackage d’import
- Repérez la valeur maximale indiquée pour chaque type de champ
- Pour la table ou les tables en sortie qui auront été créées et dont vous souhaitez récupérer toute l’information, ouvrez le gestionnaire de base de données de QGis
- Effacez tous les enregistrements de la table (DELETE FROM nom_de la table)
- Allez dans le menu Table->Modifiez une table…
Sélectionnez dans la liste des champs le champ qui vous concerne, puis cliquez sur Editer une colonne - Rentrez la valeur souhaitée dans le champ Longueur, puis sur OK
- Faites la même opération pour tous les autres champs qui vous intéressent
- Effacez du geopackage d’import toutes les autres tables qui ne posent pas de problème
- Ré-exécutez la commande ogr2ogr d’import.
Workflow pour le chargement d’un seul fichier S57
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…
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:
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='POINT' or OGR_GEOMETRY='MULTIPOINT'" -oo SPLIT_MULTIPOINT=ON -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -oo ADD_SOUNDG_DEPTH=ON -mapFieldType StringList=String,IntegerList=String C:/testgpkg/pointsENC.gpkg C:\testgpkg\US1GC09M\US1GC09M.000
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='LINESTRING' or OGR_GEOMETRY='MULTILINESTRING'" -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -mapFieldType StringList=String,IntegerList=String C:/testgpkg/LinesENC.gpkg C:\testgpkg\US1GC09M\US1GC09M.000
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='POLYGON' or OGR_GEOMETRY='MULTIPOLYGON'" -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -mapFieldType StringList=String,IntegerList=String C:/testgpkg/PolysENC.gpkg C:\testgpkg\US1GC09M\US1GC09M.000
Vous devrez modifier le texte
C:/testgpkg/pointsENC.gpkg C:\testgpkg\US1GC09M\US1GC09M.000
par l’emplacement des geopackages d’import que vous souhaitez et par le chemin et nom de fichier du fichier S547 à charger.
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.
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 haut, 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:
from osgeo import ogrdef delete_empty_tables_in_geopackages(geopackage_paths, valeur_par_defaut_enc, valeur_par_defaut_scale):
for geopackage_path in geopackage_paths:
delete_empty_tables(geopackage_path, valeur_par_defaut_enc, valeur_par_defaut_scale)
print(f"Fin du traitement")
def delete_empty_tables(geopackage_path, valeur_par_defaut_enc, valeur_par_defaut_scale):
# Ouvrir le GeoPackage en mode édition
geopackage = ogr.Open(geopackage_path, 1)if geopackage is not None:
# Récupérer le nombre de tables dans le GeoPackage
num_tables = geopackage.GetLayerCount()
# Liste des tables à supprimer
tables_to_delete = []
tables_to_update = []
# Identifier les tables à supprimer
for i in range(num_tables):
table = geopackage.GetLayerByIndex(i)
# Vérifier si la table est vide (aucun enregistrement)
if table.GetFeatureCount() == 0:
tables_to_delete.append(table.GetName())
else:
tables_to_update.append(table.GetName())
# Supprimer les tables
for table_name in tables_to_delete:
geopackage.DeleteLayer(table_name)
print(f"Table supprimée dans {geopackage_path}: {table_name}")
# Mettre à jour les tables restantes
for table_name in tables_to_update:
table = geopackage.GetLayerByName(table_name)
# Ajouter le premier champ texte
champ1 = ogr.FieldDefn('enc_chart', ogr.OFTString)
champ1.SetWidth(50)
champ1.SetNullable(True)
champ1.SetDefault(valeur_par_defaut_enc)
table.CreateField(champ1)
# Ajouter le deuxième champ texte
champ2 = ogr.FieldDefn('scale', ogr.OFTReal)
champ2.SetWidth(10)
champ2.SetPrecision(2)
champ2.SetDefault(valeur_par_defaut_scale)
table.CreateField(champ2)
# Mettre à jour les valeurs dans les champs après leur création
for feature in table:
feature.SetField('enc_chart', valeur_par_defaut_enc)
feature.SetField('scale', valeur_par_defaut_scale)
table.SetFeature(feature)
print(f"ENC & scale ajoutés à {geopackage_path}: {table_name}")
# Fermer le GeoPackage
geopackage = None
else:
print(f"Impossible d'ouvrir le GeoPackage {geopackage_path} en mode édition.")#Lignes à éditer avant l'exécution du scriptgeopackage_paths = ["c:/testgpkg/pointsENC.gpkg", "c:/testgpkg/LinesENC.gpkg", "c:/testgpkg/PolysENC.gpkg"]
valeur_par_defaut_enc = 'US1GC09M'
valeur_par_defaut_scale = '2160000'
delete_empty_tables_in_geopackages(geopackage_paths, valeur_par_defaut_enc, valeur_par_defaut_scale)
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:
Workflow pour le chargement de plusieurs fichiers S57 en même temps
Les lignes de commandes des paragraphes précédents 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.
Contrairement au chargement d’un seul fichier, lors du chargement d’un lot de fichiers il faut intégrer l’ajout du nom du fichier et de l’échelle à la boucle de transcodage S57->geopackage. Pour cela, nous utiliserons un script python qui ajoute les attributs enc_chart et scale aux tables créées par la commande ogr2ogr et qui les renseigne avec le nom du fichier correspondant et une échelle définie par le producteur. Le script efface aussi les tables vides présentes dans le geopackage d’import.
Merci à Janez AVSEC pour l’idée d’utiliser la table DSID pour récupérer l’échelle de compilation des données.
Les fichiers .bat doivent être exécutés à partir d’une fenêtre OSGeo4W pour pouvoir bénéficier de l’environnement Python de QGis.
L’appel de chaque fichier .bat est de la forme:
.\fichier.bat « répertoire contenant les fichiers S57 à charger » « nom du geopackage d’import »
Par exemple:
.\mass_load_s57_polys.bat « c:/enc_db/enc000 » « c:/enc_db/polys.gpkg »
Les fichiers .bat pour réaliser cette tâche sont les suivants:
Pour importer les tables de type point:
@echo off
setlocal enabledelayedexpansionREM Assurez-vous qu'au moins 2 arguments ont été fournis
if "%~2"=="" (
echo Usage: %0 directory000 output_geopackage scale
exit /b 1
)REM Récupérer les arguments de la ligne de commande
set "directory=%~1"
set "output_geopackage=%~2"REM Itérez sur tous les fichiers .000 dans le répertoire
for /r "%directory%" %%i in (*.000) do (
echo Traitement du fichier: %%~ni
set "file=%%~ni"
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='POINT' or OGR_GEOMETRY='MULTIPOINT'" -oo SPLIT_MULTIPOINT=ON -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -oo ADD_SOUNDG_DEPTH=ON -mapFieldType StringList=String,IntegerList=String "%output_geopackage%" "%%i"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "DSID"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "C_AGGR"
python c:/testgpkg/update_geopackage_dsid.py "%output_geopackage%" "!file!"
)echo Traitement terminé.
pause
Vous remarquerez que, en plus de la commande ogr2ogr pour importer les couches de points, on importe deux tables sans géométrie: DSID et C_AGGR.
Pour importer les tables de type ligne:
@echo off
setlocal enabledelayedexpansionREM Assurez-vous qu'au moins 2 arguments ont été fournis
if "%~2"=="" (
echo Usage: %0 directory000 output_geopackage scale
exit /b 1
)REM Récupérer les arguments de la ligne de commande
set "directory=%~1"
set "output_geopackage=%~2"
set "scale=%~3"REM Itérez sur tous les fichiers .000 dans le répertoire
for /r "%directory%" %%i in (*.000) do (
echo Traitement du fichier: %%~ni
set "file=%%~ni"
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='LINESTRING' or OGR_GEOMETRY='MULTILINESTRING'" -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -mapFieldType StringList=String,IntegerList=String "%output_geopackage%" "%%i"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "DSID"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "C_AGGR"
python c:/testgpkg/update_geopackage_dsid.py "%output_geopackage%" "!file!"
)echo Traitement terminé.
pausePour importer les tables de type polygone:
@echo off
setlocal enabledelayedexpansionREM Assurez-vous qu'au moins 2 arguments ont été fournis
if "%~2"=="" (
echo Usage: %0 directory000 output_geopackage scale
exit /b 1
)REM Récupérer les arguments de la ligne de commande
set "directory=%~1"
set "output_geopackage=%~2"
set "scale=%~3"REM Itérez sur tous les fichiers .000 dans le répertoire
for /r "%directory%" %%i in (*.000) do (
echo Traitement du fichier: %%~ni
set "file=%%~ni"
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='POLYGON' or OGR_GEOMETRY='MULTIPOLYGON'" -oo RETURN_LINKAGES=ON -oo LNAM_REFS=ON -mapFieldType StringList=String,IntegerList=String "%output_geopackage%" "%%i"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "DSID"
ogr2ogr -f GPKG -skipfailures -append -update "%output_geopackage%" "%%i" "C_AGGR"
python c:/testgpkg/update_geopackage_dsid.py "%output_geopackage%" "!file!"
)echo Traitement terminé.
pause
Les fichiers .bat contiennent la commande ogr2ogr ainsi qu’une ligne pour exécuter un script Python, update_geopackage_dsid.py. Faites attention à modifier le chemin de ce script dans le fichier .bat pour qu’il corresponde à celui de votre choix. Ce script ajoute à toutes les tables du geopackage les deux attributs enc_chart et scale.
Le code python utilisé est le suivant:
import sys
from osgeo import ogr
def get_default_scale(geopackage_path, default_enc):
# Ouvrir le GeoPackage en mode lecture seule
geopackage = ogr.Open(geopackage_path, 0)
default_scale = ‘0’
default_enc = default_enc +’.000′
if geopackage is not None:
# Récupérer la couche DSID
dsid_layer = geopackage.GetLayerByName('DSID')
# Vérifier si la couche DSID existe
if dsid_layer:
# Rechercher la valeur de DSID_CSCL correspondant à DSID_DSNM
dsid_layer.SetAttributeFilter(f"DSID_DSNM = '{default_enc}'")
feature = dsid_layer.GetNextFeature()
if feature:
default_scale = feature.GetField('DSPM_CSCL')
return default_scaledef delete_empty_tables_in_geopackages_dsid(geopackage_path, default_enc):
# Récupérer la valeur de scale à partir de la table DSID
default_scale = get_default_scale(geopackage_path, default_enc)
# Continuer avec le reste du code
if default_scale is None:
print("La valeur par défaut de l'échelle n'a pas été trouvée dans la table DSID.")
return
default_scale_str = str(default_scale)
# Ouvrir le GeoPackage en mode édition
geopackage = ogr.Open(geopackage_path, 1)
if geopackage is not None:
# Récupérer le nombre de tables dans le GeoPackage
num_tables = geopackage.GetLayerCount()
# Liste des tables à supprimer
tables_to_delete = []
tables_to_update = []
# Identifier les tables à supprimer
for i in range(num_tables):
table = geopackage.GetLayerByIndex(i)
# Vérifier si la table est vide (aucun enregistrement)
if table.GetFeatureCount() == 0:
tables_to_delete.append(table.GetName())
else:
tables_to_update.append(table.GetName())
# Supprimer les tables
for table_name in tables_to_delete:
geopackage.DeleteLayer(table_name)
print(f"Table supprimée dans {geopackage_path}: {table_name}")
# Mettre à jour les tables restantes
for table_name in tables_to_update:
table = geopackage.GetLayerByName(table_name)
# Vérifier si le champ 'enc_chart' existe
enc_chart_index = table.FindFieldIndex('enc_chart', 1)
if enc_chart_index < 0:
# Ajouter le champ 'enc_chart' s'il n'existe pas
champ1 = ogr.FieldDefn('enc_chart', ogr.OFTString)
champ1.SetWidth(50)
champ1.SetNullable(True)
champ1.SetDefault(default_enc)
table.CreateField(champ1)
# Vérifier si le champ 'scale' existe
scale_index = table.FindFieldIndex('scale', 1)
if scale_index < 0:
# Ajouter le champ 'scale' s'il n'existe pas
champ2 = ogr.FieldDefn('scale', ogr.OFTReal)
champ2.SetWidth(10)
champ2.SetPrecision(2)
champ2.SetDefault(default_scale_str) # Convertir en nombre flottant
table.CreateField(champ2)
# Mettre à jour les valeurs dans les champs seulement si les champs n'ont pas déjà de valeurs
for feature in table:
enc_chart_value = feature.GetField('enc_chart')
scale_value = feature.GetField('scale')
# Vérifier si les champs n'ont pas déjà de valeurs définies
if enc_chart_value is None or scale_value is None:
feature.SetField('enc_chart', default_enc)
feature.SetField('scale', default_scale)
table.SetFeature(feature)
# Fermer le GeoPackage
geopackage = None
else:
print(f"Impossible d'ouvrir le GeoPackage {geopackage_path} en mode édition.")if name == « main« :
if len(sys.argv) != 3:
print(« Usage: python script.py geopackage_path default_enc »)
sys.exit(1)
geopackage_path = sys.argv[1]
default_enc = sys.argv[2]
delete_empty_tables_in_geopackages_dsid(geopackage_path, default_enc)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.
Vous pouvez télécharger les fichiers .bat et le code Python en cliquant ici.
Suite quel qu’il soit le nombre de fichiers chargés
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:
from osgeo import ogrdef clone_or_append_tables_with_prefix():
input_geopackages = ["c:/testgpkg/pointsENC.gpkg", "c:/testgpkg/LinesENC.gpkg", "c:/testgpkg/PolysENC.gpkg"]
prefixes = ["pt_", "li_", "pl_"]for i in range(len(input_geopackages)):
# Ouvrir le GeoPackage source en mode lecture
input_gpkg = ogr.Open(input_geopackages[i], 0)
if input_gpkg is not None:
# Ouvrir le GeoPackage de destination en mode écriture
output_gpkg = ogr.Open("c:/testgpkg/ENC.gpkg", 1)
prefix = prefixes[i]
if output_gpkg is not None:
# Récupérer le nombre de tables dans le GeoPackage source
num_tables = input_gpkg.GetLayerCount()
# Ajouter le contenu des tables avec le préfixe approprié
for j in range(num_tables):
input_table = input_gpkg.GetLayerByIndex(j)
output_table_name = f"{prefix}{input_table.GetName()}"
input_table_name = input_table.GetName()
# Vérifier si la table de destination existe déjà
output_table = output_gpkg.GetLayerByName(output_table_name)
if output_table is None:
# Créer une nouvelle couche de destination si elle n'existe pas
output_table = output_gpkg.CreateLayer(output_table_name, geom_type=input_table.GetGeomType(), options=["OVERWRITE=YES"])
# Copiez la structure de la table d'origine vers la table temporaire
input_layer = input_gpkg.GetLayerByName(input_table_name)
input_layer_defn = input_layer.GetLayerDefn()
for k in range(input_layer_defn.GetFieldCount()):
field_defn = input_layer_defn.GetFieldDefn(k)
output_table.CreateField(field_defn)
# Copier les entités de la table source vers la table de destination
for feature in input_table:
output_feature = ogr.Feature(output_table.GetLayerDefn())
output_feature.SetGeometry(feature.GetGeometryRef())
# Copier les valeurs des champs
for field_index in range(feature.GetFieldCount()):
output_feature.SetField(field_index, feature.GetField(field_index))
output_table.CreateFeature(output_feature)
output_feature = None # Libérer la mémoire
print(f"Contenu ajouté de {input_geopackages[i]}.{input_table.GetName()} vers ENC.gpkg.{output_table_name}")
# Fermer les GeoPackages
input_gpkg = None
output_gpkg = None
else:
print(f"Impossible d'ouvrir le GeoPackage {output_geopackage} en mode écriture.")
else:
print(f"Impossible d'ouvrir le GeoPackage {input_geopackages[i]} en mode lecture.")#Appel de la fonction append_tables_with_prefixclone_or_append_tables_with_prefix()
Veillez à modifier la ligne
input_geopackages = [« c:/testgpkg/pointsENC.gpkg », « c:/testgpkg/LinesENC.gpkg », « c:/testgpkg/PolysENC.gpkg »]
avec les geopackages d’import que vous avez choisi, ainsi que la ligne
output_gpkg = ogr.Open(« c:/testgpkg/ENC.gpkg », 1)
avec le geopackage de base de données de votre choix.
Vous pouvez télécharger le fichier .py avec ce lien: clone_or_append_tables_with_prefix.py.
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.
Cette opération correspond au filtrage d’une couche, mais… que faire quand vous en avez 130 couches affichées? Pour vous éviter de faire 130 fois cette opération, voici un code Python qui vous permet de filtrer toutes les couches du projet:
#Importer le module QGISfrom qgis.core import QgsProject#Définir la valeur de l'attribut "scale" à filtrervaleur_echelle = 50000#Obtenir le projet actifprojet = QgsProject.instance()#Obtenir la liste des couches chargées dans le projetcouches = projet.mapLayers().values()#Parcourir toutes les couchesfor couche in couches:
# Vérifier si la couche a un champ nommé "scale"
if couche.fields().indexFromName('scale') != -1:
# Définir le filtre sur l'attribut "scale" égal à 50000
filtre = f'"scale" = {valeur_echelle}'
couche.setSubsetString(filtre)
print(f"Filre défini pour la couche {couche.name()}")
else:
print(f"La couche {couche.name()} n'a pas d'attribut 'scale'")
Et un autre pour enlever le filtre à toutes les couches affichées:
# Importer le module QGIS
from qgis.core import QgsProject
# Obtenir le projet actif
projet = QgsProject.instance()
# Obtenir la liste des couches chargées dans le projet
couches = projet.mapLayers().values()
# Parcourir toutes les couches
for couche in couches:
# Vérifier si la couche est une couche vectorielle
if couche.type() == QgsMapLayer.VectorLayer:
# Supprimer le filtre de la couche
couche.setSubsetString('')
print(f"Filtre supprimé pour la couche {couche.name()}")
Télécharger les codes python ici.
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 LNAM. 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:
- aux tables des feux (LIGHTS) et des sondes bathymétriques (SOUNDG) où l’attribut LNAM n’est pas utile. En effet, chaque secteur de feu a un même LNAM et
name_rcidmais correspond à une couleur est un secteur géographique différents. De même, le LNAM etname_rciddes 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êmename_rcid, et pour les sondes nous devons trouver les sondes ayant les mêmes XY dans leur géométrie. - aux tables SEAARE, LNDRGN et ADMARE qui couvrent des vastes zones et qui sont découpées selon l’emprise de la carte, tout en gardant le même LNAM. Dans ce cas même si le LNAM est le même, les polygones n’ont pas la même forme et surface.
- les tables DSID et C_AGGR qui n’ont pas de LNAM et qui, par prinipe, ne peuvent pas avoir des doublons.
from osgeo import ogr
def delete_doublons():
# Ouvre le GeoPackage en mode mise à jour
geopackage_path = « c:/enc_db/enc2.gpkg »
geopackage = ogr.Open(geopackage_path, update=True)
if geopackage is not None:
# Récupère la liste de noms de toutes les tables du GeoPackage
table_names = [geopackage.GetLayerByIndex(i).GetName() for i in range(geopackage.GetLayerCount())]
# Boucle sur les tables et les traite individuellement
for table_name in table_names:
if table_name not in ["pl_SEAARE","pl_LNDRGN","pl_ADMARE","layer_styles","natsurf","DSID","C_AGGR"]:
# Crée une table temporaire en ajoutant un suffixe "_temp"
temp_table_name = table_name + "_temp"
temp_table = geopackage.CreateLayer(temp_table_name)
# Copie la structure de la table d'origine vers la table temporaire
input_layer = geopackage.GetLayerByName(table_name)
input_layer_defn = input_layer.GetLayerDefn()
for i in range(input_layer_defn.GetFieldCount()):
field_defn = input_layer_defn.GetFieldDefn(i)
temp_table.CreateField(field_defn)
if table_name not in ["pt_LIGHTS", "pt_SOUNDG"]:
attribut = "LNAM"
geopackage.ExecuteSQL(f"INSERT INTO {temp_table_name} SELECT * FROM {table_name} GROUP BY {attribut};")
elif table_name == "pt_SOUNDG":
attribut = "name_rcid,depth, ST_X(geom),ST_Y(geom)"
geopackage.ExecuteSQL(f"INSERT INTO {temp_table_name} SELECT * FROM {table_name} GROUP BY {attribut};")
elif table_name == "pt_LIGHTS":
attribut = "name_rcid, colour, sectr1, sectr2"
geopackage.ExecuteSQL(f"INSERT INTO {temp_table_name} SELECT * FROM {table_name} GROUP BY {attribut};")
geopackage.ExecuteSQL(f"DROP TABLE {table_name}")
# Renomme la table temporaire pour qu'elle ait le même nom que la table d'origine
geopackage.ExecuteSQL(f"ALTER TABLE {temp_table_name} RENAME TO {table_name}")
print(f"alter table {table_name}")
# Ferme le GeoPackage
geopackage = None
else:
print(f"Impossible d'ouvrir le GeoPackage {geopackage} en mode écriture.")#Appel de la fonction delete_doublonsdelete_doublons()
Veillez à modifier la ligne
geopackage_path = « c:/testgpkg/ENC.gpkg »
avec votre geopackage de base de données.
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 [« pl_SEAARE », »pl_LNDRGN », »pl_ADMARE », »layer_styles », »natsurf », »DSID », »C_AGGR »]:
En principe, vous aurez que les tables préfixées et les tables layer_styles et natsurf. Cette ligne permet de ne pas traiter pl_SEAARE, pl_LNDRGN et pl_ADMARE, mais surtout, elle empêche de vider les tables layer_styles et natsurf que nous utilisons pour les symbologies. 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 à cette ligne d’exceptions.
Le projet de base de données S57 avec Geopackage est arrivé à son terme. Nous entamons maintenant la mise en place d’une procédure équivalente en utilisant une base de données PostgreSQL/Postgis. Aidez-nous à mener à terme ce projet!
