Não re-inventar a roda

Recentemente eu escrevi um post sobre o processamento de camadas raster com python. A ideia principal é que, usando o terminal de python do QGIS, podemos extrair de uma camada raster a matriz de valores (um objeto numpy ndarray) e processar essa matriz como bem entender. Depois do processamento é preciso retornar para o formato georreferenciado, como .tif ou .asc.

Mas essa abordagem é mais apropriada em situações que você deseja realmente fazer operações customizadas sobre os seus dados. Ou garantir que tudo está completamente sob controle, sem margens para dúvidas. Nesse caso, em certo sentido, corre-se o risco de estar re-inventando a roda.

Em outras ocasiões, talvez até mais frequentes, a demanda consiste em apenas processar as camadas raster de forma mais automatizada com os algoritmos disponívels, que nós confiamos cegamente. Isso a programação viabiliza com mais flexibilidade que os tradicionais processos em lote. Para isso, não é preciso realmente extrair a matriz de valores e fazer operações minuciosas em cima dos dados. Deixe isso para as funções já embutidas (built-in) da biblioteca GDAL e OGR.

Como temos trabalhando intensamente com modelos digitais de elevação (MDE) no blog, uma função embutida particularmente útil é a gdal.DEMProcessing(). Nesse post, vou tecer alguns comentários e dicas sobre essa ferramenta.

A função gdal.DEMProcessing()

A função gdal.DEMProcessing(), como o nome em inglês já anuncia, processa modelos digitais de elevação (MDE). Ela é uma função genérica pois executa mais de um tipo de processamento:

• Sombreamento (hillshade);
• Declividade (slope);
• Aspecto (aspect);
• Relevo colorido (color-relief);
• TRI – índice de robustez do terreno;
• TPI – índice de posição topográfica; e
• Rugosidade (roughness).

É preciso fornecer pelo menos três parâmetros obrigatórios para a função

• Parâmetro destName: string com o nome do arquivo de saída (caminho + nome + extensão);
• Parâmetro srcDS: o objeto Dataset do MDE;
• Parâmetro processing: string que aponta para qual processamento deverá ser feito, entre as opções: "hillshade", "slope", "aspect", "color-relief", "TRI", "TPI", "Roughness".

Claro que, dependendo do processamento, teremos parâmetros adicionais. A lista completa de opções pode ser acessada na documentação oficial da biblioteca gdal. < https://gdal.org/python/osgeo.gdal-module.html#DEMProcessing>. Vamos ver aqui uns exemplos.

Alguns exemplos de uso no terminal python

Uma vez aberto o terminal, a primeira coisa é importar a biblioteca gdal:

import gdal

A segunda coisa é apontar para uma variável o nome do arquivo do MDE e a seguir criar o objeto Dataset com a função gdal.Open():

# arquivo de entrada
file_input = 'C:/Data/mde.tif'

# Dataset de entrada
raster_input = gdal.Open(file_input, 0)

Agora podemos nos aventurar com os processos disponíveis!

Aspecto

Vamos primeiro no aspecto. É preciso definir um nome para o arquivo de saída e convocar a função gdal.DEMProcessing():

# arquivo de saida aspecto
file_output = 'C:/Data/aspecto.tif'

# convocar a função gdal.DEMProcessing()
gdal.DEMProcessing(destName=file_output, srcDS=ds_input, processing='aspect', trigonometric=False, zeroForFlat=True, band=1, computeEdges=True, alg='Horn')

Note que incluímos alguns parâmetros adicionais. A maioria desses parâmetros também estão presentes nos outros algoritmos. Os específicos do Aspecto são trigonometric e zeroForFlat, e ambos recebem booleanos (True ou False). Por padrão, ambos são definidos como False. O aspecto é retornado em graus de azimute se trigonometric=False, caso contrário é retornado em graus trigonométricos onde o leste é zero e o norte é 90. O parâmetro zeroForFlat, quando verdadeiro, define como 0 as áreas planas, em vez de nulo (null).

Na verdade, essa expressão em código consiste na versão backend (por trás dos panos) do formulário que preenchemos toda ver que convocamos o algoritmo pelo interface gráfica do QGIS. A mesma coisa ocorre nos outros processamentos.

Declividade

Novamente, definir o nome do arquivo de saída e convocar a função.

# arquivo de saida declividade
file_output = folder + 'slope.tif'

# convocar a função gdal.DEMProcessing()
gdal.DEMProcessing(destName=file_output, srcDS=ds_input, processing='slope', slopeFormat='degree', scale=1, band=1, computeEdges=True, alg='Horn')

Aqui, a novidade é o parâmetro slopeFormat, que pode ser entre 'degree' e 'percent', e o parâmetro scale, que escalona a razãol vertical para horizontal.

Sombreamento

Mesma coisa, mas agora com novos parâmetros: zFactor, azimuth, altitude.

# arquivo de saida sombreamento
file_output = folder + 'hillshade.tif'

# convocar a função gdal.DEMProcessing()
gdal.DEMProcessing(destName=file_output, srcDS=ds_input, processing='hillshade', zFactor=1, azimuth=45, altitude=20, scale=1, band=1, computeEdges=True, alg='Horn')

E daí? Como automatizar?

Se você não tem muito contato com programação talvez esteja achando muito mais simples preencher os formulários da inteface gráfica ou executar os processos em lote do QGIS. E é verdade que em muitas circunstâncias mundanas pode ser realmente exagero programar os processos.

Mas a automatização começa com os ingredientes clássicos da programação, como os laços iterativos, manejo de strings e estruturas condicionais. É o feijão com arroz da programação que permite automatizar as tarefas. Eu vou dar um exemplo (mas existem outras possibilidades!)

Observe que a função hillshade possui dois parâmetros importantes: o azimute do sol (parâmetro azimuth) e a altitude do sol (parâmetro altitude), ambos ângulos em graus. Considere que você quer avaliar 10 configurações de cada ângulo contra 10 configurações do outro ângulo. Isso vai resultar em 100 mapas! Até mesmo um processo em lote clássico irá ficar tedioso e gastar seu tempo nessa circunstância.

Com python a coisa muda de figura. Tudo que você precisa é definir uma sequência de 10 valores para cada parâmetro e manejar a string do arquivo de saída para cada arquivo ter um nome específico. Insira isso dentro de dois laços aninhados e terá uma metralhadora de camadas raster:

import gdal

# definir pasta de destino
folder = 'C:/Dados/'

# arquivo de entrada
file_input = 'C:/Dados/mde.tif' 
# Dataset de entrada
ds_input = gdal.Open(file_input, 0)

# sequencia tupla de azimutes
azm = (0, 2, 4, 6, 8, 10, 12, 14, 16, 18)
# seququencia tupla de altitudes
alt = (0, 10, 20, 30, 40, 50, 60, 70, 80, 90)

# lacos iterativos aninhados
for i in range(0, len(azm)):
      for j in range(0, len(alt)):  
              # obter os parametros locais
              lcl_azm = azm[i] # azimute local
              lcl_alt = alt[j] # altitude local
              # arquivo de saida declividade (concatenacao de strings)
              file_output = folder + 'hillshade_' + '_azm' + str(lcl_azm) + '_alt' + str(lcl_alt) + '.tif'
              # convocar a funcao gdal.DEMProcessing()
              gdal.DEMProcessing(destName=file_output, srcDS=ds_input, processing='hillshade', zFactor=1, azimuth=lcl_azm, altitude=lcl_alt, scale=1, band=1, computeEdges=True, alg='Horn')