Following the publication of the first part of our work on Geopackages, our work with PostgreSQL/Postgis has enabled us to make a series of further advances. In this article, you’ll find an updated procedure for importing S57 files, and in a later article, updates to S57 symbologies in QGis for Geopackages.
S57 file import update: “purpose” attribute added
In the original article, we discussed the problem of duplicates in database management. These duplicates can be of two different types:
- true” duplicates, where all the attributes are the same. These are rare, as they result from the superimposition of two chart overlays of very similar scale. On the other hand, they can also be the result of an error in reloading a duplicate S57 file. This type of duplicate is handled by the script provided in the original article, which is not updated.
- false” duplicates, where the layer information is duplicated without the record identifiers necessarily being identical; these occur when areas mapped at different scales are loaded into the same database. This type of duplication should not be deleted, but managed.
In the script provided in the initial article, we add the file name and scale to the geopackage tables. Using the scale to manage “false” duplicates turns out to be quite complicated:
In the following example, here’s the result of displaying the raw data for the Goulet de la Rade de Brest on two maps:
For the Rade de Brest, the two maps overlap. The scale of one is 10000 and that of the other 22500.
In addition to the scale, there’s one piece of information that can be very useful: the purpose of the map. This is a value between 1 and 6, corresponding to the main objective of the map:
- 1 : Overview
- 2 : General
- 3 : Coastal
- 4 : Approach
- 5 : Port
- 6: Docking
The two overlapping maps in our example have two different purposes, one with purpose=5 (Port), the other with purpose=6 (Berthing).
If we filter with purpose=5, we get:
And if we filter with purpose=6:
In the PostgreSQL/Postgis project database, we imported 670 S57 files.
The min and max values for the scale of all maps for goal 5 are 3000 and 30000. The min and max scale values for objective 6 are 6000 and 15000. We can see that the scale values of the most detailed maps are found within the type 5 maps.
We therefore propose to add, in addition to the file name and scale, the purpose, as an attribute of the geopackage tables. This attribute is found in the same DSID table used to retrieve the scale, with the name DSID_INTU.
The import .bat files remain the same, except that the Python script called must be the following:
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)You can download the script here
Copy the script to the same location as the scripts in the original article.
You can either replace the lines in the .bat files
python c:/testgpkg/update_geopackage_dsid.py “%output_geopackage%” “!file!”
by
python c:/testgpkg/update_geopackage_dsid_prp.py “%output_geopackage%” “!file!”
Either delete the update_geopackage_dsid.py script and rename update_geopackage_dsid_prp.py to update_geopackage_dsid.py, which saves you having to modify the .bat files.
Add a Python script to filter layers by purpose
To filter the display of all layers loaded in the QGis project, you can use the following script:
# 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")
You can download the script here.
To eliminate Purpose filters, use the same script as for scaling (unfilter.py).
The S57 database project with Geopackage has come to an end. We’re now working on setting up an equivalent procedure using a PostgreSQL/Postgis database. Help us bring this project to a successful conclusion
