Como el progreso del proyecto ha sido objeto de varios artículos y actualizaciones, ponemos a su disposición un documento resumen (PDF) en el siguiente enlace
En él podréis encontrar una visión completa y coherente de todos los pasos necesarios para crear y cargar una base de datos ENC en QGis en formato geopackage y todos los elementos necesarios (.bat, .py, .svg) para configurar una simbología automática adecuada.
El proyecto de integración de los mapas ENC en QGis en forma de Geopackage ha finalizado. Este artículo completa los dos artículos publicados anteriormente sobre la configuración y gestión de la base de datos Geopackage.
En el primer artículo: Creación de una base de datos Geopackage para mapas ENC (parte 1: configuración de la base de datos), encontrará los aspectos esenciales del trabajo con geopackages para integrar datos de archivos S57 (ENC). Debería empezar por este artículo para comprender el procedimiento de importación de archivos S57 con ogr2ogr, así como los distintos pasos necesarios para formatear los datos. También podrá descargar todas las herramientas (archivos de comandos .bat, scripts python) necesarias para completar estos pasos.
Tras la publicación de este artículo, seguimos trabajando en un procedimiento equivalente pero utilizando una base de datos PostgreSQL/postgis. Esto nos permitió encontrar varias mejoras en el flujo de trabajo inicial.
Un segundo artículo, Actualización(1) :Base de datos ENC con Geopackage en QGis, actualiza el procedimiento de importación (archivos .bat con comandos ogr2ogr) para facilitar la gestión de la información duplicada. Estas duplicaciones se producen cuando se integran en la misma base de datos mapas a diferentes escalas. La adición de un atributo «propósito» al mapa ENC significa que sólo puede mostrarse información única para cada propósito del mapa.
Como el objetivo del proyecto era producir una simbología en QGis equivalente a la de los mapas en papel, quedaba un punto por resolver: tener en cuenta el parámetro «calidad» de la posición de los objetos geográficos. Este es el principal punto abordado en el presente artículo.
Los tres artículos son complementarios en cuanto a contenido, pero en cuanto a los ficheros puestos a disposición, descargue únicamente los ficheros sugeridos en este artículo. Terminan con ‘V3’ y asegúrese de que está utilizando la última versión.
El formato S57 y el SIG
Si va a empezar a trabajar con ficheros S57 en un SIG, hay algunas cosas que debe tener en cuenta para navegar sin problemas.
En primer lugar, la estructura de los ficheros S57 no se corresponde con las estructuras adoptadas en los SIG.
En un SIG tienes un objeto geográfico que está representado por una tabla con dos tipos de información: la geometría de las entidades del objeto y los atributos de estas entidades.
Si tiene otros objetos con geometrías idénticas, la información geométrica se duplica, una vez en cada tabla.
En el formato S57, el objetivo principal es optimizar el almacenamiento de la información y, por tanto, no duplicarla. Si un objeto tiene una entidad puntual, se creará un punto. Si otros objetos tienen entidades situadas en este punto, se utilizará la referencia del punto ya creado. De este modo, un punto sólo se describe una vez en el fichero. Lo mismo ocurre con las polilíneas y las superficies. Por tanto, un fichero S57 tendrá una serie de tablas que contienen información geométrica, conocidas como «primitivas»:
- IsolatedNode (puntos)
- ConnectedNode (puntos)
- Edge (polilíneas)
- Face (polígonos)
La tabla de atributos de los distintos objetos S57 sólo contiene los atributos de los objetos.
Lo que complica la tarea es que hay dos atributos que se refieren a las geometrías: posacc (la precisión estimada de la posición, un valor cuantitativo) y quapos (la calidad de la posición, una variable cualitativa).
Estos dos atributos se encuentran en las tablas primitivas.
Para pasar de la estructura S57 a una estructura SIG (shapefile, geopackage, postgis) utilizamos la biblioteca GDAL y su comando ogr2ogr.
Este comando creará tablas SIG a partir de la estructura S57, creando una tabla por cada objeto S57, asignando la geometría correspondiente de las tablas de primitivas a cada entidad y añadiendo los atributos del objeto S57 a cada entidad. La traza de las primitivas utilizadas para la geometría de cada entidad puede encontrarse en el campo NAME_RCID de las tablas SIG, siempre que se hayan añadido las opciones -oo «RETURN_LINKAGES=ON» -oo «LNAM_REFS=ON» a la línea de comandos de ogr2ogr.
La siguiente figura muestra una capa de tipo punto. El valor indicado en el campo NAME_RCID es el RCID del punto utilizado en la tabla IsolatedNode.
La siguiente figura muestra un ejemplo de capa de tipo lineal. Los valores indicados en el campo NAME_RCID son los de los RCID de las polilíneas utilizadas en la tabla Edge.
La siguiente figura muestra un ejemplo de capa de tipo poligonal. Los valores indicados en el campo NAME_RCID son los de las polilíneas utilizadas en la tabla Edge.
Para recuperar los atributos QUAPOS y POSACC de cada entidad en las tablas de tipo punto, necesitamos recuperar los valores del punto IsolatedNode y asignarlos a las tablas de los distintos objetos ENC.
Si los identificadores fueran directamente los RCID de las tablas ENC, podríamos hacer un join entre cada tabla (soundg, obstrn,…) y IsolatedNode. Pero como se puede ver en las imágenes anteriores, el atributo NAME_RCID es de tipo stringlist, lo que bloquea esta solución. Por ello hemos desarrollado un script en python que hace el trabajo cuando se cargan los datos de S57->geopackage.
Descarga de archivos V3
Puede descargar todos los archivos .bat y .py mencionados en este artículo desde este enlace.
Modificación de los comandos de ogr2ogr para importar archivos S57
Importación de capas «puntuales
Para importar primitivas desde ficheros S57, es necesario añadir la opción -oo RETURN_PRIMITIVES = ON a los comandos de importación. En el caso de los puntos, también es necesario añadir la opción -nlt MULTIPOINT, ya que de lo contrario la tabla IsolatedNode se creará como un punto y las capas multipunto no serán compatibles.
El archivo .bat para importar capas de tipo punto será entonces :
@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é.
pausePara ejecutarlo, abra una ventana OSGeo4W Shell, vaya al directorio que contiene los archivos .bat e introduzca el comando:
.\mass_load_s57_points_V3.bat directorio_archivos_000 ruta/pointsENC.gpkg
En primer lugar, editar el archivo para cambiar las líneas:
python c:/testgpkgV3/update_geopackage_dsid_prp_prim.py «%output_geopackage%» «!file!»
y
python c:/testgpkgV3/add_posacc_quapos_to_pointsV3.py
Para introducir las rutas que ha elegido para almacenar los scripts de python.
El primer script, update_geopackage_dsid_V3.py, se ha modificado para añadir los atributos POSACC y QUAPOS a todas las tablas de puntos, además de enc_chart, scale y purpose, cuando se realiza el cambio S57->geopackage.
El segundo script, add_posacc_quapos_to_pointsV3.py, una vez cargados todos los ficheros S57 en el directorio ENC, actualiza los valores de estos atributos para cada entidad en las tablas de puntos.
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)
En el siguiente ejemplo, hemos cargado 7 mapas ENC en el geopackage import pointsENC.gpkg.
Puedes ver que las tablas del geopackage contienen los atributos POSACC y QUAPOS y que el último script ha rellenado los valores correspondientes para estos campos.
En la segunda parte de esta versión final, dedicada a la simbología, verás la nueva gestión de la capa «naturaleza del fondo marino». Esta simbología se basa en la creación de un atributo «Label» en la capa SBDARE y su relleno con un script python Label_sbdare.py . Esta operación puede realizarse sobre la base de datos completa o aquí, sobre el geopackage pointsENC.gpkg. Como se tarda algún tiempo en calcular, le aconsejamos que lo haga en esta fase. Consulte el artículo Mapas del geopackage ENC en QGis : versión final parte 2 para obtener instrucciones.
Importar capas de ‘línea
Para importar primitivas desde archivos S57, necesita añadir la opción -oo RETURN_PRIMITIVES = ON a los comandos de importación.
El archivo .bat de importación de capas de línea será entonces :
@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
Para ejecutarlo, abra una ventana OSGeo4W Shell, vaya al directorio que contiene los archivos .bat e introduzca el comando:
.\mass_load_s57_lines_V3.bat directorio_archivos_000 ruta/puntosENC.gpkg
En primer lugar, editar el archivo para cambiar las líneas:
python c:/testgpkgV3/update_geopackage_dsid_prp_prim.py «%output_geopackage%» «!file!»
y
python c:/testgpkgV3/add_posacc_quapos_to_linesV3.py
para introducir las rutas que ha elegido para almacenar los scripts python.
El segundo script, add_posacc_quapos_to_linesV3.py, una vez cargados todos los ficheros S57 en el directorio ENC, actualiza los valores de estos atributos para cada entidad en las tablas de líneas.
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)
En el caso de las tablas de tipo punto, el vínculo entre el NAME_RCID de las tablas y el RCID de la tabla IsolatedNode no plantea ninguna duda. Sin embargo, en el caso de las polilíneas, para una entidad geográfica, la lista de polilíneas concierne a varias entidades de la capa Edge. Nos preguntamos si para una característica, por ejemplo una característica DEPCNT, era posible que las polilíneas que la componen tuvieran valores de QUAPOS diferentes.
Realizamos la prueba en una base de datos que contenía 650 mapas ENC a todas las escalas. La capa DEPCNT contenía 162.585 elementos. Sólo 38.731 tenían RCID en la lista NAME_RCID de la tabla Edges con un valor QUAPOS distinto de cero. Encontramos 193 entidades DEPCNT que tenían diferentes valores QUAPOS en su lista NAME_RCID. Dado que esto representa el 0,1%, decidimos no tener en cuenta esta posibilidad, ya que la complejidad de la deconstrucción de polilíneas es extrema.
Resulta que la tabla Edge resultante del controlador ogr2ogr no plantea ningún problema a nivel de atributos (la unión entre NAME_RCID y RCID es buena), pero que las geometrías de la tabla contienen casi la mitad de geometrías no válidas (líneas con un solo punto o geometrías nulas). Sin embargo, las polilíneas de las entidades S57 no presentan errores.
Por tanto, sería necesario reparar y reconstruir todas las geometrías Edge para poder deconstruir las geometrías DEPCNT y no tener QUAPOS diferentes en cada caso.
En el siguiente ejemplo, hemos cargado 7 mapas ENC en el geopackage de importación linesENC.gpkg
Importación de capas poligonales
Aunque el formato S57 permite utilizar polígonos como primitivas, esta opción no se utiliza. Por lo tanto, no es necesario modificar la importación de datos poligonales. Las primitivas (Face) se solicitan en la línea de comandos ogr2ogr, y los campos POSACC y QUAPOS se añaden a las tablas para su posible uso futuro, pero no es necesario añadir un script para rellenar los valores.
Clonación y/o adición de tablas importadas al geopackage ENC
Para la primera carga del S57, el siguiente script de Python crea todas las tablas en los tres geopacks de importación. Para las cargas posteriores, si la tabla ya existe en el geopackage ENC, los registros de la tabla en el geopackage de importación se añaden a los registros ya presentes en ENC. Si, por el contrario, la tabla aún no existe, se crea y se rellena con los registros del geopackage de importación. En el archivo del geopackage ENC, las tablas llevan los prefijos pt_, li_ y pl_ para indicar su tipo de geometría.
Se han introducido los siguientes cambios en el script en comparación con la V2:
- El SRC de las tablas del geopackage se rellena durante la creación. Se asigna el código EPSG:4326 a las geometrías.
# 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)
- Forzar la creación del índice espacial en todas las tablas de paquetes geográficos
# Créer l'index spatial
output_gpkg.ExecuteSQL(f"CREATE SPATIAL INDEX ON {output_table_name}")La clonación se realiza mediante el código python: clone_or_append_tables_with_prefixV3.py
Asegúrate de cambiar la línea que indica los nombres de los geopackages que creaste con los archivos .bat anteriores:
input_geopackages = [» c:/testgpkg/pointsENC.gpkg «, » c:/testgpkg/LinesENC.gpkg «, » c:/testgpkg/PolysENC.gpkg «].
con los geopackages de importación que hayas elegido, y la línea
output_gpkg = ogr.Open(«c:/testgpkg/ENC.gpkg», 1)
con el geopackage de base de datos que ha creado siguiendo las instrucciones en Creación de una base de datos Geopackage para mapas ENC (parte 1: creación de la base de datos).
Una vez ejecutado este script, su base de datos se actualizará con los nuevos archivos S57.
Actualización de una base de datos de geopacks V1 o V2 a la versión 3
El procedimiento que acabamos de describir puede utilizarse para crear una nueva base de datos de geopacks. Pero si ya tienes una base de datos creada, no hace falta que vuelvas a hacer todo el trabajo. Aquí le damos el procedimiento y las herramientas que necesita para actualizar su estructura (POSAC y QUAPOSQ) y rellenar los valores de estos atributos para las tablas existentes.
Descargue aquí los archivos del procedimiento.
Las distintas etapas de esta actualización son las siguientes:
1- Desde un directorio que contenga todos los archivos .000 que ya se han cargado en la base de datos, ejecute un archivo .bat, upgrade_V2_to_V3.bat, que crea un geopackage update.gpkg que contiene las primitivas (IsolatedNode, ConnectedNode,Edge) de los archivos S57.
.\ upgrade_V2_to_V3.ba directorio_archivos_000 ruta/update.gpkg
2- Ejecuta el script append_tables_upgrade.py en la consola Python de QGis. Asegúrese de cambiar la ruta para el archivo update.gpkg
input_geopackages = [«c:/testgpkgV2/update.gpkg»]
y para su base de datos ENC
output_gpkg = ogr.Open(«c:/enc_db/ENC4.gpkg», 1)
Este script añade las primitivas a la base de datos existente.
3- En la consola Python de QGis, ejecuta el script add_fields_upgrade.py. Asegúrate de cambiar la ruta al archivo .gpkg de tu base de datos
geopackage = ogr.Open(«c:/enc_db/ENC4.gpkg», 1)
Este script crea los atributos POSACC y QUAPOS en todas las tablas de la base de datos.
4-En la consola Python de QGis, ejecute el script add_posacc_quapos_upgrade_points.py. Asegúrese de cambiar la ruta al archivo .gpkg de su base de datos y añada tablas para excluir si ha creado otras tablas en la base de datos.
#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
Este script rellenará los campos POSACC y QUAPOS en todas las tablas de tipo Punto con los valores contenidos en la primitiva IsolatedNode.
5-En la consola Python de QGis, ejecute el script add_posacc_quapos_upgrade_lines.py. Asegúrese de cambiar la ruta al archivo .gpkg de su base de datos y añada tablas a excluir si ha creado otras tablas en la base de datos.
#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
Este script rellenará los campos POSACC y QUAPOS en todas las tablas de tipo Línea con los valores contenidos en la primitiva Edge.
