La primera parte del proyecto Financiación Colaborativa para la Integración de Datos Marinos en QGIS ha llegado a su fin, gracias a las contribuciones de ORANGE Marine, Geoconceptos Uruguay, Janez Avzec y Gaetan Mas. Nuestro más sincero agradecimiento a todos ellos. Por ello, publicamos el contenido final de esta parte del trabajo.
Dado el volumen del resultado, lo publicamos en dos partes: la primera se refiere a la creación y gestión de la base de datos, la segunda a la simbología S57 en QGis.
El formato S57 es complejo y no se ajusta plenamente a las normas de formato SIG. Por lo tanto, su conversión con GDAL es compleja o incompleta. En este artículo entraremos en detalle en la transformación completa de S57 en tablas de archivos Geopackage.
Tipos de geometría
Las «capas» del S57 son clases de objetos. Por ejemplo, las diferentes zonas terrestres se codifican en la clase de objeto LNDARE.
La definición de esta clase de objeto es
Geo object: Land area (LNDARE) (P,L,A)
Attributes: CONDTN OBJNAM NOBJNM STATUS INFORM NINFOM
Aunque todos los objetos LNDARE tienen los mismos atributos, no ocurre lo mismo con el tipo de geometría. La información (P,L,A) indica que en esta clase de objeto pueden encontrarse puntos, líneas y polígonos. Contrariamente a las normas SIG, los tres tipos de geometría coexisten dentro de la misma clase de objeto.
Sea cual sea el formato que elijamos para integrar en QGis, tendremos que crear una capa para cada tipo de geometría. Si no lo hacemos, GDAL creará el tipo de capa basándose en la primera entidad encontrada durante la conversión. Si es de tipo punto, la capa creada será de tipo punto y las entidades líneas y polígonos serán ignoradas. Del mismo modo, si la primera entidad es de tipo línea, se ignorarán los puntos y polígonos.
También hay que tener en cuenta que el formato S57 no tiene restricciones sobre las clases de objetos: si no hay ninguna clase de objeto en el área cubierta por el fichero S57, no habrá nada sobre él en el fichero (ninguna capa vacía). Del mismo modo, si sólo hay un tipo de geometría presente, aunque los tres tipos sean posibles, no habrá rastro de los otros tipos de geometría.
Por lo tanto, el tratamiento de un fichero S57 con ogr2ogr debe dividirse en tres etapas, una para cada tipo de geometría. Las siguientes opciones permiten procesar cada clase de objeto S57 seleccionando un solo tipo de geometría:
-where « OGR_GEOMETRY=’POINT’ or OGR_GEOMETRY=’MULTIPOINT’ »
-where « OGR_GEOMETRY=’LINESTRING’ or OGR_GEOMETRY=’MULTILINESTRING’ »
-where « OGR_GEOMETRY=’POLYGON’ or OGR_GEOMETRY=’MULTIPOLYGON’ »
Para determinados tipos de controlador, GDAL permite crear prefijos en las tablas de salida. En este caso, podrías crear todas las tablas (puntos, líneas, polígonos) en un único geopackage, prefijándolas con pt_, li_ y pl_, por ejemplo. El problema es que el controlador S57 de GDAL no permite esta opción. Por lo tanto, tenemos que crear tres geopaquetes distintos, en los que se crearán las tablas según el tipo de geometría. Cada geopackage contendrá entonces una tabla con el mismo nombre pero con una geometría diferente. Aunque ésta es la solución más sencilla, no es la más elegante ni la más práctica. Aquí vamos a utilizar tres geopackages de importación para los comandos ogr2ogr en los que importaremos tablas desde un fichero S57. Los llamaremos PointsENC, LinesENC y PolysENC. También crearemos un único geopackage llamado ENC para toda la base de datos. Para cada mapa ENC importaremos su contenido en los tres geopackages de importación usando ogr2ogr, luego ejecutaremos scripts Python para actualizar la base de datos del geopackage ENC con las nuevas tablas importadas. Tenemos que utilizar scripts de Python en lugar de consultas y funciones SQL porque la versión SQL de SQLite, utilizada para los geopaquetes, es muy limitada y no permite la creación de procedimientos o funciones.
Si tienes un único archivo S57 para importar, el proceso es el siguiente:
- Creación del geopackage ENC vacío.
- Carga de las clases de objetos Point en el geopackage PointsENC mediante ogr2ogr
- Carga de las clases de objetos Lines en el geopackage LinesENC utilizando ogr2ogr
- Cargar las clases de objeto Polygon en el geopackage PolysENC utilizando ogr2ogr
- Eliminación de las tablas vacías de los tres geopackages y adición de los atributos map_enc y scale…
- Clonación de las tablas de los tres geopackages de importación en el geopackage ENC, prefijando las tablas pointsENC con pt_, las tablas LinesENC con li_ y las tablas PolysENC con pl_.
Si tienes un lote de archivos S57 para cargar simultáneamente, el proceso es el siguiente:
- Cargar las clases de objeto Point en el geopackage pointsENC con ogr2ogr, borrando las tablas vacías y añadiendo el nombre del fichero enc y la escala en los atributos de todas las tablas.
- Cargar las clases de objeto Líneas en el geopackage LinesENC con ogr2ogr, borrando las tablas vacías y añadiendo el nombre del fichero enc y la escala en los atributos de todas las tablas.
- Cargar las clases de objeto Polygon en el geopackage PolysENC con ogr2ogr, borrando las tablas vacías y añadiendo el nombre del archivo enc y la escala en los atributos de todas las tablas.
- Actualización de las tablas existentes en el geopackage ENC a partir de los tres geopackages de importación ET.
- Clonación de tablas de los tres geopacks de importación que no estaban presentes en el geopackage ENC,
- Detección y supresión de cualquier duplicado resultante de la actualización de tablas en el geopaquete ENC (opcional).
Tratamiento de sondas batimétricas
Las sondas batimétricas plantean varios problemas a la hora de convertir formatos. El primero es que el valor de profundidad no está contenido en un atributo, sino como Z en la geometría (XYZ). El segundo es que las sondas no son de tipo Punto sino Multipunto. Para obtener el valor de las sondas directamente en un campo de atributo, es necesario añadir dos parámetros a la línea de comandos de ogr2ogr:
-oo SPLIT_MULTIPOINT=ON -oo ADD_SOUNDG_DEPTH=ON
Recuperación de identificadores «padre
Algunos objetos pueden agruparse por un objeto padre no geográfico. Por ejemplo, los sectores de incendio se codifican con un sector por registro en la tabla «LUCES». Los distintos sectores de un mismo incendio no tienen ningún atributo particular que pueda indicar que corresponden a la misma entidad. Esta información figura en un registro «padre». Para recuperar este identificador como atributo de la tabla «LIGTHS», es necesario añadir una opción a la línea de comandos:
-oo RETURN_LINKAGES=O
Esta opción crea un atributo name_rcid común a todos los sectores de un mismo semáforo, pero también crea una serie de campos (name_rcnm,ORNT,USAG,MASK).
Tratamiento de tipo «lista»
Además de los tipos de campo tradicionales (entero, real, texto), el formato S57 utiliza un tipo especial: las listas de cadenas de caracteres o enteros. Estos tipos de campo se encuentran en PostgreSQL, pero no en shapefiles y geopackage.
En ausencia de opciones específicas, aparecerán mensajes de advertencia indicando que el formato no está soportado y el contenido de estos campos será ignorado. Para evitar perder el contenido, es necesario añadir la extensión :
-mapFieldType StringList=String,IntegerList=String
Esta opción toma una lista {..,…} y genera una cadena de texto con un primer número que indica el número de elementos de la lista, dos puntos y, a continuación, los elementos de la lista separados por comas:
{3} -> ‘(1:3)’
{3,1} -> ‘(2:3,1)’
Esto complica inevitablemente un poco el tratamiento de estos campos, pero la información no se pierde. Por otro lado, el controlador GDAL presenta un problema para gestionar el formato String resultante.
En primer lugar, no hay ningún problema con los atributos de los objetos S57. El campo de la tabla geopackage es de tipo ‘Texto’, lo que permite crear una cadena de caracteres de cualquier longitud.
Sin embargo, para los «identificadores padre», la definición utilizada por el controlador es «Text()», es decir, una cadena de caracteres con una longitud máxima. Esta longitud viene determinada a priori por el contenido del primer registro procesado. En otras palabras, sólo por casualidad tendrá un campo con la longitud necesaria para incluir todos los valores del fichero.
Verá entonces mensajes del tipo
WARNING : Value of field ‘NAME_RCID’ has 371 characters, whereas maximum allowed is 30
La solución no es complicada, pero implica bastante trabajo. Partimos de la base de que para la visualización clásica de las tarjetas ENC, esto no representa ningún problema y puedes simplemente ignorar estos mensajes de advertencia.
Sin embargo, si necesitas un procesamiento más avanzado que incluya esta información, esta es la única manera que he encontrado hasta ahora:
- Ejecuta el comando ogr2ogr para crear el geopackage de importación
- Localiza el valor máximo indicado para cada tipo de campo
- Para la tabla o tablas de salida que se han creado y de las que se quiere recuperar toda la información, abre el gestor de base de datos QGis
- Elimina todos los registros de la tabla (DELETE FROM nombre_tabla)
- Vé al menú Tabla->Modificar una tabla…
- Selecciona el campo que desees de la lista de campos, luego haz clic en Editar una columna
- Introduce el valor deseado en el campo Longitud, luego haz clic en OK
- Haz lo mismo para todos los demás campos que te interesen
- Elimina todas las demás tablas del geopackage de importación que no planteen problemas
- Vuelve a ejecutar el comando de importación ogr2ogr.
Flujo de trabajo para cargar un único archivo S57
Creación del geopackage ENC
Vamos a crear un geopackage vacío que recibirá las tablas importadas, y luego las sucesivas actualizaciones cuando se carguen otros archivos S57.
En el panel Explorador de QGis, haz clic con el botón derecho del ratón en Geopackages y selecciona Crear una base de datos
En la ventana que se abre, introduce el nombre del geopackage, deja la opción por defecto para el nombre de la tabla y selecciona Sin geometría en Tipo de geometría.
Haz clic en Aceptar para que el nuevo geopackage vacío se añada a las conexiones del geopackage de QGis. Es aconsejable borrar la tabla ENC vacía creada al crear el geopackage. Para ello, abre el gestor de bases de datos de QGis desde el menú principal, sitúate sobre la tabla ENC del geopackage ENC.gpkg y haz clic con el botón derecho del ratón para seleccionar Borrar…
comandos ogr2ogr para crear geopaquetes de importación
Con todos estos elementos en mente, he aquí las tres líneas de comandos ogr2ogr para crear las tablas de los geopaquetes de importación:
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
Para ejecutar estas líneas de comando, simplemente abre la ventana shell de OSGeo4W:
Se abre la ventana Shell
Introduce las líneas de comando, teniendo cuidado de introducir la ruta y el nombre del archivo con la extensión .000.
El resultado será una serie de tablas creadas en cada geopackage de importación. Sin embargo, algunas tablas estarán completamente vacías. Si el tipo de geometría es posible para una clase de objeto, la tabla se creará, pero si no hay ocurrencias en el archivo S57 procesado, la tabla no tendrá registros.
Uso de scripts Python para gestionar la base de datos
Utilizaremos la consola Python de QGis para todas las operaciones de gestión de la base de datos.
Para abrirla, vaya al menú Extensiones -> Consola Python. Se abrirán los paneles de la consola Python:
Para ejecutar un script Python, carga el código en el panel Editor. Puedes copiar y pegar o utilizar el icono Directorio para apuntar al archivo .py. La sangría es crucial para los scripts Python, así que elige la segunda opción.
Elimine las tablas vacías y añada el nombre y la escala del mapa
Como hemos visto antes, al ejecutar los comandos ogr2ogr, se crean tablas vacías para las capas en las que no se solicita geometría en la línea de comandos. Por tanto, eliminaremos estas tablas vacías y aprovecharemos el script de Python para añadir dos atributos a todas las tablas: el nombre del mapa ENC y su escala.
En cuanto al nombre, puede ser útil en el futuro mantenimiento de la base de datos para poder seleccionar las entidades correspondientes a un fichero fuente S57. Esta información no figura en ningún atributo del fichero S57. Por lo tanto, hay que introducirla manualmente.
Lo mismo ocurre con la escala del mapa. Aunque en el S57 se proporciona un atributo para introducir la escala del mapa «de papel», el atributo CSCALE de la tabla M_CSCL, es muy raro que esta tabla y su atributo estén presentes en los ficheros S57.
Si bien el nombre del mapa tiene una utilidad relativa, no puede decirse lo mismo de la escala, ya que se trata de un criterio esencial en la búsqueda de duplicados. En efecto, cuando se mezclan entidades procedentes de varios ficheros S57, los datos de diferentes escalas coexistirán más o menos felizmente.
El siguiente código Python suprime las tablas vacías de un esquema de importación y añade dos atributos (enc_chart y scale) a todas las tablas:
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)
Ten cuidado de respetar la sangría del código Python. Para simplificar la tarea, descarga el archivo .py directamente desde este enlace: delete_empty_tables_in_geopackages.py.
Para ejecutarlo, asegúrate de cambiar los parámetros de entrada del script:
Sólo tendrás que introducir las geopackage_paths la primera vez que ejecutes el script. Por otro lado, el nombre del mapa ENC y la escala tendrán que cambiarse cuando se importe cada nuevo archivo S57.
Al final de esta ejecución tendrás los nuevos atributos en cada tabla
Flujo de trabajo para cargar varios archivos S57 al mismo tiempo
Las líneas de comando de los párrafos anteriores suponen que tienes un único archivo S57 para cargar. Pero puedes cargar varios archivos simultáneamente. Utilizando las opciones -append -update, puedes procesar todos los archivos .000 de un directorio en un bucle. Entonces tendrás todos los archivos cargados en los geopacks de importación.
A diferencia de la carga de un único archivo, cuando cargamos un lote de archivos necesitamos incluir la adición del nombre del archivo y la escala en el bucle de transcodificación S57->geopackage. Para ello, utilizaremos un script python que añade los atributos enc_chart y scale a las tablas creadas por el comando ogr2ogr y las rellena con el nombre del archivo correspondiente y una escala definida por el productor. El script también elimina cualquier tabla vacía en el geopackage de importación.
Gracias a Janez AVSEC por la idea de utilizar la tabla DSID para recuperar la escala de compilación de datos.
Los archivos .bat deben ser ejecutados desde una ventana de OSGeo4W para poder beneficiarse del entorno Python de QGis.
La llamada a cada archivo .bat es de la forma
.\fichier.bat «directorio que contiene los archivos S57 a cargar» «nombre del geopackage de importación»
Por ejemplo:
.\mass_load_s57_polys.bat «c:/enc_db/enc000» «c:/enc_db/polys.gpkg»
Los archivos .bat para esta tarea son los siguientes:
Para importar tablas de tipo de punto:
@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
Observará que, además del comando ogr2ogr para importar capas de puntos, importamos dos tablas sin geometría: DSID y C_AGGR. La tabla DSID contiene dos atributos con el nombre del archivo S57 y la escala de compilación de datos. Esta tabla será utilizada por el script de Python (update_geopackage_dsid.py) para introducir el archivo S57 original y la escala de compilación en cada tabla de geopackage.
Para importar tablas de tipo línea:
@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é.
pausePara importar tablas de polígonos:
@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
Los archivos .bat contienen el comando ogr2ogr y una línea para ejecutar un script de Python, update_geopackage_dsid.py. Asegúrate de cambiar la ruta de este script en el archivo .bat por la de tu elección. Este script añade los dos atributos enc_chart y scale a todas las tablas del geopackage.
El código python utilizado es el siguiente:
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)Podrías hacer un único archivo .bat para los tres procesos, pero necesitas comprobar los resultados para cada tipo de geometría.
Puedes descargar los archivos .bat y el código Python haciendo clic aquí.
Continúa independientemente del número de archivos cargados
Clonar y/o añadir tablas importadas al geopackage ENC
Para la primera carga del S57, el siguiente script de Python crea todas las tablas presentes en los tres geopacks de importación. Para las cargas siguientes, 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.
La clonación se realiza mediante el siguiente código:
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()
Puedes descargar el archivo .py aquí: clone_or_append_tables_with_prefix.py.
Gestión de duplicados
Cuando se superponen dos mapas ENC, las entidades presentes en la zona de superposición se insertarán tantas veces como superposiciones haya. Para evitar duplicados y más, la base de datos debe limpiarse después de cada importación de un fichero S57.
Pero hay dos tipos diferentes de duplicación. Por un lado, las zonas en las que se solapan dos mapas de escala similar producirán duplicados «reales» que habrá que limpiar. Por otro lado, los mapas superpuestos de escalas muy diferentes deben tratarse de forma diferente. El mapa menos detallado (escala pequeña) tendrá una selección de datos del mapa más detallado (escala grande). En este segundo caso, no es aconsejable eliminar esta información.
El código Python siguiente eliminará los duplicados «reales». Para otros duplicados, recomendamos utilizar filtros QGis. Dependiendo de sus necesidades, puedes utilizar el atributo «escala» que añadimos al eliminar las tablas vacías para definir qué datos se muestran en su mapa y cuáles no.
En este ejemplo, sólo se mostrarán en el mapa los elementos correspondientes a una escala comprendida entre 12500 y 45000.
Esta operación corresponde a filtrar una capa, pero… ¿qué haces cuando tienes 130 capas desplegadas? Para ahorrarte tener que realizar esta operación 130 veces, aquí tienes un código Python que te permite filtrar todas las capas del proyecto:
#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'")
Y otra para eliminar el filtro de todas las capas mostradas:
# 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()}")
Descarga los códigos python aquí.
El siguiente código crea una tabla temporal (temp_table) para cada tabla del geopackage ENC. A continuación, rellena esta tabla con los registros separados basados en LNAM. A continuación, borra la tabla original y renombra la tabla temporal con el nombre de la tabla original.
Este tratamiento se aplica a todas las tablas S57 excepto:
- las tablas de luces (LIGHTS) y de sondeos batimétricos (SOUNDG), en las que el atributo LNAM no es útil. En efecto, cada sector de incendio tiene el mismo LNAM y name_rcid pero corresponde a un color y a un sector geográfico diferentes. Del mismo modo, la LNAM y el name_rcid de los sondeos batimétricos son los mismos para todos los puntos de sondeo de la entidad multipunto original. En el primer caso tenemos que buscar colores y sectores duplicados para el mismo name_rcid, y para los sondeos tenemos que encontrar los sondeos con el mismo XY en su geometría.
- Tablas SEAARE, LNDRGN y ADMARE que cubren grandes áreas y que se dividen según la huella del mapa, manteniendo el mismo LNAM. En este caso, aunque la LNAM sea la misma, los polígonos no tienen la misma forma y área.
- las tablas DSID y C_AGGR, que no tienen LNAM y que, por principio, no pueden tener duplicados.
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()
Puedes descargar el archivo .py desde este enlace: supprime_doublons.py. No olvides cambiar las rutas y los nombres de archivo para que coincidan con los tuyos. Fíjate también en la línea :
if table_name not in [«pl_SEAARE», «pl_LNDRGN», «pl_ADMARE», «layer_styles», «natsurf»,»DSID»,»C_AGGR»]:
En principio, sólo tendrás las tablas prefijadas y las tablas layer_styles y natsurf. Esta línea evita que se procesen pl_SEAARE, pl_LNDRGN y pl_ADMARE, pero sobre todo evita que se vacíen las tablas layer_styles y natsurf que utilizamos para las simbologías. Si necesitas crear tablas distintas a las del fichero s57 del geopackage, debes tener en cuenta los nombres de estas tablas en este script, añadiéndolos a esta línea de excepciones.
El proyecto de base de datos S57 con Geopackage ha llegado a su fin. Ahora estamos empezando a implementar un procedimiento equivalente utilizando una base de datos PostgreSQL/Postgis. Por favor, ¡ayúdanos a completar este proyecto!
