As the progress of this project has been the subject of several articles and updates, the following link provides you with a summary document (PDF)
This will give you a complete and coherent overview of all the steps involved in creating and loading an ENC database in QGis in geopackage format and all the elements required (.bat, .py, .svg) to set up an appropriate automatic symbology.
The project to integrate ENC maps into QGis in the form of Geopackage is now complete. This article completes the two previously published articles on setting up and managing the Geopackage database.
In the first article: Creating a Geopackage database for ENC maps (Part 1: Setting up the database), you’ll find the essentials of working with Geopackages to integrate data from S57 (ENC) files. You should start with this article to understand the procedure for importing S57 files with ogr2ogr, as well as the various steps involved in formatting the data. You’ll also be able to download all the tools (.bat command files, python scripts) needed to complete these steps.
Following publication of this article, we continued to work on an equivalent procedure, but using a PostgreSQL/postgis database. This enabled us to find several improvements to the initial workflow.
A second article: Update(1) :ENC database with Geopackage under QGis, presents an update of the import procedure (.bat files with ogr2ogr commands) to enable easy management of duplicate information. Such duplication occurs when maps at different scales are integrated into the same database. The addition of a “purpose” attribute to the ENC map allows only unique information to be displayed for each map purpose.
As the aim of the project was to produce in QGis a symbology equivalent to that of paper maps, there was one point left to resolve: the “quality” parameter of the position of geographical objects. This is the main point addressed in this article.
The three articles are complementary in terms of content, but for the files made available, please download only the files proposed in this article. They end with ‘V3’ and ensure that you’re using the latest version.
The S57 format and GIS
If you’re about to start working with S57 files in a GIS, there are a few things you need to be aware of to navigate without pitfalls.
First of all, the structure of S57 files does not correspond to the structures adopted in GIS.
In a GIS, you have a geographic object represented by a table with two types of information: the geometry of the object’s entities and the attributes of these entities.
If you have other objects with identical geometries, the geometric information is duplicated, once in each table.
For the S57 format, the main objective is to optimize information storage and therefore not to duplicate information. If an object has a point entity, a point will be created. If other objects have entities located on this point, we use the reference of the point that has already been created. In this way, a point is described only once in the file. The same applies to polylines and surfaces. An S57 file will therefore have a series of tables containing geometric information, called “primitives”:
- IsolatedNode (points)
- ConnectedNode (points)
- Edge (polylines)
- Face (polygons)
The attribute table for the various S57 objects contains only the object attributes.
What complicates the task is that there are two attributes that refer to geometries: posacc (the estimation of positional accuracy, a quantitative value) and quapos (positional quality, a qualitative variable).
These two attributes can be found in the primitive tables.
To switch from the S57 structure to a GIS structure (shapefile, geopackage, postgis) we use the GDAL library and its ogr2ogr command.
This command will create GIS tables from the S57 structure, creating one table per S57 object, assigning the corresponding geometries from the primitive tables to each entity and adding the attributes of the S57 object to each entity. The trace of the primitives used for the geometry of each entity can be found in the NAME_RCID field of the GIS tables, provided the options -oo “RETURN_LINKAGES=ON” -oo “LNAM_REFS=ON” have been added to the ogr2ogr command line.
The following figure shows a point-type layer. The value indicated in the NAME_RCID field is the RCID of the point used in the IsolatedNode table.
The following figure shows an example of a linear layer. The values indicated in the NAME_RCID field are the RCID values of the polylines used in the Edge table.
The following figure shows an example of a polygon layer. The values indicated in the NAME_RCID field are those of the polylines used in the Edge table.
In order to retrieve the QUAPOS and POSACC attributes of each entity in the point-type tables, we need to retrieve the values of the IsolatedNode point and assign them to the tables of the various ENC objects.
If the identifiers were directly the RCIDs in the ENC tables, we could make a join between each table (soundg, obstrn,…) and IsolatedNode. But as you can see in the previous images, the NAME_RCID attribute is of stringlist type, which blocks this solution. So we’ve developed a python script that does the job when loading S57->geopackage data.
Downloading V3 files
You can download all the .bat and .py files mentioned in this article from this link.
Modification of ogr2ogr commands for importing S57 files
Importing “point” layers
To import primitives from S57 files, the -oo RETURN_PRIMITIVES = ON option must be added to the import commands. In the case of points, it is also necessary to add the option -nlt MULTIPOINT, otherwise the IsolatedNode table will be created as a point and multipoint layers will not be supported.
The .bat file for importing point-type layers will then be :
@echo off
setlocal enabledelayedexpansionREM Assurez-vous qu'au moins 2 arguments ont été fournis
if "%~2"=="" (
echo Usage: %0 directory000 output_geopackage
exit /b 1
)REM Récupérer les arguments de la ligne de commande
set "directory=%~1"
set "output_geopackage=%~2"REM Compter le nombre total de fichiers à traiter
set /a total_files=0
for /r "%directory%" %%a in (*.000) do (
set /a total_files+=1
)REM Initialiser le compteur de fichiers traités
set /a processed_files=0
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 /a processed_files+=1
echo Traitement en cours : !processed_files! sur !total_files!
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 -oo RETURN_PRIMITIVES=ON -nlt MULTIPOINT -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:/testgpkgV3/update_geopackage_dsid_V3.py "%output_geopackage%" "!file!"
)
echo Import terminé. Ajout de posass et quapos.
python c:/testgpkgV3/add_posacc_quapos_to_pointsV3.py "%output_geopackage%"
echo Traitement terminé.
pauseTo run it, open an OSGeo4W Shell window, go to the directory containing the .bat files and enter the command:
.\mass_load_s57_points_V3.bat files_path_000 chemin/pointsENC.gpkg
First, edit the file to modify the lines:
python c:/testgpkgV3/update_geopackage_dsid_prp_prim.py “%output_geopackage%” “!file!”
and
python c:/testgpkgV3/add_posacc_quapos_to_pointsV3.py
To enter the paths you’ve chosen to store your python scripts.
The first script, update_geopackage_dsid_V3.py, has been modified to add, in addition to enc_chart, scale and purpose, the POSACC and QUAPOS attributes to all point tables, at the time of the S57->geopackage switch.
The second script, add_posacc_quapos_to_pointsV3.py, once all S57 files in the ENC directory have been loaded, updates the values of these attributes for each entity in the points tables.
import sys
from osgeo import ogr
def add_posacc_quapos_to_pointsV3 (geopackage_path):
# Chemin vers le GeoPackage
# Liste des tables à exclure
tables_exclues = ["IsolatedNode", "ConnectedNode","DSID","C_AGGR","C_ASSO","layer_styles"]
# Ouvrir le GeoPackage
driver = ogr.GetDriverByName("GPKG")
geopackage = driver.Open(geopackage_path, 1)
# Récupérer la table IsolatedNode
isolated_node_table = geopackage.GetLayerByName("IsolatedNode")
# Parcourir toutes les tables du GeoPackage
for i in range(geopackage.GetLayerCount()):
table = geopackage.GetLayerByIndex(i)
table_name = table.GetName()
# Vérifier si la table n'est pas dans la liste des tables exclues et si elle n'est pas IsolatedNode
if table_name not in tables_exclues and table_name != "IsolatedNode" :
print(f"Traitement table {table_name}")
# Récupérer les noms des champs des tables
table_defn = table.GetLayerDefn()
field_names = [table_defn.GetFieldDefn(j).GetName() for j in range(table_defn.GetFieldCount())]
# Index des champs rcid et enc_chart dans la table IsolatedNode
isolated_node_defn = isolated_node_table.GetLayerDefn()
rcid_index_isolated_node = isolated_node_defn.GetFieldIndex("RCID")
enc_chart_index_isolated_node = isolated_node_defn.GetFieldIndex("enc_chart")
# Parcourir les enregistrements de la table
table.ResetReading()
for feature in table:
rcid_full = feature.GetField("NAME_RCID")
rcid = rcid_full.split(":")[1]
rcid = rcid[:-1]
enc_chart = feature.GetField("enc_chart")
# Rechercher le rcid et enc_chart dans la table IsolatedNode
isolated_node_table.SetAttributeFilter("RCID = '{}' AND enc_chart = '{}'".format(rcid, enc_chart))
isolated_node_feature = isolated_node_table.GetNextFeature()
if isolated_node_feature:
# Mettre à jour les attributs POSACC et QUAPOS de la table en cours avec ceux d'IsolatedNode
posacc_index = table_defn.GetFieldIndex("POSACC")
quapos_index = table_defn.GetFieldIndex("QUAPOS")
if posacc_index >= 0:
posacc = isolated_node_feature.GetField("POSACC")
feature.SetField(posacc_index, posacc)
if quapos_index >= 0:
quapos = isolated_node_feature.GetField("QUAPOS")
feature.SetField(quapos_index, quapos)
# Sauvegarder les modifications
table.SetFeature(feature)
# Fermer le GeoPackage
geopackage = None
if __name__ == "__main__":
if len(sys.argv) != 2:
print("Usage: python script.py geopackage_path ")
sys.exit(1)
geopackage_path = sys.argv[1]
add_posacc_quapos_to_pointsV3 (geopackage_path)
In the following example, we have loaded 7 ENC maps into the pointsENC.gpkg import geopackage.
You can see that the geopackage tables contain the POSACC and QUAPOS attributes and that the last script has filled in the corresponding values for these fields.
In the second part of this final version, dedicated to symbology, you’ll see the new management of the “nature of bottom” layer. This symbology is based on the creation of a “Label” atribute in the SBDARE layer and its filling with a python script Label_sbdare.py . This operation can be performed on the complete database or, here, on the pointsENC.gpkg geopackage. As it takes some time to calculate, we advise you to do it at this stage. Please refer to the article ENC Geopackage maps in QGis Final version: part 2 for instructions.
Importing line layers
To import primitives from S57 files, you need to add the -oo RETURN_PRIMITIVES = ON option to the import commands.
The line layer import .bat file will then be :
@echo off
setlocal enabledelayedexpansion
REM Assurez-vous qu'au moins 2 arguments ont été fournis
if "%~2"=="" (
echo Usage: %0 directory000 output_geopackage
exit /b 1
)
REM Récupérer les arguments de la ligne de commande
set "directory=%~1"
set "output_geopackage=%~2"
REM Compter le nombre total de fichiers à traiter
set /a total_files=0
for /r "%directory%" %%a in (*.000) do (
set /a total_files+=1
)
REM Initialiser le compteur de fichiers traités
set /a processed_files=0
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 /a processed_files+=1
echo Traitement en cours : !processed_files! sur !total_files!
set "file=%%~ni"
ogr2ogr -f GPKG -skipfailures -append -update -where "OGR_GEOMETRY='LINESTRING' or OGR_GEOMETRY='MULTILINESTRING'" -oo RETURN_PRIMITIVES=ON -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:/testgpkgV3/update_geopackage_dsid_V3.py "%output_geopackage%" "!file!"
)
echo Import terminé. Ajout de posacc et quapos.
python c:/testgpkgV3/add_posacc_quapos_to_linesV3.py "%output_geopackage%"
echo Traitement terminé.
pause
To run it, open an OSGeo4W Shell window, go to the directory containing the .bat files and enter the command:
.\mass_load_s57_lines_V3.bat répertoire_fichiers_000 chemin/pointsENC.gpkg
First, edit the file to modify the lines:
python c:/testgpkgV3/update_geopackage_dsid_prp_prim.py “%output_geopackage%” “!file!”
and
python c:/testgpkgV3/add_posacc_quapos_to_linesV3.py
to enter the paths you’ve chosen to store your python scripts.
The second script, add_posacc_quapos_to_linesV3.py, once all the S57 files in the ENC directory have been loaded, updates the values of these attributes for each entity in the line tables.
import sys
from osgeo import ogr
def add_posacc_quapos_to_linesV3 (geopackage_path):
# Liste des tables à exclure
tables_exclues = ["Edge", "ConnectedNode","DSID","C_AGGR","C_ASSO","Edge","FACE"] # Ajoutez les noms des tables que vous souhaitez exclure
# Ouvrir le GeoPackage
driver = ogr.GetDriverByName("GPKG")
geopackage = driver.Open(geopackage_path, 1)
# Récupérer la table Edge
isolated_node_table = geopackage.GetLayerByName("Edge")
# Parcourir toutes les tables du GeoPackage
for i in range(geopackage.GetLayerCount()):
table = geopackage.GetLayerByIndex(i)
table_name = table.GetName()
# Vérifier si la table n'est pas dans la liste des tables exclues et si elle n'est pas Edge
if table_name not in tables_exclues and table_name != "Edge":
print(f"Traitement table {table.name()}")
# Récupérer les noms des champs des tables
table_defn = table.GetLayerDefn()
field_names = [table_defn.GetFieldDefn(j).GetName() for j in range(table_defn.GetFieldCount())]
# Index des champs rcid et enc_chart dans la table Edge
isolated_node_defn = isolated_node_table.GetLayerDefn()
rcid_index_isolated_node = isolated_node_defn.GetFieldIndex("RCID")
enc_chart_index_isolated_node = isolated_node_defn.GetFieldIndex("enc_chart")
# Parcourir les enregistrements de la table
table.ResetReading()
for feature in table:
rcid_full = feature.GetField("NAME_RCID")
rcid = rcid_full.split(":")[1]
rcid = rcid[:-1]
enc_chart = feature.GetField("enc_chart")
# Rechercher le rcid et enc_chart dans la table Edge
isolated_node_table.SetAttributeFilter("RCID = '{}' AND enc_chart = '{}'".format(rcid, enc_chart))
isolated_node_feature = isolated_node_table.GetNextFeature()
if isolated_node_feature:
# Mettre à jour les attributs POSACC et QUAPOS de la table en cours avec ceux d'Edge
posacc_index = table_defn.GetFieldIndex("POSACC")
quapos_index = table_defn.GetFieldIndex("QUAPOS")
if posacc_index >= 0:
posacc = isolated_node_feature.GetField("POSACC")
feature.SetField(posacc_index, posacc)
if quapos_index >= 0:
quapos = isolated_node_feature.GetField("QUAPOS")
feature.SetField(quapos_index, quapos)
# Sauvegarder les modifications
table.SetFeature(feature)
# Fermer le GeoPackage
geopackage = None
if __name__ == "__main__":
if len(sys.argv) != 2:
print("Usage: python script.py geopackage_path ")
sys.exit(1)
geopackage_path = sys.argv[1]
add_posacc_quapos_to_linesV3 (geopackage_path)
In the case of point tables, the link between the NAME_RCID of the tables and the RCID of the IsolatedNode table poses no problem. In the case of polylines, on the other hand, the list of polylines for a given geographic entity concerns several entities in the Edge layer. We wondered whether, for an entity such as DEPCNT, it was possible for the polylines making it up to have different QUAPOS values.
We carried out the test on a database containing 650 ENC maps at all scales. The DEPCNT layer contained 162,585 features. Only 38,731 had RCIDs in the NAME_RCID list of the Edges table with a non-zero QUAPOS value. We found 193 DEPCNT entities that had different QUAPOS values in their NAME_RCID list. Since this represents 0.1%, we decided not to take this possibility into account, as the complexity of polyline deconstruction is extreme.
It turns out that the Edge table resulting from the ogr2ogr driver poses no problem at attribute level (the join between NAME_RCID and RCID is good), but that the table’s geometries contain almost half invalid geometries (lines with a single point or null geometries). Yet the polylines of S57 entities have no errors.
It would therefore be necessary to repair and rebuild all Edge geometries in order to be able to deconstruct DEPCNT geometries so as not to have different QUAPOS in any case.
In the following example, we have loaded 7 ENC maps into the import geopackage linesENC.gpkg
Vous pouvez constater que les tables du geopackage contiennent les attributs POSACC et QUAPOS et que le dernier script a bien renseigné les valeurs correspondantes de ces champs.
Importing polygon layers
Although the S57 format allows polygons as primitives, this option is not used. It is therefore not necessary to modify the import of polygon data. The primitives (Face) are requested in the ogr2ogr command line, and the POSACC and QUAPOS fields are added to the tables for later use, but there is no need to add a script to fill in the values.
Cloning and/or adding imported tables to the ENC geopackage
For the first S57 load, the following Python script creates all 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.
Compared with V2, the changes made to the script are as follows:
- SRC information for geopackage tables on creation. The EPSG:4326 code is assigned to geometries.
# Créer une nouvelle couche de destination si elle n'existe pas
output_srs = osr.SpatialReference()
output_srs.ImportFromEPSG(4326) # Définir EPSG 4326
output_table = output_gpkg.CreateLayer(output_table_name, geom_type=input_table.GetGeomType(), options=["OVERWRITE=YES"], srs=output_srs)
- Force spatial index creation on all geopackage tables
# Créer l'index spatial
output_gpkg.ExecuteSQL(f"CREATE SPATIAL INDEX ON {output_table_name}")Cloning is performed with the python code: clone_or_append_tables_with_prefixV3.py
Be sure to modify the line indicating the names of the geopackages you created with the previous .bat files:
input_geopackages = [“c:/testgpkg/pointsENC.gpkg”, “c:/testgpkg/LinesENC.gpkg”, “c:/testgpkg/PolysENC.gpkg”]
with the import geopackages you’ve chosen, and the line
output_gpkg = ogr.Open(“c:/testgpkg/ENC.gpkg”, 1)
with the geopackage database you’ve created as described in Creating a Geopackage database for ENC maps (Part 1: Creating the database).
Once this script has been run, your database will be updated with the new S57 files.
Updating a geopackage V1 or V2 database to version 3
The procedure described above allows you to create a new geopackage database. But if you already have a geopackage database, there’s no need to go through the whole process again. Here you’ll find the procedure and tools you need to update the structure (POSAC and QUAPOSQ) and set the attribute values for existing tables.
Download the files for the procedure here.
The various stages of this “upgrade” are:
1- From a directory containing all the .000 files that have already been loaded into the database, execute a .bat file, upgrade_V2_to_V3.bat, which creates an update.gpkg geopackage containing the primitives (IsolatedNode, ConnectedNode,Edge) of the S57 files.
.\ upgrade_V2_to_V3.ba répertoire_fichiers_000 chemin/update.gpkg
2- Run append_tables_upgrade.py in the QGis Python console. Be sure to modify the path for the update.gpkg file
input_geopackages = [“c:/testgpkgV2/update.gpkg”]
and for your ENC database
output_gpkg = ogr.Open(“c:/enc_db/ENC4.gpkg”, 1)
This script adds the primitives to the existing database.
3- In the QGis Python console, run the add_fields_upgrade.py script. Be sure to change the path to your database’s .gpkg file
geopackage = ogr.Open(“c:/enc_db/ENC4.gpkg”, 1)
This script creates the POSACC and QUAPOS attributes in all database tables.
4-In the QGis Python console, run the script add_posacc_quapos_upgrade_points.py. Be sure to change the path to your database’s .gpkg file and add tables to exclude if you’ve created other tables in the database.
#Chemin vers le GeoPackage
geopackage_path = “c:/enc_db/ENC4.gpkg”
#Liste des tables à exclure
tables_exclues = [“IsolatedNode”, “ConnectedNode”,”DSID”,”C_AGGR”,”C_ASSO”,”layer_styles”] # Ajoutez les noms des tables que vous souhaitez exclure
This script will populate the POSACC and QUAPOS fields in all Point type tables with the values contained in the IsolatedNode primitive.
5-In the QGis Python console, run the script add_posacc_quapos_upgrade_lines.py. Be sure to change the path to your database’s .gpkg file and to add tables to exclude if you have created other tables in the database.
#Chemin vers le GeoPackage
geopackage_path = “c:/enc_db/ENC4.gpkg”
#Liste des tables à exclure
tables_exclues = [“IsolatedNode”, “ConnectedNode”,”DSID”,”C_AGGR”,”C_ASSO”,”layer_styles”] # Ajoutez les noms des tables que vous souhaitez exclure
This script will populate the POSACC and QUAPOS fields in all Line type tables with the values contained in the Edge primitive.
Thanks for the helpful tutorial!
Is there any chance you can provide additional bash scrips for the batch scripts or maybe port all the scripts to python?