Após a publicação da primeira parte do trabalho sobre Geopacotes, o trabalho com o PostgreSQL/Postgis nos permitiu fazer uma série de avanços adicionais. Neste artigo, você encontrará uma atualização sobre o procedimento de importação de arquivos S57 e, em um artigo posterior, encontrará atualizações sobre as simbologias S57 no QGis para Geopackages.
Atualização da importação de arquivos S57: adição do atributo “purpose” (finalidade)
No artigo original, discutimos o problema das duplicatas ao gerenciar o banco de dados. Essas duplicatas podem ser de dois tipos diferentes:
- Duplicatas “verdadeiras”, em que todos os atributos são os mesmos. Elas são raras porque resultam da sobreposição de duas sobreposições de gráficos de escala muito semelhante. No entanto, elas também podem ser o resultado de um erro ao recarregar um arquivo S57 duplicado. Esse tipo de duplicata é tratado usando o script fornecido no artigo original, que não é atualizado.
- Duplicatas “falsas”, em que as informações da camada são duplicadas sem que os identificadores de registro sejam necessariamente idênticos; elas ocorrem quando áreas mapeadas em escalas diferentes são carregadas no mesmo banco de dados. Esse tipo de duplicação não deve ser excluído, mas gerenciado.
No script fornecido no artigo inicial, adicionamos o nome do arquivo e a escala às tabelas do geopackage. Usar a escala para gerenciar duplicatas “falsas” é bastante complicado:
No exemplo a seguir, aqui está o resultado da exibição dos dados brutos do Goulet de la Rade de Brest, que está em dois mapas:
Para o Rade de Brest, os dois mapas se sobrepõem. A escala de um é 10000 e a do outro é 22500.
Além da escala, há uma informação que pode ser muito útil: a finalidade do mapa. Esse é um valor entre 1 e 6 e corresponde ao objetivo principal do mapa:
- 1: Visão geral
- 2: Geral
- 3: Litoral
- 4: Aproximação
- 5: Porto
- 6: Atracação
Os dois mapas sobrepostos em nosso exemplo têm duas finalidades diferentes: um com finalidade=5 (Porto) e o outro com finalidade=6 (Atracação).
Se filtrarmos com finalidade=5, obteremos:
E se filtrarmos com propósito=6:
Importamos 670 arquivos S57 para o banco de dados do projeto usando o PostgreSQL/Postgis.
Os valores mínimo e máximo da escala de todos os mapas para o objetivo 5 são 3000 e 30000. Os valores mínimo e máximo de escala para o objetivo 6 são 6000 e 15000. É evidente que os valores de escala dos mapas mais detalhados encontram-se nos mapas do tipo 5.
Portanto, propomos que, além do nome do arquivo e da escala, o objetivo seja adicionado como um atributo às tabelas do geopacote. Esse atributo é encontrado na mesma tabela DSID usada para recuperar a escala, com o nome DSID_INTU.
Os arquivos .bat de importação permanecem os mesmos, exceto pelo fato de que o script Python chamado deve ser o seguinte:
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)Você pode fazer o download do script aqui
Copie o script para o mesmo local que os scripts do artigo original.
Você pode substituir as linhas nos arquivos .bat
python c:/testgpkg/update_geopackage_dsid.py “%output_geopackage%” “!file!”
por
python c:/testgpkg/update_geopackage_dsid_prp.py “%output_geopackage%” “!file!”
Como alternativa, exclua o script update_geopackage_dsid.py e renomeie update_geopackage_dsid_prp.py para update_geopackage_dsid.py, o que evita que você modifique os arquivos .bat.
Adicionar um script Python para filtrar camadas por finalidade
Para filtrar a exibição de todas as camadas carregadas no projeto QGis, você pode usar o seguinte 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")
Você pode baixar o script aqui.
Para eliminar os filtros Purpose, você pode usar o mesmo script usado para dimensionamento (unfilter.py).
O projeto do banco de dados S57 com o Geopackage chegou ao fim. Agora estamos trabalhando na configuração de um procedimento equivalente usando um banco de dados PostgreSQL/Postgis. Ajude-nos a concluir esse projeto
