Après la publication de la première partie du travail sur les Geopackages, le travail avec PostgreSQL/Postgis nous a permis une série d’avancées complémentaires. Vous trouverez dans cet article la mise à jour de la procédure d’import des fichiers S57 et dans un article postérieur les mises à jour des symbologies S57 sous QGis pour les Geopackages.
Mise à jour de l’import des fichiers S57 : ajout de l’attribut « purpose »
Dans l’article original on aborde le problème des doublons lors de la gestion de la base de données. Ces doublons peuvent être de deux types différents:
- des « vrais » doublons où tous les attributs sont les mêmes. Ils sont rares car ils proviennent de la superposition de deux emprises de cartes marines avec une échelle très proche. Par contre ils peuvent aussi résulter d’une erreur de rechargement en double d’un fichier S57. Ce type de doublon est traité avec le script fourni dans l’article d’origine qui n’est pas mis à jour.
- des « faux » doublons où l’information de la couche se retrouve en double sans que les identifiants des enregistrements soient forcément identiques.Ils apparaissent quand des zones cartographiées à différentes échelles sont chargées dans la même base de données. Ce type de doublon ne doit pas être supprimé, mais géré.
Dans le script fourni dans l’article initial nous ajoutons le nom du fichier et l’échelle aux tables geopackage. L’utilisation de l’échelle pour gérer les « faux » doublons s’avère assez compliquée:
Dans l’exemple qui suit, voici le résultat de l’affichage des données brutes pour le Goulet de la Rade de Brest qui se trouve sur deux cartes:
Pour la Rade de Brest, les deux cartes se chevauchent. L’échelle de l’une est de 10000 et celle de l’autre de 22500.
Outre l’échelle, il y aune information qui peut s’avérer très utile , la finalité (purpose) de la carte. C’est une valeur comprise entre 1 et 6 et qui correspond à l’objectif principal de la carte:
- 1 : Vue d’ensemble
- 2 : Généralités
- 3 : Côtière
- 4 : Approche
- 5 : Port
- 6 : Accostage
Les deux cartes qui se chevauchent dans notre exemple ont deux finalités différentes, l’une a une finalité=5 (Port), l’autre une finalité=6 (Accostage).
Si on filtre avec purpose=5, on obtient:
Et si on filtre avec purpose=6:
Dans la base de données du projet avec PostgreSQL/Postgis, nous avons importé 4500 fichiers S57.
Les valeurs min et max de l’échelle de toutes les cartes pour l’objectif 5 sont 3000 et 60000. Les valeurs min et max de l’échelle pour l’objectif 6 sont 2500 et 15000. On voit bien que les valeurs d’échelle des cartes les plus détaillées se retrouvent à l’intérieur des cartes de type 5.
Voici le résultat pour l’ensemble des finalités:
| Purpose | min_scale | max_scale |
|---|---|---|
| 1 | 325000 | 10000000 |
| 2 | 100000 | 1534076 |
| 3 | 50000 | 600000 |
| 4 | 12500 | 150000 |
| 5 | 3000 | 60000 |
| 6 | 2500 | 15000 |
Nous proposons donc d’ajouter, en plus du nom de fichier et de l’échelle, la finalité, comme attribut des tables geopackage. Cet attribut se retrouve dans la même table DSID utilisée pour récupérer l’échelle avec le nom de DSID_INTU.
Les fichiers .bat d’import restent les mêmes, sauf que le script Python appelé doit être le script 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_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_scale
def get_default_purpose(geopackage_path, default_enc):
# Ouvrir le GeoPackage en mode lecture seule
geopackage = ogr.Open(geopackage_path, 0)
default_purpose = '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_INTU correspondant à DSID_DSNM
dsid_layer.SetAttributeFilter(f"DSID_DSNM = '{default_enc}'")
feature = dsid_layer.GetNextFeature()
if feature:
default_purpose = feature.GetField('DSID_INTU')
return default_purpose
def 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)
default_purpose = get_default_purpose(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)
default_purpose_str = str(default_purpose)
# 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)
# Vérifier si le champ 'purpose' existe
purpose_index = table.FindFieldIndex('purpose', 1)
if purpose_index < 0:
# Ajouter le champ 'purpose' s'il n'existe pas
champ3 = ogr.FieldDefn('purpose', ogr.OFTReal)
champ3.SetWidth(10)
champ3.SetPrecision(2)
champ3.SetDefault(default_purpose_str) # Convertir en nombre flottant
table.CreateField(champ3)
# 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')
purpose_value = feature.GetField('purpose')
# 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 or purpose_value is None:
feature.SetField('enc_chart', default_enc)
feature.SetField('scale', default_scale)
feature.SetField('purpose', default_purpose)
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)Vous pouvez télécharger le script ici.
Copiez le script au même emplacement que les scripts de l’article original.
Vous pouvez soit remplacer les lignes des fichiers .bat
python c:/testgpkg/update_geopackage_dsid.py « %output_geopackage% » « !file! »
par
python c:/testgpkg/update_geopackage_dsid_prp.py « %output_geopackage% » « !file! »
Soit effacer le script update_geopackage_dsid.py et renommer update_geopackage_dsid_prp.py en update_geopackage_dsid.py ce qui vous évite de modifier les fichiers .bat.
Ajout d’un script Python pour filtrer les couches par « purpose »
Pour filtrer l’affichage de toutes les couches chargées dans le projet QGis, vous pouvez utiliser le script suivant:
# Importer le module QGIS
from qgis.core import QgsProject, QgsMapLayerType
# Définir le valeur de l'attribut "purpose" à filtrer
valeur_purpose = 3
# 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 de type vecteur
if couche.type() == QgsMapLayerType.VectorLayer:
# Vérifier si la couche a un champ nommé "purpose"
if couche.fields().indexFromName('purpose') != -1:
# Vérifier si le nom de la couche commence par l'un des préfixes spécifiés
if couche.name().startswith(('pt_', 'li_', 'pl_')):
# Définir le filtre sur l'attribut "purpose" égal à la valeur spécifiée
filtre = f'"purpose" = {valeur_purpose}'
couche.setSubsetString(filtre)
print(f"Filtre défini pour la couche {couche.name()}")
else:
print(f"La couche {couche.name()} ne commence pas par 'pt_', 'li_', ou 'pl_'")
else:
print(f"La couche {couche.name()} n'a pas d'attribut 'purpose'")
else:
print(f"La couche {couche.name()} n'est pas une couche vecteur")
Vous pouvez télécharger le script ici.
Pour éliminer les filtres sur Purpose vous pouvez utiliser le même script que pour l’échelle (unfilter.py)
Le projet de base de données S57 avec Geopackage est arrivé à son terme. Nous travaillons maintenant sur la mise en place d’une procédure équivalente en utilisant une base de données PostgreSQL/Postgis. Aidez-nous à mener à terme ce projet!
