Processando camadas raster com rasterio e colab notebooks

Uma alternativa relativamente interessante para processar camadas raster em Python é com a biblioteca rasterio. Essa biblioteca é construída em cima da boa e velha GDAL, mas a vantagem é que muito mais “pytônica” (já que a GDAL é escrita em C originalmente). Outra vantagem (e na minha opinião a melhor vantagem) é que é possível usar essa biblioteca diretamente na plataforma de programação em nuvem google colab. Para isso execute a primeira célula com o comando de instalação:

# [1] -- install rasterio
!pip install rasterio

Obviamente, importe a biblioteca (também aproveite para importar outras bibliotecas extras):

# [2]
import rasterio
# some extra libs:
import numpy as np

Na próxima célula informe o caminho do seu arquivo raster de entrada (ou os caminhos, se for mais de uma camada):

# [3] -- define input file path
s_filepath_input = "/content/input_raster_layer.tif" 

A receita para carregar uma camada raster é a seguinte:

# [4] -- load raster layer
rst_input_layer = rasterio.open(s_filepath_input)

E para se acessar a matriz numpy na primeira banda da camada:

# [5] -- access numpy 2d array
grd_input_layer = rst_input_layer.read(1)

Agora sim: podemos nos divertir processamento essa matriz. Um exemplo típico seria se obter uma cena booleana da camada:

# [6] -- process numpy array:
grd_output_layer = 1.0 * (grd_input_layer == 10)  # boolean map for 10 values

A exportação é uma operação um pouco mais complexa. Vamos considerar o caso simples em que a camada de saída é uma réplica (em tamanho e extensão da camada de entrada. Primeiro informe o caminho do arquivo de saída:

# [7] -- define output file path
s_filepath_output = "/content/output_raster_layer.tif" 

A seguir aplique a receita:

# [8] -- export layer
# open file
with rasterio.open(
s_filepath_output,
'w', # write mode
driver='GTiff', # driver
height=grd_output_layer.shape[0], # sizes
width=grd_output_layer.shape[1], # sizes
count=1,
dtype=np.int8, # change for different output data types (eg, np.float32, np.uint16)
crs='EPSG:4326',
transform=rst_input_layer.transform, # take the same transform from input
) as rst_output_layer:
    # write file
    rst_output_layer.write(grd_output_layer, 1)
# close file
rst_output_layer.close() 

E está feito.