Processamento em lotes – automatizando o uso de ferramentas do GDAL e do SAGA GIS por meio de scripts Shell do tipo Bash

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

Bem-vindos ao meu blog!

Gosto de começar com um aviso - posso ser uma pesquisadora desta área, mas isso não significa que tudo que eu faço ou escrevo aqui funcionará para você, no seu computador e com suas configurações e versões de pacotes. Me eximo de qualquer responsabilidade se você perder dados ou estragar sua instalação. Também não autorizo cópias do meu conteúdo.

Já faz algum tempo desde que eu apresentei para vocês a versão de linha de comando do SAGA GIS. Uma das vantagens de usar SAGA ou qualquer outro software na interface de linha de comando (CLI) é a facilidade de fazer processamentos em lotes. Além do SAGA, outro software de SIG que eu normalmente uso pela CLI é o GDAL.

Hoje irei mostrar para vocês como usar um script do tipo Bourne Again Shell (Bash) para fazer processamento em lotes de arquivos, fornecendo um exemplo em que o GDAL e o SAGA GIS são usados. Existem outras maneiras de fazer processamento em lotes, por exemplo, através do Python, mas não vou falar delas no post de hoje. Estou, portanto, mostrando uma das maneiras mais simples e primitivas de fazer geoprocessamento em lotes, uma que não requer interpretador de linguagem Python.

Certifique-se de que você pode executar um script Bash em sua máquina

Uma das maneiras mais diretas, e por outro lado, uma das mais complicadas, de fazer processamento em lotes é usando scripts Bash. Bash é um tipo de Shell, o que significa que é executado diretamente em uma interface de usuário que acessa os serviços do sistema operacional. Então, isso significa que scripts Bash podem ser executados diretamente no Prompt de Comando do Windows? Na verdade, não. Veja abaixo como usar Bash no Windows.

Usuários Windows

Eu atualmente uso o Windows e uso scripts Bash, então dá para fazer isso, sim. Para conseguir usar Bash no Windows, instalei o Git para Windows (https://gitforwindows.org/) que vem com a funcionalidade Git Bash Here. Mas existem outras maneiras. Uma delas é usar o Windows Subsystem for Linux - subsistema do Windows para Linux (https://docs.microsoft.com/pt-br/windows/wsl/install-win10).

No Windows 10, se você tem Git para Windows, você pode clicar com o botão direito dentro de uma pasta no Windows Explorer, e deve aparecer a opção “Git Bash Here”. Se você clicar, um prompt semelhante a um terminal do Ubuntu será aberto, e a maioria dos comandos do terminal do Ubuntu também funcionará lá. Caso você não esteja muito familiarizado com esses comandos, alguns deles são:

• ls: lista os arquivos na pasta

• ls *.sh: lista todos os arquivos com extensão .sh

• cd ..: vai para a pasta acima

• cd nome_da_pasta: vai até a pasta

• ./nome_arquivo.sh: executa este arquivo executável

• bash ./nome_arquivo.sh: executa este arquivo executável como Bash

• rm nome_arquivo: remove (deleta) o arquivo

Se você tiver um Shell habilitado para Bash funcionando, você pode começar a construir seu script de Shell para geoprocessamento em lotes.

Usuários Linux (Ubuntu, Fedora, Debian, etc.)

Nenhuma ação necessária.

Faça um teste

Antes de construir um script inteiro e descobrir que na verdade você não precisava daquele comando, mas de outro, ou que o aplicativo de que você precisa não está na pasta que você pensava, você deve primeiro testar os comandos separadamente.

Os scripts do Bash executam cada linha separadamente sem “validar” ou compilar o código, portanto, executar cada linha separadamente em um terminal (ou usando o Git Bash Here) deve ter um resultado semelhante.

Esquematize

Você provavelmente quer fazer processamento em lotes porque tem muitos arquivos para processar. Leve o primeiro deles para outra pasta para testar os seus comandos nele antes de fazer para todos. Em seguida, esquematize quais são os comandos que você vai precisar executar.

Exemplo:

Preciso calcular a declividade do terreno, baseado em um Modelo Digital de Elevação (MDE), dentro dos limites de um município, que está delimitado por um polígono salvo em um shapefile.

Isso significa que tenho um raster grande que preciso reprojetar, cortar com base em um vetor de máscara, e calcular a declividade.

Para reprojetar, você pode usar a ferramenta GDAL Warp. Para cortá-lo com base em um vetor de máscara, você também pode usar GDAL Warp, com a opção -cutlines, se seu GDAL for de uma das versões mais recentes. Para calcular a declividade, você pode usar o Módulo SAGA Slope, Aspect, Curvature.

No entanto, essa não é a melhor ordem para os cálculos, infelizmente. É melhor primeiro reprojetar para um sistema de coordenadas projetadas, calcular a declividade e, em seguida, recortar pela camada de máscara. Desta forma, garantimos que a região das divisas do município não terá valores de declividade espúrios. Se a projeção original do MDE é um sistema de coordenadas projetadas e ele está sendo projetado em um sistema de coordenadas geográficas (do tipo lat long), devemos calcular a declividade antes de reprojetar, para evitar erros do cálculo na declividade.

As ferramentas de que vou precisar são:

• GDAL Warp

• Módulo SAGA GIS Slope, Aspect, Curvature

• GDAL Warp (novamente)

Encontre as ferramentas a serem usadas em seu sistema, ou baixe-as

Por exemplo, se você tem o QGIS instalado, o GDAL geralmente está dentro da pasta em que o QGIS está instalado, na pasta “bin”. Encontre as ferramentas GDAL a serem usadas lá dentro.

Entre na pasta do GDAL e verifique a versão do mesmo digitando:

./gdalinfo.exe --version

Para referência, a minha versão é:

GDAL 3.1.4, released 2020/10/20

No Windows, a versão da linha de comando SAGA GIS pode ser encontrada onde quer que ela tenha sido baixada e descompactada. No Linux, o SAGA GIS pode ser baixado no repositório relevante ou compilado pelo usuário. Expliquei como usar a interface de linha de comando SAGA GIS em um post anterior.

Por exemplo, em um computador Windows:

GDAL Warp é gdalwarp.exe dentro da pasta QGIS 3.18/bin.

SAGA GIS Module Slope, Aspect, Curvature pode ser chamado dentro da pasta em que saga_cmd.exe está, pelo comando ./saga_cmd ta_morphometry 0. Veja mais em SAGA-GIS Module Library Documentation.

Crie um código de teste

Digamos que o raster original é original_raster.tif e é salvo em caminho_arquivo_original. Também salvo em caminho_arquivo_original está municipio_mascara.shp

caminho_gdal é o caminho onde GDAL pode ser encontrado e caminho_saga é o caminho onde saga_cmd pode ser encontrado.

Para o exemplo, pode-se construir algo como o disposto abaixo.

Para entrar na pasta onde o GDAL fica:

cd caminho_gdal

Para executar GDAL Warp (verifique as entradas possíveis aqui https://gdal.org/programs/gdalwarp.html)

./gdalwarp -s_srs original_CRS -t_srs destino_CRS -r bilinear caminho_arquivo_original/original_raster.tif caminho_arquivo_original/reprojetado.tif

Agora, devemos calcular a declividade do raster gerado. Primeiro, entre na pasta do SAGA e, em seguida, execute o comando.

cd caminho_saga
saga_cmd ta_morphometry 0 -ELEVATION caminho_arquivo_original/reprojetado.tif -SLOPE caminho_arquivo_original/calc_slope.tif -ASPECT apagar_este_arquivo_depois.tif

Observe que a unidade padrão para a declividade é radianos. Para alterar para graus, adicione a opção “-UNIT_SLOPE 1” ao final do comando.

Volte para a pasta do GDAL e execute GDAL Warp novamente:

cd caminho_gdal
./gdalwarp -r bilinear -cutline caminho_arquivo_original/municipio_mascara.shp -cl municipio_mascara -crop_to_cutline caminho_arquivo_original/calc_slope.tif caminho_arquivo_original/slope_mascara.tif

E, para este caso, slope_mascara.tif é o arquivo final. Abra-o em um programa de geoprocessamento e verifique se é o que você deseja. Se for, prossiga.

Envolva os comandos em um loop Bash

É hora de construir este código em loop para que ele seja um processamento em lote.

Crie o script Bash

Crie um novo documento de texto no Bloco de Notas. Salve-o como qualquer_nome.sh

Se você estiver no Windows, inicie o código com:

#!/bin/bash

Assim, você pode chamar o executável no Windows. No Linux, para poder executar seu script de Shell, você precisará alterar as permissões do arquivo para permitir que ele seja executado.

Nas demais linhas, “#” denota o início de um comentário.

Um loop do tipo “for” no Bash se parece com isto:

for i in {1..100}
do
	# insira os comandos aqui
done

Nesse exemplo, você está fazendo um loop na variável i.

Para imprimir algo na tela do terminal, o comando é “echo”. Teste isto:

for i in {1..2}
do
   echo "---"
   echo "iteracao $i"
   echo "---"
done

E execute-o simplesmente escrevendo no Git Bash Here

./seu_script.sh 

Ele deve mostrar as iterações 1 e 2 na tela.

---
iteracao 1
---
---
iteracao 2
---

Se isso der certo, então você sabe que está indo na direção certa.

Construa os nomes dos arquivos.

Normalmente, os arquivos para processamento em lote são gerados automaticamente e por isso possuem todos o mesmo nome seguido de um número, o que facilita esse processo.

ATUALIZAÇÃO (17.09.22): se os seus arquivos não estão nomeados desta maneira, mas estão todos na mesma pasta, preparei um snippet que pode lhe ajudar. Encontre-o no fim deste post.

Com base em seu script escrito anteriormente, adicione o caminho para seus rasters de entrada no loop:

for i in {1..2}
do
  echo "---"
   echo "iteracao $i"
   echo "---"
   folder_orig="pasta_arquivo_original/orig_raster"
   extensao=".tif"
   compose="$folder_orig$i$extensao"
   echo $compose
done

E veja se ele mostra corretamente no terminal o caminho para seus arquivos de entrada. Faça ajustes, se necessário. Para o exemplo acima, a saída é:

---
iteracao 1
---
pasta_arquivo_original/orig_raster1.tif
---
iteracao 2
---
pasta_arquivo_original/orig_raster2.tif

Adicionando os comandos e pré-definições de pastas

Crie os nomes das saídas da mesma forma que fez com a entrada.

As declarações de variáveis que não mudam podem estar acima do loop:

pasta_intermediarios="pasta_arquivo_original/intermediarios/"
arquivo_warp1="warp1_"
arquivo_slope1="slope1_"
arquivo_final="final_"

E os caminhos completos podem ser compostos dentro do loop:

composicao="$pasta_intermediarios$arquivo_warp1$i$extensao"

Adicione isso e os comandos de geoprocessamento ao seu script:

pasta_intermediarios="pasta_arquivo_original/intermediarios/"
arquivo_warp1="warp1_"
arquivo_slope1="slope1_"
arquivo_final="final_"
for i in {1..2}
do
   echo "---"
   echo "iteracao $i"
   echo "---"
   folder_orig="pasta_arquivo_original/orig_raster"
   extensao=".tif"
   caminho_arquivo_original="$folder_orig$i$extensao"
   caminho_arquivo_warp="$pasta_intermediarios$arquivo_warp1$i$extensao"
   caminho_arquivo_slope="$pasta_intermediarios$arquivo_slope1$i$extensao"
   caminho_arquivo_final="$pasta_intermediarios$arquivo_final$i$extensao"
   #GDAL Warp
   cd caminho_gdal
   ./gdalwarp -s_srs sistema_coord_original -t_srs sistema_coord_destino -r bilinear $caminho_arquivo_original $caminho_arquivo_warp
   #SAGA Modulos
   cd caminho_saga
   # o arquivo da orientacao da encosta (aspect) pode ser reescrito 
   # em todos os loops porque ele nao e importante para esse processamento.
   ./saga_cmd ta_morphometry 0 -ELEVATION $caminho_arquivo_warp -SLOPE $caminho_arquivo_slope -ASPECT delete_this_file_later.tif
   # GDAL Warp
   cd caminho_gdal
   ./gdalwarp -r bilinear -cutline pasta_arquivo_original/municipality_mask.shp -cl municipality_mask -crop_to_cutline $caminho_arquivo_slope $caminho_arquivo_final
done

Teste inicialmente com poucas iterações

Antes de executar o processamento em lotes para seus 100 (1000?) rasters, execute-o para os dois primeiros. Se funcionar sem problemas, aí sim, execute o script de shell para o restante dos rasters. É melhor descobrir um erro depois de 30 segundos, do que depois de duas horas. Ajuste as iterações no loop “for”.

Para executar o script, digite no terminal:

./seu_script.sh 

Substitua seu_script pelo nome que você deu ao seu script.

Se não funcionar, você pode tentar:

bash seu_script.sh

E é assim que eu faço processamento em lotes de geoprocessamento usando Bash!

Exemplos, para o caso apresentado, de como os arquivos devem ficar em cada etapa do processamento

Testei esse script com sucesso com dados de elevação do CIAT-CSI SRTM, uma versão do Modelo Digital de Elevação da missão da NASA Shuttle Radar Topographic Mission (SRTM), com vazios preenchidos. O polígono utilizado foi desenhado no QGIS.

Para o exemplo apresentado, os dados originais usados foram:

Após a projeção, o MDE não muda significativamente:

Declividade calculada:

Após o corte pela máscara:

Essa foi a saída que obtive no terminal:

---
iteracao 17
---
Creating output file that is 5875P x 6281L.
Processing initial.tif [1/1] : 0Using internal nodata values (e.g. -32768) for image initial.tif.
Copying nodata values from source initial.tif to destination results_warp1_17.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
____________________________

   #####   ##   #####    ##
  ###     ###  ##       ###
   ###   # ## ##  #### # ##
    ### ##### ##    # #####
 ##### #   ##  ##### #   ##
____________________________

SAGA Version: 7.9.0 (64 bit)

____________________________
library path: ...
library name: ta_morphometry
library     : ta_morphometry
tool        : Slope, Aspect, Curvature
identifier  : 0
author      : O.Conrad (c) 2001
processors  : 4 [4]
____________________________

loading: results_warp1_17

100%

Parameters


Grid System: 89.069219; 5875x 6281y; 81361.296189x 7229119.462180y
Elevation: results_warp1_17
Slope:
Aspect:
General Curvature: <not set>
Profile Curvature: <not set>
Plan Curvature: <not set>
Tangential Curvature: <not set>
Longitudinal Curvature: <not set>
Cross-Sectional Curvature: <not set>
Minimal Curvature: <not set>
Maximal Curvature: <not set>
Total Curvature: <not set>
Flow Line Curvature: <not set>
Method: 9 parameter 2nd order polynom (Zevenbergen & Thorne 1987)
Unit: radians
Unit: radians


101%Saving grid: results_slope1_17.tif...
Export GeoTIFF


Parameters


Grid System: 89.069219; 5875x 6281y; 81361.296189x 7229119.462180y
Grid(s): 1 object (results_slope1_17)
File: results_slope1_17.tif
Creation Options:

Band 1

100%
100%okay
Saving grid: delete_this_file_later.tif...
Export GeoTIFF


Parameters


Grid System: 89.069219; 5875x 6281y; 81361.296189x 7229119.462180y
Grid(s): 1 object (delete_this_file_later)
File: delete_this_file_later.tif
Creation Options:

Band 1

100%
100%okay
Creating output file that is 3633P x 4339L.
Processing results_slope1_17.tif [1/1] : 0Using internal nodata values (e.g. -99999) for image results_slope1_17.tif.
Copying nodata values from source results_slope1_17.tif to destination results_final_17.tif.
...10...20...30...40...50...60...70...80...90...100 - done.

Extras:

• Por que não usar Python?

Se você está familiarizado com Python e prefere usá-lo, pode usar, sim! No entanto, o objetivo da minha postagem de hoje é mostrar como fazer processamento em lote com scripts Bash!

• Tenho uma parte do código em Python, que preciso executar, mas gostaria de fazer o processamento em lote com o Bash, é possível?

Sim, é possível! Você pode chamar seu código Python personalizado como faria normalmente no terminal:

python arquivo_python.py argumento1 argumento2

Seu arquivo de código Python deve definir ações para quando o arquivo for chamado. Você pode fazer isso inserindo este snippet no final do seu arquivo Python:

if __name__ == '__main__':
    nome_da_funcao(sys.argv[1], sys.argv[2])

lembre-se de adicionar ao início do seu arquivo:

import sys

• Ok, gerei mais de 200 arquivos. Como posso fazer para abri-los no QGIS?

Eu escrevi um post sobre isso na semana passada! Nele, eu explico como eu faço para abrir muitos arquivos ao mesmo tempo usando o Console Python no QGIS 3.18. Confira!

• Posso excluir os arquivos intermediários automaticamente, dentro do script?

Sim, você pode! Às vezes, o processamento em lote pode ocupar muito espaço em disco com os arquivos auxiliares. Exclua-os como faria em um terminal normal, mas use as variáveis com os seus nomes, dentro do loop, depois de usar o arquivo para o processamento seguinte.

rm $caminho_arquivo_warp 
rm $caminho_arquivo_slope

• ATUALIZAÇÃO (17.09.22) E se os meus arquivos não estão numerados, e não seguem nenhum padrão na nomenclatura? Ainda posso processá-los em bateladas utilizando o Bash?

Sim, se os seus arquivos estão todos na mesma pasta, você pode!

Primeiramente, você deve listar os arquivos com a extensão desejada (geralmente .tif) dentro da pasta. Em seguida, realizar o loop nesta lista de arquivos.

Exemplo (baseado nesta resposta do Stack Overflow):

folder_intermediate="D:/.../intermediate/"
file_warp1="warp1_"
yourfilenames=`ls *.tif`
count=1
for entry in $yourfilenames
do
   echo "---"
   echo "Iteration $count"
   echo "---"
   echo $entry
   folder_orig="D:/..."
   extension=".tif"
   path_file_orig="$folder_orig$entry"
   path_file_warp="$folder_intermediate$file_warp1$count$extension"
   count=`expr $count + 1`
   #GDAL Warp
   cd  path_gdal
   ./gdalwarp -s_srs the_original_CRS -t_srs the_target_CRS -r bilinear $path_file_orig $path_file_warp
done

Também adicionei um contador, para que os arquivos de saída tenham nomes diferentes e não sobrescrevam-se a cada iteração.