Tras la publicación de la primera parte del trabajo sobre Geopackages, el trabajo con PostgreSQL/Postgis nos ha permitido realizar una serie de avances adicionales. En este artículo encontrará una actualización sobre el procedimiento de importación de archivos S57, y en un artículo posterior encontrará actualizaciones sobre las simbologías S57 en QGis para Geopackages.
Actualización de la importación de ficheros S57: adición del atributo «finalidad
En el artículo original hablábamos del problema de los duplicados al gestionar la base de datos. Estos duplicados pueden ser de dos tipos diferentes:
- duplicados «verdaderos» en los que todos los atributos son iguales. Estos duplicados son poco frecuentes porque son el resultado de la superposición de dos gráficos de escala muy similar. Sin embargo, también pueden ser el resultado de un error al recargar un archivo S57 duplicado. Este tipo de duplicado se gestiona mediante el script proporcionado en el artículo original, que no se actualiza.
- Duplicados «falsos», en los que la información de la capa se duplica sin que los identificadores del registro sean necesariamente idénticos; se producen cuando se cargan en la misma base de datos zonas cartografiadas a diferentes escalas. Este tipo de duplicación no debe eliminarse, sino gestionarse.
En el script proporcionado en el artículo inicial, añadimos el nombre del archivo y la escala a las tablas del geopackage. Utilizar la escala para gestionar los duplicados «falsos» es bastante complicado:
En el siguiente ejemplo, se muestra el resultado de visualizar los datos brutos del Goulet de la Rade de Brest, que se encuentra en dos mapas:
Para la Rade de Brest, los dos mapas se superponen. La escala de uno es 10000 y la del otro 22500.
Además de la escala, hay un dato que puede ser muy útil: la finalidad del mapa. Se trata de un valor comprendido entre 1 y 6 y corresponde al objetivo principal del mapa:
- 1: Visión general
- 2: General
- 3: Costero
- 4: Aproximación
- 5: Puerto
- 6: Atraque
Los dos mapas superpuestos de nuestro ejemplo tienen dos propósitos diferentes, uno con propósito=5 (Puerto), el otro con propósito=6 (Atraque).
Si filtramos con propósito=5, obtenemos:
Y si filtramos con propósito=6:
Hemos importado 670 archivos S57 en la base de datos del proyecto utilizando PostgreSQL/Postgis.
Los valores mínimo y máximo de la escala de todos los mapas para el objetivo 5 son 3000 y 30000. Los valores mínimo y máximo de escala para el objetivo 6 son 6000 y 15000. Es evidente que los valores de escala de los mapas más detallados se encuentran dentro de los mapas de tipo 5.
Por lo tanto, proponemos que, además del nombre del archivo y la escala, se añada el objetivo como atributo a las tablas de geopaquetes. Este atributo se encuentra en la misma tabla DSID utilizada para recuperar la escala, con el nombre DSID_INTU.
Los archivos .bat de importación siguen siendo los mismos, excepto que el script Python llamado debe ser el siguiente script:
import sys
from osgeo import ogr
def get_default_scale(geopackage_path, default_enc):
# Ouvrir le GeoPackage en mode lecture seule
geopackage = ogr.Open(geopackage_path, 0)
default_enc = default_enc +'.000'
if geopackage is not None:
# Récupérer la couche DSID
dsid_layer = geopackage.GetLayerByName('DSID')
# Vérifier si la couche DSID existe
if dsid_layer:
# Rechercher la valeur de DSID_CSCL correspondant à DSID_DSNM
dsid_layer.SetAttributeFilter(f"DSID_DSNM = '{default_enc}'")
feature = dsid_layer.GetNextFeature()
if feature:
default_scale = feature.GetField('DSPM_CSCL')
return default_scale
def get_default_purpose(geopackage_path, default_enc):
# Ouvrir le GeoPackage en mode lecture seule
geopackage = ogr.Open(geopackage_path, 0)
default_purpose = '0'
default_enc = default_enc +'.000'
if geopackage is not None:
# Récupérer la couche DSID
dsid_layer = geopackage.GetLayerByName('DSID')
# Vérifier si la couche DSID existe
if dsid_layer:
# Rechercher la valeur de DSID_INTU correspondant à DSID_DSNM
dsid_layer.SetAttributeFilter(f"DSID_DSNM = '{default_enc}'")
feature = dsid_layer.GetNextFeature()
if feature:
default_purpose = feature.GetField('DSID_INTU')
return default_purpose
def delete_empty_tables_in_geopackages_dsid(geopackage_path, default_enc):
# Récupérer la valeur de scale à partir de la table DSID
default_scale = get_default_scale(geopackage_path, default_enc)
default_purpose = get_default_purpose(geopackage_path, default_enc)
# Continuer avec le reste du code
if default_scale is None:
print("La valeur par défaut de l'échelle n'a pas été trouvée dans la table DSID.")
return
default_scale_str = str(default_scale)
default_purpose_str = str(default_purpose)
# Ouvrir le GeoPackage en mode édition
geopackage = ogr.Open(geopackage_path, 1)
if geopackage is not None:
# Récupérer le nombre de tables dans le GeoPackage
num_tables = geopackage.GetLayerCount()
# Liste des tables à supprimer
tables_to_delete = []
tables_to_update = []
# Identifier les tables à supprimer
for i in range(num_tables):
table = geopackage.GetLayerByIndex(i)
# Vérifier si la table est vide (aucun enregistrement)
if table.GetFeatureCount() == 0:
tables_to_delete.append(table.GetName())
else:
tables_to_update.append(table.GetName())
# Supprimer les tables
for table_name in tables_to_delete:
geopackage.DeleteLayer(table_name)
print(f"Table supprimée dans {geopackage_path}: {table_name}")
# Mettre à jour les tables restantes
for table_name in tables_to_update:
table = geopackage.GetLayerByName(table_name)
# Vérifier si le champ 'enc_chart' existe
enc_chart_index = table.FindFieldIndex('enc_chart', 1)
if enc_chart_index < 0:
# Ajouter le champ 'enc_chart' s'il n'existe pas
champ1 = ogr.FieldDefn('enc_chart', ogr.OFTString)
champ1.SetWidth(50)
champ1.SetNullable(True)
champ1.SetDefault(default_enc)
table.CreateField(champ1)
# Vérifier si le champ 'scale' existe
scale_index = table.FindFieldIndex('scale', 1)
if scale_index < 0:
# Ajouter le champ 'scale' s'il n'existe pas
champ2 = ogr.FieldDefn('scale', ogr.OFTReal)
champ2.SetWidth(10)
champ2.SetPrecision(2)
champ2.SetDefault(default_scale_str) # Convertir en nombre flottant
table.CreateField(champ2)
# Vérifier si le champ 'purpose' existe
purpose_index = table.FindFieldIndex('purpose', 1)
if purpose_index < 0:
# Ajouter le champ 'purpose' s'il n'existe pas
champ3 = ogr.FieldDefn('purpose', ogr.OFTReal)
champ3.SetWidth(10)
champ3.SetPrecision(2)
champ3.SetDefault(default_purpose_str) # Convertir en nombre flottant
table.CreateField(champ3)
# Mettre à jour les valeurs dans les champs seulement si les champs n'ont pas déjà de valeurs
for feature in table:
enc_chart_value = feature.GetField('enc_chart')
scale_value = feature.GetField('scale')
purpose_value = feature.GetField('purpose')
# Vérifier si les champs n'ont pas déjà de valeurs définies
if enc_chart_value is None or scale_value is None or purpose_value is None:
feature.SetField('enc_chart', default_enc)
feature.SetField('scale', default_scale)
feature.SetField('purpose', default_purpose)
table.SetFeature(feature)
# Fermer le GeoPackage
geopackage = None
else:
print(f"Impossible d'ouvrir le GeoPackage {geopackage_path} en mode édition.")
if __name__ == "__main__":
if len(sys.argv) != 3:
print("Usage: python script.py geopackage_path default_enc")
sys.exit(1)
geopackage_path = sys.argv[1]
default_enc = sys.argv[2]
delete_empty_tables_in_geopackages_dsid(geopackage_path, default_enc)Puedes descargar el script aquí
Copia el script en la misma ubicación que los scripts del artículo original.
Puedes reemplazar las líneas en los archivos .bat
python c:/testgpkg/update_geopackage_dsid.py «%output_geopackage%» «!file!»
en
python c:/testgpkg/update_geopackage_dsid_prp.py «%output_geopackage%» «!file!»
O si no, elimina el script update_geopackage_dsid.py y cambia el nombre de update_geopackage_dsid_prp.py por update_geopackage_dsid.py, lo que te ahorrará modificar los archivos .bat.
Añadir un script de Python para filtrar capas por propósito
Para filtrar la visualización de todas las capas cargadas en el proyecto QGis, puedes utilizar el siguiente script:
# Importer le module QGIS
from qgis.core import QgsProject, QgsMapLayerType
# Définir le valeur de l'attribut "purpose" à filtrer
valeur_purpose = 3
# Obtenir le projet actif
projet = QgsProject.instance()
# Obtenir la liste des couches chargées dans le projet
couches = projet.mapLayers().values()
# Parcourir toutes les couches
for couche in couches:
# Vérifier si la couche est de type vecteur
if couche.type() == QgsMapLayerType.VectorLayer:
# Vérifier si la couche a un champ nommé "purpose"
if couche.fields().indexFromName('purpose') != -1:
# Vérifier si le nom de la couche commence par l'un des préfixes spécifiés
if couche.name().startswith(('pt_', 'li_', 'pl_')):
# Définir le filtre sur l'attribut "purpose" égal à la valeur spécifiée
filtre = f'"purpose" = {valeur_purpose}'
couche.setSubsetString(filtre)
print(f"Filtre défini pour la couche {couche.name()}")
else:
print(f"La couche {couche.name()} ne commence pas par 'pt_', 'li_', ou 'pl_'")
else:
print(f"La couche {couche.name()} n'a pas d'attribut 'purpose'")
else:
print(f"La couche {couche.name()} n'est pas une couche vecteur")
Puede descargar el script aquí.
Para eliminar los filtros Purpose, puede utilizar el mismo script que para el escalado (unfilter.py).
El proyecto de base de datos S57 con Geopackage ha llegado a su fin. Ahora estamos trabajando en la creación de un procedimiento equivalente utilizando una base de datos PostgreSQL/Postgis. Ayúdenos a completar este proyecto
