Esse post objetiva trazer mais uma contribuição no tópico de processamento de camadas raster no terminal de Python do QGIS. O primeiro post dessa série aborda de forma geral como acessar a matriz de dados bruta das camadas raster e processá-la como bem entender. O segundo post aborda que podemos usar funções embutidas da biblioteca gdal para realizar operações, como a função gdal.DEMProcessing() para processamento de Modelos Digitais de Elevação. Agora, trago algumas considerações e “receitas” de como usar a função gdal.Translate(), uma função de grande utilidade no manejo de camadas raster.
A função gdal.Translate()
Manejar camadas raster é uma realidade para quem trabalha com geoprocessamento. Nada mais comum, por exemplo, do que precisar recortar uma extensão de uma camada raster. Ou converter para outro formato. Para isso, e outras funcionalidades, existe a função gdal.Translate(). Essa função está disponível na caixa de ferramentas do QGIS como “Translate”, em inglês, e “Converter”, em português. Mas não somente: quando usamos a calculadora raster, ela está lá, nos bastidores. E também quando usamos a ferramenta de “recortar raster pela extensão”.
Opções, muitas opções
As funções gdal tendem a ser amplas, e por isso possuem muitos parâmetros que precisam ser informados. Por isso, existe uma maneira mais “humana” de programar elas: aplicando as opções em uma única variável e depois a passando a variável de opções para um parâmetro chamado de options. O trabalho é praticamente o mesmo, mas o código sem dúvida fica mais legível por seres humanos (normais).
Para a função gdal.Translate(), temos a função auxiliar gdal.TranslateOptions(), que retorna um objeto de opções. A situação é análoga para outras funções da biblioteca gdal. De fato, são muitas opções para apresentar em um único post, mas elas estão documentadas (não muito bem, é verdade) aqui: https://gdal.org/python/osgeo.gdal-module.html#TranslateOptions .
Recortando uma extensão
O exemplo mais prático da função gdal.Translate() é o recorte de uma extensão retangular na camada raster original, exportando uma nova camada menor. Isso requer que seja informado como opção a “janela” de recorte pelas coordenadas X e Y do canto superior esquerdo e X e Y do canto inferior direito. Essas coordenadas precisam ser passadas como uma lista de quatro números para o parâmetro projWin, na ordem X, Y, superior, inferior. Além disso, é preciso dizer qual é a banda (parâmetro bandList) e o sistema de referência da extensão de recorte (parâmetro projWinSRS). A receita abaixo deixa tudo mais óbvio:
import gdal
# abrir camada de entrada
src_f = 'C:/dados/../camada_raster_entrada.tif'
src_ds = gdal.Open(src_f, 0) # carrega a camada raster em um objeto "dataset"
# definir o nome da camada de saída
out_f = 'C:/dados/../camada_raster_saida.tif' # incluir a extensão de saída
'''
[opcional] usar a função auxiliar gdal.TranslateOptions() para configurar as opções de tradução
ver em https://gdal.org/python/osgeo.gdal-module.html#TranslateOptions todas as opções
'''
# aqui será extraído um recorte da seguinte extensão do SRC:
ulx = 342505 # exemplo de coordenada X projetada do canto superior esquerdo
uly = 6736694 # exemplo de coordenada Y projetada do canto superior esquerdo
lrx = 356509 # exemplo de coordenada X projetada do canto inferior direito
lry = 6730179 # exemplo de coordenada Y projetada do canto inferior direito
# texto padronizado para definir o sistema de projeção:
srs_proj = 'EPSG:31982, SIRGAS 2000 / UTM zone 22S'
# definição das opções em um objeto de opções
options = gdal.TranslateOptions(format='GTiff', bandList=[1], projWin=[ulx, uly, lrx, lry], projWinSRS=srs_proj)
# convocar a função Translate e passar o objeto opts no parâmetro 'options'
out_ds = gdal.Translate(destName=out_f, srcDS=src_ds, options=options)
# fechar os datasets:
src_ds = None
out_ds = None
Recortando a matriz
Uma alternativa ao recorte da camada pela extensão é o recorte da matriz de dados. Com isso, o recorte é feito por uma extensão retangular definida em píxeis, não no sistema de referência de coordenadas. A “janela”, assim, é definida no parâmetro srcWin pela coordenada I, J na matriz do píxel superior esquerdo e pela largura e altura em píxeis (o que também resulta em uma lista de quatro números).
import gdal
# abrir camada de entrada
src_f = 'C:/dados/../camada_raster_entrada.tif'
src_ds = gdal.Open(src_f, 0) # carrega a camada raster em um objeto "dataset"
# definir o nome da camada de saída
out_f = 'C:/dados/../camada_raster_saida.tif' # incluir a extensão de saída
'''
[opcional] usar a função auxiliar TranslateOptions para configurar as opções de tradução
ver em https://gdal.org/python/osgeo.gdal-module.html#TranslateOptions todas as opções
'''
# aqui será extraído um recorte de 100 x 100 píxeis a partir da origem da matriz (canto superior esquerdo)
uli = 0 # exemplo de coordenada I (em píxel) do canto superior esquerdo
ulj = 0 # exemplo de coordenada J (em píxel) do canto superior esquerdo
width = 100 # exemplo de largura em píxel
height = 100 # exemplo de altura em píxel
# definição das opções:
options = gdal.TranslateOptions(format='GTiff', bandList=[1], srcWin=[uli, ulj, width, height])
# convocar a função Translate e passar o objeto 'options'
out_ds = gdal.Translate(destName=out_f, srcDS=src_ds, options=options)
# fechar os datasets:
src_ds = None
out_ds = None
Exportando para outro formato
O que realmente confere o nome “translate” (traduzir, em inglês) é que a gdal.Translate() realizar as conversões entre formatos de arquivos de camadas. Não tenho como apresentar todas as opções, mas o formato deve ser informado por uma string (texto) reservada para o parâmetro format. Fizemos isso nos exemplos acima explicitamente para manter formato em GTiff (geotiff). O formato .asc, já mencionado no primeiro post dessa série, é codificado pela string AAIGrid:
import gdal
# abrir camada de entrada
src_f = 'C:/dados/../camada_raster_entrada.tif'
src_ds = gdal.Open(src_f, 0) # carrega a camada raster em um objeto "dataset"
# definir o nome da camada de saída
out_f = 'C:/dados/../camada_raster_saida.asc' # incluir nova extensão de saída
# aqui a camada será convertida para o formato "ASC"
options = gdal.TranslateOptions(format='AAIGrid', bandList=[1])
# convocar a função Translate e passar o objeto 'options'
out_ds = gdal.Translate(destName=out_f, srcDS=src_ds, options=options)
# fechar os datasets:
src_ds = None
out_ds = None
Exportando as bandas de uma camada
É possível acessar a matriz de cada banda em uma camada raster e exportar elas separadamente usando a receita que eu apresentei no primeiro post da série. Mas aquele método é melhor caso você realmente precise processar a matriz das bandas. Se o problema é só extrais em arquivos separados, a função gdal.Translate() deixa isso muito mais simples:
import gdal
# abrir camada de entrada
src_f = 'C:/dados/../camada_raster_entrada.tif'
src_ds = gdal.Open(src_f, 0) # carrega a camada raster em um objeto "dataset"
# iterar ao longo das bandas (aqui 3 bandas):
for i in range(1, 4):
# definir o nome da camada de saída
out_f = 'C:/dados/../camada_raster_saida_banda{}.tif'.format(str(i)) # nome customizado com o índice da banda
# selecionar a banda pelo índice
options = gdal.TranslateOptions(format='GTiff', bandList=[i])
# convocar a função Translate e passar o objeto 'options'
out_ds = gdal.Translate(destName=out_f, srcDS=src_ds, options=options)
# fechar os datasets:
out_ds = None # dentro do laço for
src_ds = None # fora do laço for
Brinde: processamento de uma pasta inteira
O verdadeiro poder da programação está na automatização do fluxo de trabalho de uma forma incomparável, especialmente com as estruturas de iteração como o laço for. E um aspecto especial do Python são suas funções embutidas que fazem tudo muito mais simples. Por exemplo, com a função nativa os.listdir() podemos carregar em uma lista o nome de todos os arquivos em uma pasta. O caminho da pasta deve ser passado como o parâmetro da função. Assim, podemos realizar uma batelada da função gdal.Translate() sobre uma pasta inteira. No exemplo abaixo, a operação é um recorte em matriz. Nesse caso, é preciso garantir manualmente que os arquivos da pasta sejam apenas camadas raster com o tamanho apropriado para a operação. Se você já domina Python, poderá elaborar um fluxo um pouco mais sofisticado para “validar” as camadas de entrada, mas isso deixo como um desafio.
import gdal
import os
# definir as pastas de insumos e produtos
src_dir = 'C:/dados/../pasta_insumos' # a pasta precisa conter apenas as camadas raster para processamento
out_dir = 'C:/dados/../pasta_produtos' # é preciso criar essa pasta
# carregar em uma lista o nome de todas as as camadas
insumos = os.listdir(src_dir)
# iterar ao longo das bandas (aqui 3 bandas):
for rst in insumos:
# abrir camada de entrada
src_f = src_dir + '/' + rst # concatenar o caminho inteiro
src_ds = gdal.Open(src_f, 0) # carrega a camada raster em um objeto "dataset"
# definir o nome da camada de saída
out_f = out_dir + '/' + rst # nome customizado com o sufixo 'out_'
# definir as opções de tradução (aqui será o recorte de 100 x 100)
options = gdal.TranslateOptions(format='GTiff', bandList=[1], srcWin=[400, 1400, 300, 300])
# convocar a função Translate e passar o objeto 'options'
out_ds = gdal.Translate(destName=out_f, srcDS=src_ds, options=options)
# fechar os datasets:
out_ds = None
src_ds = None
Mapas Abertos 