Resolvido: Python Console do QGIS 3.20 mostra um erro ao executar a ferramenta Clip Vector by Extent

Colegas pesquisadores e entusiastas de SIG de código-livre,

Bem-vindos ao meu blog!

Gostaria de começar com um aviso - posso ser uma pesquisadora desta área, mas isso não significa que tudo o que faço ou escrevo aqui funcionará para você, em suas próprias configurações de desktop e versões de packages. Não tenho nenhuma responsabilidade se você perder dados ou estragar sua instalação. Eu também não autorizo nenhuma cópia do meu conteúdo.

Hoje, estou escrevendo sobre um erro que obtive ao usar o Console Python do QGIS 3.20. Eu estava tentando recortar uma camada vetorial preenchida com pontos que tinha acabado de criar. Estava cortando ela com base em uma extensão do tipo $xmin, xmax, ymin, ymax$, por meio da ferramenta Clip Vector by Extent, do GDAL.

O erro

Primeiro, executei a seguinte linha, para definir os parâmetros da ferramenta:

paramgdal= {'INPUT': points_layer, 'EXTENT':extent_final, 'OUTPUT':output_path}

As variáveis “extent_final” e “output_path” foram previamente setadas. “points_layer” é a camada vetorial (no caso, um shapefile) em que o universo de pontos dos quais quero recortar uma parte estão armazenados. Ele foi carregado anteriormente, usando o seguinte comando:

points_layer = QgsVectorLayer(path_to_file, 'points_loaded', 'ogr')

Até aqui, tudo bem.

Mas, depois de definir os parâmetros, executei:

processing.run("gdal:clipvectorbyextent", paramgdal)

e aí, este erro apareceu:

ERROR 1: Attempt to write non-multipoint (POINT) geometry to multipoint shapefile. 
ERROR 1: Unable to write feature 0 from layer points_layer. 
ERROR 1: Terminating translation prematurely after failed translation of layer points_layer (use -skipfailures to skip errors)

E o layer resultante não foi carregado (porque não foi criado corretamente).

Interpretação do problema

Relendo a mensagem de erro, percebi que o GDAL, por algum motivo, tenta salvar a camada gerada como um shapefile do tipo Multipoint, ao invés de, simplesmente, um shapefile do tipo Point. Por esse motivo, ele não consegue gravar nenhum dos pontos que estavam no layer original dentro do novo shapefile, porque os tipos não correspondem.

Solução

Verificando a documentação do QGIS, podemos notar que existem quatro possíveis parâmetros para executar gdal:clipvectorbyextent. Um deles, “OPTIONS”, é opcional, e se refere às opções de criação do GDAL. Como o algoritmo Clip Vector by Extent é baseado em um utilitário do GDAL chamado ogr2ogr, as opções disponíveis para o algoritmo original também podem ser utilizadas para construir o novo layer, no caso, no Python Console.

Uma das opções é particularmente útil, neste cenário. A opção “nlt” do ogr2ogr define a geometria do arquivo gerado, o que é exatamente o que estávamos procurando. No caso mostrado aqui, simplesmente adicionei

'OPTIONS': '-nlt point'

ao dicionário Python chamado de $gdalparams$. Isso diz ao GDAL que o vetor a ser criado é um vetor do tipo Points. E isso resolveu meu problema completamente! O fragmento de código mostrado no início desta postagem foi adaptado por mim como:

paramgdal= {'INPUT': points_layer, 'EXTENT':extent_final, 'OPTIONS': '-nlt point', 'OUTPUT':output_path}
processing.run("gdal:clipvectorbyextent", paramgdal)

E quando executei o fragmento acima, a camada resultante foi gerada corretamente, sem erros ou warnings.

Ou seja, uma solução simples!

Extras

• Por que você simplesmente não usa o GDAL direto por linha de comando, ou clica na ferramenta e a usa diretamente na interface gráfica? Por que você quer usar o Console Python do QGIS?

Porque, às vezes, o console Python é a maneira mais prática de realizar um cálculo. Este trecho de código mostrado aqui é uma pequena parte de um código muito maior. Por esse motivo, não vale a pena mudar toda a interface ou linguagem na qual estamos desenvolvendo o código.

• Carreguei meu arquivo vetorial no Python Console, mas ele não apareceu na lista de camadas do QGIS! O que aconteceu?

Isto não é um problema. O trecho de código usado para carregar a camada vetorial não lista a camada no QGIS, nem a mostra a visualização da camada no mapa. Para fazer isso, você deve executar

QgsProject.instance().addMapLayers([points_layer])

• Não sei o que devo importar antes de usar a ferramenta “gdal:clipvectorbyextent” no Console Python do QGIS.

Você deve importar os algoritmos do QGIS e a ferramenta de processamento. Veja o fragmento de código abaixo:

from qgis.core import *
import processing