Como usar o GDAL para salvar arquivos raster no Python

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

Bem-vindos ao meu blog!

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

Hoje, eu escrevo sobre como salvar rasters usando GDAL em Python, sem usar o pacote rasterio.

Eu tenho um post anterior sobre leitura e salvamento de dados raster em Python, que pode ser útil para quem está lendo o post atual.

Soluções usando rasterio funcionam, todos nós sabemos disso. No entanto, instalar o rasterio e o GDAL no mesmo ambiente do conda pode ser um desafio… Então, geralmente, eu tento não instalar os dois juntos.

Vou mostrar para vocês a maneira que eu acho que é a mais fácil para salvar dados raster usando o GDAL no Python.

Por aqui, estou usando Python 3.10, GDAL 3.0 e a minha IDE de escolha é o Spyder 5.3.3, executado em uma máquina com Windows 10.

Importar o pacote do GDAL no Python

from osgeo import gdal

Importar a referência espacial OGR

import osr

Escolher salvar o raster como .tif

driver = gdal.GetDriverByName('GTiff')

Definir o local de salvamento do seu raster, com a extensão

path_output = 'C:/path_to_save.tif'

Definir o número de pixels no seu raster

Geralmente, isso é definido pela array que você quer salvar como raster.

size1,size2=img_as_array.shape

Crie o espaço onde seu conjunto de dados será colocado

Atenção: as opções do seu driver podem ser diferentes dessas.

Especificamente, fique atento ao tipo do dado que você está salvando como raster (o quinto argumento). Verifique todos os tipos disponíveis na documentação GDAL.

O quarto argumento, aqui simplesmente listado como “1”, é o número de bandas que seu raster de saída deve ter.

dataset_output=driver.Create(path_output,size2,size1,1,gdal.GDT_Float64)

Também é possível anexar outras opções em forma de strings, após esses cinco argumentos obrigatórios.

Leia mais sobre GDAL Create na documentação do GDAL.

Defina a transformação geográfica

Normalmente, essa informação será baseada em outro raster.

Use o trecho abaixo se precisar extrair uma GeoTransform de outro raster.

orig_file = gdal.Open(filename, gdal.GA_ReadOnly)
GT_output = orig_file.GetGeoTransform()

Defina a GeoTransform:

dataset_output.SetGeoTransform(GT_output)

Defina a referência espacial do seu GeoTIFF

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')

Obs: Substitua “WGS84” por outro sistema de referência geográfica, quando for o caso.

dataset_output.SetProjection(srs.ExportToWkt())

Finalmente, neste passo salvaremos a nossa array dentro do raster

export= dataset_output.GetRasterBand(1).WriteArray(img_as_array)

Libere a memória que estava sendo utilizada

export = None

Pronto!

Estes passos devem ter resultado em um raster devidamente salvo usando Python.

Extras

E se eu preferir usar rasterio (também conhecido como “rio”)?

Você pode continuar usando seus pacotes preferidos. Este post simplesmente mostra que é possível salvar rasters no Python usando o pacote do GDAL.

Por que devemos definir a variável “export” como None?

Acredito que o GDAL tenha uma maneira explícita de fechar os conjuntos de dados, GDALClose. Definindo os dados exportados com o valor None, sinalizamos que não estamos mais usando esse espaço de memória. Isso também ajuda a evitar erros, como tentar abrir um arquivo no QGIS ou ArcMap e receber uma mensagem de erro informando que o arquivo está sendo usado por outro programa ou que o raster não pôde ser aberto.

Por que devemos pesquisar o driver GTiff pelo nome?

Esta é a maneira mais simples de selecionar o driver certo para salvar seu arquivo com a extensão desejada.

A lista dos drivers disponíveis pode ser encontrada na documentação do GDAL.