A função GDAL Translate

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 

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair /  Alterar )

Foto do Google

Você está comentando utilizando sua conta Google. Sair /  Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair /  Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair /  Alterar )

Conectando a %s

%d blogueiros gostam disto: