A primeira parte do projeto Financiamento colaborativo para a integração de dados marinhos no QGIS chegou ao fim, graças às contribuições de ORANGE Marine, Geoconceptos Uruguay, Janez Avzec e Gaetan Mas. Nossos mais sinceros agradecimentos a eles. Portanto, estamos publicando o conteúdo final dessa parte do trabalho.
Devido ao volume do resultado, estamos publicando-o em duas partes: a primeira diz respeito à criação e ao gerenciamento do banco de dados, a segunda à simbologia S57 no QGis.
O formato S57 é complexo e não está totalmente em conformidade com os padrões de formato GIS. Sua conversão com o GDAL é, portanto, complexa ou incompleta. Neste artigo, entraremos em detalhes sobre a transformação completa do S57 em tabelas de arquivos do Geopackage.
Tipos de geometria
As “camadas” do S57 são classes de objetos. Por exemplo, as diferentes áreas de terra são codificadas na classe de objeto LNDARE.
A definição dessa classe de objeto é:
Geo object: Land area (LNDARE) (P,L,A)
Attributes: CONDTN OBJNAM NOBJNM STATUS INFORM NINFOM
Embora todos os objetos LNDARE tenham os mesmos atributos, o mesmo não pode ser dito sobre o tipo de geometria. A informação (P,L,A) indica que pontos, linhas e polígonos podem ser encontrados nessa classe de objeto. Ao contrário dos padrões do GIS, os três tipos de geometria coexistem na mesma classe de objeto.
Seja qual for o formato que escolhermos para integrar ao QGis, precisaremos criar uma camada para cada tipo de geometria. Se não fizermos isso, o GDAL criará o tipo de camada com base na primeira entidade encontrada durante a conversão. Se ela for do tipo ponto, a camada criada será do tipo ponto e as entidades linhas e polígonos serão ignoradas. Da mesma forma, se a primeira entidade for do tipo linha, os pontos e polígonos serão ignorados.
Também deve ser observado que o formato S57 não tem restrições sobre classes de objetos: se não houver nenhuma classe de objeto na área coberta pelo arquivo S57, não haverá nada sobre ela no arquivo (nenhuma camada vazia). Da mesma forma, se houver apenas um tipo de geometria presente, mesmo que todos os três tipos sejam possíveis, não haverá nenhum vestígio dos outros tipos de geometria.
Portanto, o processamento de um arquivo S57 com ogr2ogr deve ser dividido em três estágios, um para cada tipo de geometria. As opções a seguir permitem que você processe cada classe de objeto S57 selecionando apenas um tipo de geometria:
-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 driver, o GDAL permite que você crie prefixos nas tabelas de saída. Nesse caso, você poderia criar todas as tabelas (pontos, linhas, polígonos) em um único geopackage, prefixando-as com pt_, li_ e pl_, por exemplo. O problema é que o driver S57 do GDAL não permite essa opção. Portanto, precisamos criar três geopacotes separados, nos quais as tabelas serão criadas de acordo com o tipo de geometria. Cada geopacote conterá uma tabela com o mesmo nome, mas com uma geometria diferente. Embora essa seja a solução mais simples, não é a mais elegante nem a mais prática. Aqui, usaremos três geopacotes de importação para os comandos ogr2ogr nos quais importaremos tabelas de um arquivo S57. Nós os chamaremos de PointsENC, LinesENC e PolysENC. Também criaremos um único geopackage chamado ENC para todo o banco de dados. Para cada mapa ENC, importaremos seu conteúdo para os três geopacotes de importação usando ogr2ogr e, em seguida, executaremos scripts Python para atualizar o banco de dados do geopacote ENC com as novas tabelas importadas. Temos de usar scripts Python em vez de consultas e funções SQL porque a versão SQL do SQLite, usada para geopacotes, é muito limitada e não permite a criação de procedimentos ou funções.
Se você tiver um único arquivo S57 para importar, o processo é o seguinte:
- Criação do geopacote ENC vazio.
- Carregar classes de objetos de pontos no geopacote pointsENC usando ogr2ogr
- Carregar as classes de objeto Lines no geopackage LinesENC usando ogr2ogr
- Carregar as classes de objeto Polygon no geopackage PolysENC usando ogr2ogr
- Remoção de tabelas vazias dos três geopacotes e adição de atributos map_enc e scale…
- Clonagem das tabelas dos três geopacotes de importação para o geopacote ENC, prefixando as tabelas pointsENC com pt_, as tabelas LinesENC com li_ e as tabelas PolysENC com pl_.
Se você tiver um lote de arquivos S57 para carregar simultaneamente, o processo é o seguinte:
- Carregue as classes de objeto Point no geopackage pointsENC com ogr2ogr, excluindo tabelas vazias e adicionando o nome do arquivo enc e a escala aos atributos de todas as tabelas.
- Carregue as classes de objeto Lines no geopackage LinesENC com ogr2ogr, excluindo tabelas vazias e adicionando o nome do arquivo enc e a escala aos atributos de todas as tabelas.
- Carregar classes de objeto Polygon no geopackage PolysENC usando ogr2ogr, excluindo tabelas vazias e adicionando o nome do arquivo enc e a escala aos atributos de todas as tabelas.
- Atualização de tabelas existentes no geopacote ENC a partir dos três geopacotes de importação ET.
- Clonagem de tabelas dos três geopacotes de importação que não estavam presentes no geopacote ENC,
- Detecção e exclusão de quaisquer duplicatas resultantes da atualização de tabelas no geopacote ENC (opcional)
Processamento de sondagens batimétricas
As sondas batimétricas apresentam vários problemas na conversão de formatos. O primeiro é que o valor da profundidade não está contido em um atributo, mas como Z na geometria (XYZ). O segundo é que as sondas não são do tipo Ponto, mas do tipo Multiponto. Para obter o valor das sondas diretamente em um campo de atributo, é necessário adicionar dois parâmetros à linha de comando do ogr2ogr:
-oo SPLIT_MULTIPOINT=ON -oo ADD_SOUNDG_DEPTH=ON
Recuperação de identificadores “pai”
Alguns objetos podem ser agrupados por um objeto pai não geográfico. Por exemplo, os setores de incêndio são codificados com um setor por registro na tabela “LIGHTS”. Os diferentes setores do mesmo incêndio não têm nenhum atributo específico que possa indicar que eles correspondem à mesma entidade. Essa informação está contida em um registro “pai”. Para recuperar esse identificador como um atributo da tabela LIGTHS, é necessário adicionar uma opção à linha de comando:
-oo RETURN_LINKAGES=O
Essa opção cria um atributo name_rcid que é comum a todos os setores no mesmo semáforo, mas também cria uma série de campos (name_rcnm,ORNT,USAG,MASK).
Processamento do tipo “Lista”
Além dos tipos de campo tradicionais (inteiro, real, texto), o formato S57 usa um tipo especial: listas de cadeias de caracteres ou inteiros. Esses tipos de campo são encontrados no PostgreSQL, mas não nos shapefiles e no geopackage.
Na ausência de opções específicas, haverá mensagens de aviso indicando que o formato não é compatível e o conteúdo desses campos será ignorado. Para evitar a perda do conteúdo, você precisa adicionar a extensão :
-mapFieldType StringList=String,IntegerList=String
Essa opção usa uma lista {…,…} e gera uma string de texto com um primeiro número indicando o número de elementos na lista, dois pontos e, em seguida, os elementos da lista separados por vírgulas:
{3} -> ‘(1:3)’
{3,1} -> ‘(2:3,1)’
Isso inevitavelmente complica um pouco o processamento desses campos, mas as informações não são perdidas. Por outro lado, o driver GDAL apresenta um problema para gerenciar o formato String resultante.
Em primeiro lugar, não há problema com os atributos dos objetos S57. O campo na tabela de geopacotes é do tipo “Text”, que permite a criação de uma cadeia de caracteres de qualquer comprimento.
No entanto, para “parent identifiers”, a definição usada pelo driver é “Text()”, ou seja, uma cadeia de caracteres com um comprimento máximo. Esse comprimento é determinado a priori pelo conteúdo do primeiro registro processado. Em outras palavras, é apenas por acaso que você tem um campo com o comprimento necessário para incluir todos os valores no arquivo.
Em seguida, você verá mensagens do tipo:
WARNING : Value of field ‘NAME_RCID’ has 371 characters, whereas maximum allowed is 30
A solução não é complicada, mas envolve uma boa quantidade de trabalho. Comece com a premissa de que, para a exibição clássica de cartões ENC, isso não representa um problema e você pode simplesmente ignorar essas mensagens de aviso.
Entretanto, se você precisar de um processamento mais avançado que inclua essas informações, essa é a única maneira que encontrei até o momento:
- Execute o comando ogr2ogr para criar o geopacote de importação
- Localize o valor máximo indicado para cada tipo de campo
- Para a tabela ou tabelas de saída que foram criadas e das quais você deseja recuperar todas as informações, abra o gerenciador de banco de dados QGis
- Exclua todos os registros da tabela (DELETE FROM table_name)
- Vá para o menu Tabela->Modificar uma tabela…
- Selecione o campo desejado na lista de campos e, em seguida, clique em Edit a column (Editar uma coluna)
- Digite o valor desejado no campo Length (Comprimento) e clique em OK
- Faça o mesmo com todos os outros campos de seu interesse
- Exclua todas as outras tabelas do geopacote de importação que não representem um problema
- Execute o comando de importação ogr2ogr novamente.
Fluxo de trabalho para fazer upload de um único arquivo S57
Criação do geopacote ENC
Vamos criar um geopackage vazio que receberá as tabelas importadas e, em seguida, atualizações sucessivas quando outros arquivos S57 forem carregados.
No painel QGis Explorer, clique com o botão direito do mouse em Geopackages e selecione Create a database (Criar um banco de dados)
Na janela que se abre, digite o nome do geopacote, deixe a opção padrão para o nome da tabela e selecione No geometry (Sem geometria) em Geometry type (Tipo de geometria).
Haga 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, abra el gestor de bases de datos de QGis desde el menú principal, sitúese sobre la tabla ENC del geopackage ENC.gpkg y haga 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 executar essas linhas de comando, basta abrir a janela do shell do OSGeo4W:
A janela Shell é aberta
Digite as linhas de comando, tomando cuidado para inserir o caminho e o nome do arquivo com a extensão .000.
O resultado será uma série de tabelas criadas em cada geopackage de importação. Entretanto, algumas tabelas ficarão completamente vazias. Se o tipo de geometria for possível para uma classe de objeto, a tabela será criada, mas se não houver ocorrências no arquivo S57 processado, a tabela não terá registros.
Uso de scripts Python para gerenciar o banco de dados
Usaremos o console QGis Python para todas as operações de gerenciamento do banco de dados.
Para abri-lo, acesse o menu Extensões -> Console Python. Os painéis do console Python serão abertos:
Para executar um script Python, carregue o código no painel Editor. Você pode copiar e colar ou usar o ícone Diretório para apontar para o arquivo .py. A indentação é crucial para scripts Python, portanto, escolha a segunda opção.
Remova as tabelas vazias e adicione o nome e a escala do mapa
Como vimos anteriormente, ao executar comandos ogr2ogr, são criadas tabelas vazias para camadas em que não há geometria solicitada na linha de comando. Portanto, removeremos essas tabelas vazias e aproveitaremos o script Python para adicionar dois atributos a todas as tabelas: o nome do mapa ENC e sua escala.
No que diz respeito ao nome, pode ser útil na manutenção futura do banco de dados poder selecionar as entidades correspondentes a um arquivo de origem S57. Essa informação não está contida em nenhum atributo do arquivo S57. Portanto, temos de inseri-la manualmente.
O mesmo se aplica à escala do mapa. Embora um atributo seja fornecido no S57 para inserir a escala do mapa “em papel”, o atributo CSCALE na tabela M_CSCL, é muito raro que essa tabela e seu atributo estejam presentes nos arquivos S57.
Embora o nome do mapa seja de uso relativo, o mesmo não pode ser dito da escala, pois esse é um critério essencial na busca por duplicatas. De fato, quando entidades de vários arquivos S57 são misturadas, os dados de diferentes escalas coexistirão de forma mais ou menos feliz.
O código Python a seguir exclui tabelas vazias de um esquema de importação e adiciona dois atributos (enc_chart e scale) a todas as tabelas:
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)
Tenha o cuidado de respeitar a indentação do código Python. Para simplificar a tarefa, baixe o arquivo .py diretamente deste link: delete_empty_tables_in_geopackages.py.
Para executá-lo, certifique-se de alterar os parâmetros de entrada do script:
Você só precisará inserir os geopackage_paths na primeira vez que executar o script. Por outro lado, o nome do mapa ENC e a escala precisarão ser alterados quando cada novo arquivo S57 for importado.
No final dessa execução, você terá os novos atributos em cada tabela
Fluxo de trabalho para carregar vários arquivos S57 ao mesmo tempo
As linhas de comando nos parágrafos anteriores pressupõem que você tenha um único arquivo S57 para carregar. Mas você pode carregar vários arquivos simultaneamente. Ao usar as opções -append -update, você pode processar todos os arquivos .000 em um diretório em um loop. Em seguida, você terá todos os arquivos carregados nos geopacotes de importação.
Diferentemente do carregamento de um único arquivo, ao carregar um lote de arquivos, precisamos incluir a adição do nome do arquivo e da escala no loop de transcodificação S57->geopackage. Para fazer isso, usaremos um script python que adiciona os atributos enc_chart e scale às tabelas criadas pelo comando ogr2ogr e as preenche com o nome do arquivo correspondente e uma escala definida pelo produtor. O script também exclui todas as tabelas vazias no geopackage de importação.
Agradecemos a Janez AVSEC pela ideia de usar a tabela DSID para recuperar a escala de compilação de dados.
A chamada para cada arquivo .bat tem o seguinte formato:
.\fichier.bat “diretório que contém os arquivos S57 a serem carregados” “nome do geopacote de importação”
Por exemplo:
.\mass_load_s57_polys.bat “c:/enc_db/enc000” “c:/enc_db/polys.gpkg”
Os arquivos .bat para essa tarefa são os seguintes:
Para importar tabelas de tipos de pontos:
@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
Você perceberá que, além do comando ogr2ogr para importar camadas de pontos, importamos duas tabelas sem geometria: DSID e C_AGGR. A tabela DSID contém dois atributos com o nome do arquivo S57 e a escala de compilação de dados. Essa tabela será usada pelo script Python (update_geopackage_dsid.py) para inserir o arquivo S57 original e a escala de compilação em cada tabela do geopackage.
Para importar tabelas do tipo linha:
@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 tabelas 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
Os arquivos .bat contêm o comando ogr2ogr e uma linha para execução de um script Python, update_geopackage_dsid.py. Certifique-se de alterar o caminho desse script no arquivo .bat para o caminho de sua escolha. Esse script adiciona os dois atributos enc_chart e scale a todas as tabelas do geopackage.
O código python usado é 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_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)Você poderia criar um único arquivo .bat para todos os três processos, mas precisaria verificar os resultados para cada tipo de geometria.
Você pode baixar os arquivos .bat e o código Python clicando aqui.
Continua, independentemente do número de arquivos carregados
Clonagem e/ou adição de tabelas importadas ao geopacote ENC
Para o primeiro carregamento do S57, o script Python a seguir cria todas as tabelas presentes nos três geopacotes de importação. Para cargas subsequentes, se a tabela já existir no geopacote ENC, os registros da tabela no geopacote de importação serão adicionados aos registros já presentes no ENC. Se, por outro lado, a tabela ainda não existir, ela será criada e preenchida com registros do geopacote de importação. No arquivo do geopacote ENC, as tabelas são prefixadas com pt_, li_ e pl_ para indicar seu tipo de geometria.
A clonagem é realizada usando o seguinte 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()
Você pode fazer o download do arquivo .py aqui: clone_or_append_tables_with_prefix.py.
Gerenciamento de duplicatas
Quando dois mapas ENC são sobrepostos, as entidades presentes na zona de sobreposição serão inseridas tantas vezes quantas forem as sobreposições. Para evitar duplicatas e mais, o banco de dados deve ser limpo após cada importação de um arquivo S57.
Mas há dois tipos diferentes de duplicação. Por um lado, as áreas em que dois mapas de escala semelhante se sobrepõem produzirão duplicatas “reais” que precisam ser limpas. Por outro lado, a sobreposição de mapas de escalas muito diferentes precisa ser tratada de forma diferente. O mapa menos detalhado (escala pequena) terá uma seleção de dados do mapa mais detalhado (escala grande). Nesse segundo caso, não é aconselhável remover essas informações.
O código Python abaixo removerá as duplicatas “reais”. Para outras duplicatas, recomendamos o uso de filtros QGis. Dependendo de suas necessidades, você pode usar o atributo “scale” (escala) que adicionamos ao excluir as tabelas vazias para definir quais dados são exibidos no seu mapa e quais não são.
Nesse exemplo, somente os recursos correspondentes a uma escala entre 12500 e 45000 serão exibidos em seu mapa.
Essa operação corresponde à filtragem de uma camada, mas… o que você faz quando tem 130 camadas exibidas? Para evitar que você tenha de executar essa operação 130 vezes, aqui está um código Python que permite filtrar todas as camadas do projeto:
#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'")
E outro para remover o filtro de todas as camadas exibidas:
# 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()}")
Faça o download dos códigos python aqui.
O código a seguir cria uma tabela temporária (temp_table) para cada tabela no geopackage ENC. Em seguida, ele preenche essa tabela com os registros separados com base no LNAM. Em seguida, ele exclui a tabela original e renomeia a tabela temporária com o nome da tabela original.
Esse tratamento se aplica a todas as tabelas S57, exceto
- as tabelas de luzes (LIGHTS) e sondagens batimétricas (SOUNDG), em que o atributo LNAM não é útil. Na verdade, cada setor de incêndio tem o mesmo LNAM e name_rcid, mas corresponde a uma cor e a um setor geográfico diferentes. Da mesma forma, o LNAM e o name_rcid das sondagens batimétricas são os mesmos para todos os pontos de sondagem da entidade multiponto original. No primeiro caso, precisamos procurar cores e setores duplicados para o mesmo name_rcid e, para as sondagens, precisamos encontrar as sondagens com o mesmo XY em sua geometria.
- Tabelas SEAARE, LNDRGN e ADMARE que cobrem vastas áreas e que são divididas de acordo com a área de cobertura do mapa, mas mantendo o mesmo LNAM. Nesse caso, mesmo que o LNAM seja o mesmo, os polígonos não têm a mesma forma e á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()
Você pode fazer download do arquivo .py neste link: supprime_doublons.py. Não se esqueça de alterar os caminhos e os nomes dos arquivos para que correspondam aos seus. Observe também a linha :
if table_name not in [“pl_SEAARE”, “pl_LNDRGN”, “pl_ADMARE”, “layer_styles”, “natsurf”,”DSID”,”C_AGGR”]:
Em princípio, você terá apenas as tabelas prefixadas e as tabelas layer_styles e natsurf. Essa linha impede que pl_SEAARE, pl_LNDRGN e pl_ADMARE sejam processadas, mas, mais importante, impede que as tabelas layer_styles e natsurf que usamos para simbologias sejam esvaziadas. Se você precisar criar tabelas diferentes daquelas do arquivo s57 no geopackage, será necessário levar em conta os nomes dessas tabelas nesse script, adicionando-os a essa linha de exceções.
O projeto do banco de dados S57 com o Geopackage chegou ao fim. Agora estamos começando a implementar um procedimento equivalente usando um banco de dados PostgreSQL/Postgis. Por favor, ajude-nos a concluir esse projeto!
