Como abrir, editar e salvar arquivos raster (.tif) usando Python

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

Bem-vindos ao meu blog!

Gostaria de começar com um aviso - posso ser 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. Me eximo de responsabilidade se você perder dados ou bagunçar sua instalação. Eu também não autorizo nenhum tipo de cópia do meu conteúdo.

Hoje, vou escrever sobre como abrir, editar e salvar arquivos raster usando Python, se você já estiver familiarizado com Python (principalmente Anaconda).

Talvez você precise executar uma função Python em cada pixel de um raster, ou talvez você deseje fazer operações personalizadas que não seriam possíveis em softwares com interfaces gráficas como o QGIS. O Python também pode ser usado para fazer processamento em lotes. Em outro post, eu mostrei outra forma de fazer processamento em lote, veja aqui.

Para o processamento de hoje, usaremos o GDAL no Python, instalado através do Conda.

Como fazer isso

A primeira coisa que você precisa fazer é instalar o GDAL no Python usando o Anaconda

conda install gdal

Li que pode ser instalado usando o pip install também, mas não testei:

pip install gdal

e carregue o pacote gdal do osgeo em seu console Python ou script. Eu costumo usar scripts no Anaconda Spyder.

from osgeo import gdal

Sabendo o caminho de seu raster a ser aberto, você pode abri-lo usando o comando a seguir para abrir o raster no modo somente leitura e salvá-lo na variável passo1.

passo1 = gdal.Open ('caminho_do_arquivo.tif', gdal.GA_ReadOnly)

Copie a transformação geográfica para uma variável, que será útil mais tarde.

GT_entrada = passo1.GetGeoTransform()

Leia as bandas de seu raster usando GetRasterBand. No caso, é um raster de somente uma banda.

passo2 = passo1.GetRasterBand(1)

Transforme a banda lida em uma matriz, para poder trabalhar com ela.

img_como_array = passo2.ReadAsArray()

Leia as dimensões do seu array

tam1,tam2 = img_como_array.shape

Crie a saída, a ser salva com as mesmas dimensões do raster de entrada

import numpy as np
saida = np.zeros(shape=(tam1, tam2))

Faça o processamento necessário e salve a saída do mesmo na variável de saída. Por exemplo, todo o raster elevado a 1.2:

for i in range(0,tam1):
    for j in range(0,tam2):
        saida[i,j]=img_como_array[i,j] ** 1.2

A seguir, podemos exportar a variável de saída para um arquivo usando o rasterio.

import rasterio

Pode ser necessário instalar o rasterio digitando no console:

conda install -c conda-forge rasterio

Crie uma variável com a projeção do seu raster. Por exemplo:

dst_crs='EPSG:32722'

Transforme sua saída para float (ponto flutuante)

saida = np.float32(saida)

E, finalmente, salve no arquivo raster de saida

with rasterio.open(
    'raster_saida.tif',
    'w',
    driver='GTiff',
    height=saida.shape[0],
    width=saida.shape[1],
    count=1,
    dtype=np.float32,
    crs=dst_crs,
    transform=GT_entrada,
) as dest_file:
    dest_file.write(saida, 1)
dest_file.close()

Este exemplo salva um arquivo chamado raster_saida.tif

Exemplo com um raster randômico

Eu gerei um raster randômico por distribuição binomial no QGIS 3.18, com valores entre 0 e 10. O intervalo real do raster gerado foi [2,7].

Depois de executar o código do exemplo dado neste post (elevar a 1.2), o seguinte raster foi gerado:

Parece igual, mas a amplitude (min e max) aparece alterada:

Mostrando que o código funcionou!

Edit (2 de agosto de 2021): se você estiver tendo problemas com transformações do GDAL no rasterio, há uma solução rápida para isso. Confira este post.

Edit (21 de março de 2023): escrevi um post explicando como podemos salvar arquivos raster usando o pacote OSGEO (GDAL) no Python. Pode ser útil para alguém que queira usar exclusivamente o GDAL para abrir e salvar arquivos raster usando Python.