The first part of the Collaborative Financing for the Integration of Marine Data in QGIS project has been completed, thanks to contributions from ORANGE Marine, Geoconceptos Uruguay, Janez Avzec and Gaetan Mas. Our warmest thanks to all of them. We hereby publish the final content of this part of the work.
Given the volume of the result, we’re publishing it in two parts: the first concerns the creation and management of the database, the second the S57 symbology in QGis.
The S57 format is complex and does not fully comply with GIS format standards. Its conversion with GDAL is therefore either complex or incomplete. In this article, we’ll go into detail about the complete transformation of S57 into Geopackage file tables.
Geometry types
S57 “layers” are object classes. For example, the various land zones are encoded in the LNDARE object class.
The definition of this object class is:
Geo object: Land area (LNDARE) (P,L,A)
Attributes: CONDTN OBJNAM NOBJNM STATUS INFORM NINFOM
While all LNDARE objects have the same attributes, the same cannot be said for the type of geometry. The information (P,L,A) indicates that points, lines and polygons can be found in this object class. Contrary to GIS standards, all three types of geometry coexist within the same object class.
Whatever format we choose to integrate into QGis, we’ll need to create one layer per geometry type. If we don’t, GDAL will create the layer type according to the first entity found during conversion. If it’s a point type, the layer created will be a point type, and the line and polygon entities will be ignored. Similarly, if the first entity is of line type, points and polygons will be ignored.
It should also be noted that the S57 format has no constraints on object classes: if there is no object class in the area covered by the S57 file, there will be nothing concerning it in the file (no empty layer). Similarly, if only one type of geometry is present, even though all three types are possible, there will be no trace of the other types of geometry.
The processing of an S57 file with ogr2ogr must therefore be broken down into three steps, one for each type of geometry. The following options allow you to process each S57 object class by selecting only one type of geometry:
-where « OGR_GEOMETRY=’POINT’ or OGR_GEOMETRY=’MULTIPOINT’ »
-where « OGR_GEOMETRY=’LINESTRING’ or OGR_GEOMETRY=’MULTILINESTRING’ »
-where « OGR_GEOMETRY=’POLYGON’ or OGR_GEOMETRY=’MULTIPOLYGON’ »
For certain types of driver, GDAL allows you to create prefixes in output tables. In this case, you could create all tables (points, lines, polygons) in a single geopackage, prefixing them with pt_, li_ and pl_, for example. The problem is that GDAL’s S57 driver does not allow this option. This means creating three separate geopackages, in which tables are created according to geometry type. Each geopackage will then contain a table with the same name but a different geometry. While this is the simplest solution, it is not the most elegant or practical. Here, we’ll use three import geopackages for ogr2ogr commands, into which we’ll import tables from an S57 file. We’ll call them PointsENC, LinesENC and PolysENC. We’ll also create a single geopackage called ENC for the entire database. For each ENC map, we’ll import its contents into the three import geopackages with ogr2ogr, then run Python scripts to update the ENC geopackage database with the new imported tables. We are obliged to use Python scripts instead of SQL queries and functions because the SQL version of SQLite, used for geopackages, is very limited and does not allow the creation of procedures or functions.
If you have a single S57 file to import, the process is as follows:
- Create an empty ENC geopackage.
- Load point object classes into pointsENC geopackage with ogr2ogr
- Load Lines object classes into LinesENC geopackage with ogr2ogr
- Load Polygon object classes into PolysENC geopackage with ogr2ogr
- Remove empty tables from the three geopackages and add map_enc and scale attributes.
- Cloning of tables from the three import geopackages into the ENC geopackage, prefixing pointENC tables with pt_, LinesENC with li_ and PolysENC with pl_.
If you have a batch of S57 files to load simultaneously, the process is as follows:
- Load Point object classes into geopackage pointsENC with ogr2ogr, remove empty tables and add enc file name and scale to all table attributes.
- Load Lines object classes into LinesENC geopackage with ogr2ogr, deleting empty tables and adding enc file name and scale to all table attributes.
- Load Polygon object classes into PolysENC geopackage with ogr2ogr, deleting empty tables and adding enc file name and scale to all table attributes.
- Update existing tables in the ENC geopackage from the three ET import geopackages.
- Cloning of tables from the three import geopackages that were not present in the ENC geopackage,
- Detection and deletion of any duplicates resulting from ENC geopackage table updates (optional)
Processing bathymetric soundings
Bathymetric soundings pose several problems during format conversion. The first is that the depth value is not contained in an attribute, but as Z in the geometry (XYZ). The second is that the sounds are not of the Point type, but of the Multipoint type. To obtain sound values directly in an attribute field, you need to add two parameters to the ogr2ogr command line:
-oo SPLIT_MULTIPOINT=ON -oo ADD_SOUNDG_DEPTH=ON
Retrieving “parent” identifiers
Some objects can be grouped by a non-geographic parent object. For example, fire sectors are coded with one sector per record in the “LIGHTS” table. The different sectors of the same fire have no particular attribute that could indicate that they correspond to the same entity. This information is contained in a “parent” record. To retrieve this identifier as an attribute from the LIGTHS table, you need to add an option to the command line:
-oo RETURN_LINKAGES=O
This option creates a name_rcid attribute that is common to all sectors of the same fire, but also creates a series of fields (name_rcnm,ORNT,USAG,MASK).
“List” type processing
In addition to the classic field types (integer, real, text), the S57 format uses a special type: lists of character strings or integers. These field types are found in PostgreSQL, but not in shapefiles and geopackages.
In the absence of specific options, there will be warning messages indicating that the format is not supported, and the content of these fields will be ignored. To avoid losing content, add the :
-mapFieldType StringList=String,IntegerList=String
This option takes a list {..,…} and generates a text string with a first number indicating the number of elements in the list, a colon, then the list elements separated by commas:
{3} -> ‘(1:3)’
{3,1} -> ‘(2:3,1)’
This inevitably complicates the processing of these fields, but the information is not lost. On the other hand, the GDAL driver presents a problem when it comes to handling the resulting String format.
First of all, there’s no problem with the attributes of S57 objects. The field in the geopackage table is of type ‘Text’, which allows the creation of a character string of any length.
However, for “parent identifiers”, the definition used by the driver is ‘Text()’, i.e. a character string with a maximum length. This length is, a priori, determined by the content of the first record processed. In other words, it’s only by chance that you have a field with the length needed to integrate all the values in the file.
You’ll then see messages like:
WARNING : Value of field ‘NAME_RCID’ has 371 characters, whereas maximum allowed is 30
The solution isn’t complicated, but it does involve a fair amount of work. Assume that for the classic display of ENC cards, this is not a problem and you can simply ignore these warning messages.
On the other hand, if you need more advanced processing including this information, here’s the only way I’ve found yet:
- Run the ogr2ogr command to create the import geopackage
- Locate the maximum value indicated for each field type
- For the output table(s) from which you wish to retrieve all information, open the QGis database manager
- Delete all table records (DELETE FROM table_name)
- Go to the Table->Modify Table… menu.
- Select the relevant field from the list of fields, then click on Edit a column
- Enter the desired value in the Length field, then click OK
- Do the same for all other fields of interest
- Delete all other tables in the import geopackage that do not pose a problem
- Re-execute the ogr2ogr import command.
Workflow for loading a single S57 file
Creating the ENC geopackage
We’re going to create an empty geopackage that will receive the imported tables, then successive updates when other S57 files are loaded.
In the QGis Explorer panel, right-click on Geopackages and select Create database
In the window that opens, enter the geopackage name, leave the default option for table name and select No geometry in Geometry type.
Click OK to have the new empty geopackage added to the QGis geopackage connections. It is advisable to delete the empty ENC table created when the geopackage was created. To do this, open the QGis database manager from the main menu, position yourself on the ENC table of the ENC.gpkg geopackage and right-click to select Delete…
ogr2ogr commands for creating import geopackages
With all this in mind, here are the three lines of ogr2ogr commands for creating tables in import geopackages:
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
o execute these command lines, simply open the OSGeo4W shell window:
The Shell window opens
Enter the command lines, taking care to enter the path and filename with the .000 extension.
The result will be a series of tables created in each import geopackage. On the other hand, some tables will be completely empty. Indeed, if the geometry type is possible for an object class, the table will be created, but if there is no occurrence in the processed S57 file, the table will have no record.
Using Python scripts for database management
We will use the QGis Python console for all database management operations.
To open it, go to Extensions -> Python console. The Python console panels open:
To execute a Python script, load the code into the Editor panel. You can either copy and paste, or use the Directory icon to point to the .py file. As indentation is crucial for Python scripts, you should use the second option.
Remove empty tables and add map name and scale
As we saw earlier, when executing ogr2ogr commands, empty tables are created for layers where the geometry requested in the command line is not present. We will therefore remove these empty tables and take advantage of the Python script to add two attributes to all tables: the name of the ENC map and its scale.
As far as the name is concerned, it may be useful in future database maintenance to be able to select the entities corresponding to an S57 source file. This information is not contained in any attribute of the S57 file. We therefore have to enter it manually.
The same applies to the map scale. Although an attribute is provided in S57 to enter the scale of the “paper” map, the CSCALE attribute in the M_CSCL table, it is very rare for this table and its attribute to be present in S57 files.
While the name of the map is of relative use, the same cannot be said of the scale, as this is an essential criterion in the search for duplicates. Indeed, when mixing entities from several S57 files, data from different scales will coexist more or less happily.
The following Python code deletes empty tables from an import schema and adds two attributes (enc_chart and scale) to all 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)
Be careful to respect the indentation of the Python code. To simplify the task, download the .py file directly from this link: delete_empty_tables_in_geopackages.py.
For execution, be sure to modify the script’s input parameters:
You will only need to enter the geopackage_paths the first time you run the script. On the other hand, the ENC map name and scale must be modified when importing each new S57 file.
At the end of this run, you’ll have the new attributes in each table
Workflow for loading several S57 files at once
The command lines in the previous paragraphs assume that you have a single S57 file to load. But you can load several files simultaneously. By using the -append -update options, you can loop through all the .000 files in a directory. You’ll then have all the files loaded in the import geopackages.
Unlike loading a single file, when loading a batch of files, the addition of the file name and scale must be integrated into the S57->geopackage transcoding loop. To do this, we’ll use a python script that adds the enc_chart and scale attributes to the tables created by the ogr2ogr command and populates them with the corresponding file name and a scale defined by the producer. The script also deletes any empty tables present in the import geopackage.
Thanks to Janez AVSEC for the idea of using the DSID table to retrieve the data compilation scale.
The .bat files must be run from an OSGeo4W window to benefit from the QGis Python environment.
The call to each .bat file is of the form:
.\fichier.bat “directory containing S57 files to be loaded” “name of import geopackage”
For example:
.\mass_load_s57_polys.bat “c:/enc_db/enc000” “c:/enc_db/polys.gpkg”
The .bat files for this task are as follows:
To import point type tables:
@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
You’ll notice that, in addition to the ogr2ogr command for importing point layers, we import two tables without geometry: DSID and C_AGGR. The DSID table contains two attributes: the S57 file name and the data compilation scale. This table will then be used by the Python script (update_geopackage_dsid.py) to fill each geopackage table with the original S57 file and compilation scale.
To import line-type tables:
@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é.
pauseTo import polygon tables:
@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
The .bat files contain the ogr2ogr command and a line to execute a Python script, update_geopackage_dsid.py. Be sure to change the path of this script in the .bat file to the one of your choice. This script adds the two attributes enc_chart and scale to all geopackage tables.
The python code used is as follows:
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)You could create a single .bat file for all three treatments, but you’ll need to check the results for each type of geometry.
You can download the .bat files and the Python code by clicking here.
Continues regardless of the number of files loaded
Cloning and/or adding imported tables to the ENC geopackage
For the first S57 load, the following Python script creates all the tables present in the three import geopackages. For subsequent loads, if the table already exists in the ENC geopackage, the records of the table in the import geopackage are added to the records already present in ENC. If, on the other hand, the table does not yet exist, it is created and populated with records from the import geopackage. In the ENC geopackage file, tables are prefixed with pt_, li_ and pl_ to indicate their geometry type.
Cloning is performed using the following code:
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()
You can download the .py file here: clone_or_append_tables_with_prefix.py.
Duplicate management
When two ENC maps overlap, the entities present in the overlap zone will be inserted as many times as there are overlaps. To avoid duplicates and more, the database must be cleaned after each import of an S57 file.
But there are two different types of duplicates. On the one hand, overlapping areas between two maps of similar scale will produce “real” duplicates that need to be cleaned up. On the other hand, overlapping maps of very different scales need to be treated differently. The less detailed map (small scale) will have a selection of data from the more detailed map (large scale). In the latter case, it is not advisable to remove this information.
The Python code below will remove the “true” duplicates. For other duplicates, we recommend working with QGis filters. Depending on your needs, you can use the “scale” attribute we added when deleting the empty tables to define which data will or will not be displayed on your map.
In this example, only entities corresponding to a scale between 12500 and 45000 will be displayed on your map.
This operation corresponds to the filtering of a layer, but… what do you do when you have 130 layers displayed? To save you having to perform this operation 130 times, here’s a Python code that lets you filter all the layers in your project:
#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'")
And another to remove the filter from all displayed layers:
# 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()}")
Download the python codes here.
The following code creates a temporary table (temp_table) for each table in the ENC geopackage. It then fills this table with the separate LNAM-based records. It then deletes the original table and renames the temporary table with the name of the original table.
This treatment applies to all S57 tables, except:
- fire (LIGHTS) and bathymetric sounding (SOUNDG) tables, where the LNAM attribute is not useful. In fact, each fire sector has the same LNAM and name_rcid, but corresponds to a different color and geographic sector. Similarly, the LNAM and name_rcid of bathymetric soundings are the same for all sounding points in the original multipoint entity. In the first case, we need to search for duplicate colors and sectors for the same name_rcid, and for probes we need to find probes with the same XY in their geometry.
- SEAARE, LNDRGN and ADMARE tables, which cover large areas and are split according to the map footprint, while keeping the same LNAM. In this case, even if the LNAM is the same, the polygons do not have the same shape and area.
- DSID and C_AGGR tables, which have no LNAM and, by principle, cannot have duplicates.
from osgeo import ogrdef 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()
You can download the .py file from this link: supprime_doublons.py. Don’t forget to change the paths and filenames to match your own. Note also the line :
if table_name not in [“pl_SEAARE”, “pl_LNDRGN”, “pl_ADMARE”, “layer_styles”, “natsurf”,”DSID”,”C_AGGR”]:
In principle, you’ll only have the prefixed tables and the layer_styles and natsurf tables. This line prevents pl_SEAARE, pl_LNDRGN and pl_ADMARE from being processed, but more importantly, it prevents the layer_styles and natsurf tables we use for symbologies from being emptied. If you need to create tables other than those from the s57 file in the geopackage, you need to take into account the names of these tables in this script, by adding them to this line of exceptions.
The S57 database project with Geopackage has come to an end. We are now in the process of setting up an equivalent procedure using a PostgreSQL/Postgis database. Help us bring this project to a successful conclusion!
