Python e matrizes

Python é a linguagem de programação do momento, especialmente para quem trabalha com geoprocessamento. O QGIS inclusive possui um terminal de python, onde você pode escrever scripts de python e ver o QGIS fazer muito mais que qualquer processo em lote ou modelo pode fazer. Honestamente, é lindo.

Se os problemas que você enfrenta são problemas envolvendo camadas raster, naturalmente estamos falando de matrizes. Afinal, uma camada raster com apenas uma banda consiste em uma matriz bidimensional (uma tabela) de números que são georreferenciados. Mais bandas só adicionam mais matrizes. O que o QGIS faz é renderizar esses dados na forma de uma imagem georreferenciada. Ou seja, uma camada raster possui dois componentes: os dados (a matriz) e os metadados (que informam a geolocalização da imagem).

A matriz de dados é o que importa para quem precisa processar uma camada raster por algum script python. Os metadados podem ficar de lado, pelo menos temporariamente. Depois de processada a matriz, aí sim, os metadados são necessários para reconstruir uma camada raster, com a mesma geolocalização da camada usada como insumo do processo.

Se estamos falando de python e matrizes, obviamente estamos falando da biblioteca numpy . O principal objeto que essa biblioteca oferece são os ndarrays, ou matrizes N-dimensionais. São objetos programados em um nível bem inferior ao da linguagem, o que aumenta muito a velocidade de processamento. Assim, se você quer se aventurar no processamento pytônico de camadas raster, seu objetivo número 1 é obter o ndarray da camada. O seu objetivo número 2 é usar os metadados da camada de entrada para (re)construir uma nova camada raster.

O caminho burocrático

O geoprocessamento surgiu antes de python virar moda. As bibliotecas de código aberto para camadas raster e vetoriais que foram as pioneiras, gdal e ogr, foram escritas em C. Não me pergunte como, mas a comunidade desenvolvedora conseguiu criar uma forma de converter essas bibliotecas para python.

Para usar essas bibliotecas dentro do QGIS, é muito fácil: tudo que você precisa fazer é abrir o editor de python dentro do próprio QGIS e programar.

Voltando ao nosso problema de de processar uma camada raster, abaixo eu apresento um código para calcular a declividade de um modelo digital de elevação. O código pode ser usado como modelo. Os passos do algoritmo são:

1. Obter os nomes das camadas de entrada (existente) e de saída (ainda não existente)
2. Obter os dados da camada de entrada (array).
3. Obter os metadados da camada de entrada.
4. Criar uma camada de saida com os metadados da camada de entrada;
5. Processar os dados de entrada;
6. Escrever os dados de saída na camada de saída;
7. Exportar a camada de saída.

import gdal 
import numpy as np 

# 1) informar os nomes dos arquivos de entrada e saida: 
fnm_input = 'C:/dados/dem.tif' 
fnm_output = 'C:/dados/slope.tif' 

# 2) Obter os dados da camada de entrada 
# abrir arquivo de entrada em modo leitura (parametro 0): 
rst_input = gdal.Open(fnm_input, 0) 

# carregar a banda 1 da camada de entrada para uma matriz numpy: 
dem = rst_input.GetRasterBand(1).ReadAsArray() 

# 3 e 4) Obter os metadados da camada de entrada e criar camada de saida 
# carregar o driver para GeoTiff 
driver_tiff = gdal.GetDriverByName('GTiff') 
# criar um raster de saida com os metadados da camada de entrada 
x = rst_input.RasterXSize y = rst_input.RasterYSize 
# camada raster de saida 
rst_out = driver_tiff.Create(fnm_output, xsize=x, ysize=y, bands=1, eType=gdal.GDT_Float32) 
# configurar a geotransformacao com base na camada de entrada 
rst_out.SetGeoTransform(rst_input.GetGeoTransform()) 
# obter o tamanho da celula 
cellsize = np.abs(rst_input.GetGeoTransform()[1]) 
# configurar a projecao com base na camada de entrada 
rst_out.SetProjection(rst_input.GetProjection()) 
# ler a banda 1 da camada de saida 
bnd1 = rst_out.GetRasterBand(1).ReadAsArray() 

# 5) PROCESSAMENTO (Calcula a declividade): 
grad = np.gradient(dem) 
gradx = grad[0] / cellsize 
grady = grad[1] / cellsize 
gradv = np.sqrt((gradx * gradx) + (grady * grady)) 
slope = np.arctan(gradv) * 360 / (2 * np.pi) 
bnd1 = slope 

# 6 e 7) EXPORTAR CAMADA DE SAIDA 
# sobrescrever a banda 1 do raster de saida 
rst_out.GetRasterBand(1).WriteArray(bnd1) 
# fechar a camada 
rst_out = None 
# carregar no QGIS 
iface.addRasterLayer(fnm_output)

Bom, até aí beleza, pois estamos programando dentro do QGIS. Mas e se você não quer programar dentro do QGIS? Suponha que você use uma IDE de programação, como o Sublime ou PyCharm, e quer trabalhar no seu próprio ambiente?

Nesse caso, é preciso instalar as bibliotecas de geoprocessamento. Para camadas vetoriais, temos a geopandas. Para camadas matriciais, rasterio, pyqgis e a própria gdal.

O problema disso tudo é que realmente não é muito simples fazer essas instalações. Para instalar a gdal fora do QGIS, por exemplo, parece ser tão burocrático que eu simplesmente desisiti. Se você gosta de manter o menor número possível de dependências (e você pode ter motivos muito bons para isso) esse é um caminho não muito favorável.

No meu caso, por exemplo, eu penso em desenvolver aplicações ou algo próximo a aplicações que sejam práticas para outras pessoas. Por isso, valorizo um grau alto de minimalismo nas dependências de bibliotecas.

O caminho da liberdade

Um caminho extremamente prático para processar camadas raster fora do QGIS que eu andei explorando consiste em converter as camadas raster para arquivos .asc. Ao que consta, o formato .asc que estamos falando foi desenvolvido pela Autodesk. Esse formato codifica uma camada raster de apenas uma banda de forma extremamente simples. O arquivo é tão simples que é possível abrir ele em um editor de texto e literalmente editar manualmente os valores dos dados e metadados.

Por ser um arquivo simples, é simples decodificar e codificar ele por scripts de python. E isso nos abre as portas para a liberdade, pelo menos para processar camadas raster fora do QGIS. Tudo que você precisa é numpy e dois códigos de codificação e decodificação, que eu apresento aqui. Nenhuma biblioteca de geoprocessamento adicional. Claro, se você quer visualizar as imagens que está processando, convém também instalar matplotlib.

Passo 1: converter geotiff para asc

Isso é facil. No QGIS, vá em GDAL/OGR >> [GDAL] Conversion >> Translate (convert format). Mantenha todos os parâmetros como padrão.

Ao definir o nome do arquivo de saída, você talvez se surpreenda com a diversidade de formatos que o QGIS consegue (ou tenta) fazer as conversões.

Passo 2: obter os dados e metadados via Python

Para fazer isso, é preciso decodificar o arquivo .asc. Eu programei o seguinte script para isso, que até agora está funcionando:

def input_asc(file):
    """
    Uma função para importar arquivo .asc
    :param file: string do caminho do arquivo, com a extensão '.asc'
    :return: metadados em um dicionário e os dados em um array 2d
    """
    def_f = open(file)
    def_lst = def_f.readlines()
    def_f.close()    
    #
    # get metadata constructor loop
    meta_lbls = ('ncols', 'nrows', 'xllcorner', 'yllcorner', 'cellsize', 'NODATA_value')
    meta_format = ('int', 'int', 'float', 'float', 'float', 'float')
    meta_dct = dict()
    for i in range(6):
        lcl_lst = def_lst[i].split(' ')
        lcl_meta_str = lcl_lst[len(lcl_lst) - 1].split('\n')[0]
        if meta_format[i] == 'int':
            meta_dct[meta_lbls[i]] = int(lcl_meta_str)
        else:
            meta_dct[meta_lbls[i]] = float(lcl_meta_str)
    # array constructor loop:
    array_lst = list()
    for i in range(6, len(def_lst)):
        lcl_lst = def_lst[i].split(' ')[1:]
        lcl_lst[len(lcl_lst) - 1] = lcl_lst[len(lcl_lst) - 1].split('\n')[0]
        array_lst.append(lcl_lst)
    def_array = np.array(array_lst, dtype='float64')
    # replace NoData value by np.nan
    ndv = float(meta_dct['NODATA_value'])
    for i in range(len(def_array)):
        lcl_row_sum = np.sum((def_array[i] == ndv) * 1)
        if lcl_row_sum > 0:
            for j in range(len(def_array[i])):
                if def_array[i][j] == ndv:
                    def_array[i][j] = np.nan
    return meta_dct, def_array

A função definida lê o arquivo, decodifica e retorna duas coisas: os metadados em um dicionário e os dados em uma matriz (array) numpy .

Passo 3: faça o que você quiser com os dados

Por exemplo, para calcular as declividades de um modelo digital de elevação (MDE), eu criei o seguinte script:

def slope(array, cellsize, degree=True):
    """
    Algoritmo para calcular as declividades de um MDE com base na função gradiente da biblioteca Numpy
    :param array: matriz 2d do MDE
    :param cellsize: float valor do tamanho da célula
    :param degree: boolean para definir as unidades de saída. True= graus. False= radianos
    :return: matriz 2d das declividades
    """
    grad = np.gradient(array)
    gradx = grad[0] / cellsize
    grady = grad[1] / cellsize
    gradv = np.sqrt((gradx * gradx) + (grady * grady))
    slope_array = np.arctan(gradv)
    if degree:
        slope_array = slope_array * 360 / (2 * np.pi)
    return slope_array

Com a biblioteca matplotlib podemos visualizar o resultado na IDE (no meu caso, PyCharm)

Passo 4: exportar de volta para um arquivo .asc

Agora precisamos fazer o caminho de volta: codificar dados e metadados em um arquivo .asc. O código, que está funcionando até agora, está abaixo:

def output_asc(array, meta, folder, filename):
    """
    Função para exportar arquivo .asc
    :param array: Matriz 2d com os dados
    :param meta: dicionário dos metadados. Exemplo:

    {'ncols': 366,
     'nrows': 434,
      'xllcorner': 559493.087150689564,
       'yllcorner': 6704832.279550871812,
        'cellsize': 28.854232957826,
        'NODATA_value': -9999}

    :param folder: string da pasta de destino
    :param filename: string do nome do arquivo
    :return: string completa do caminho do arquivo
    """
    meta_lbls = ('ncols', 'nrows', 'xllcorner', 'yllcorner', 'cellsize', 'NODATA_value')
    ndv = float(meta['NODATA_value'])
    exp_lst = list()
    for i in range(len(meta_lbls)):
        line = '{}    {}\n'.format(meta_lbls[i], meta[meta_lbls[i]])
        exp_lst.append(line)
    # data constructor loop:
    for i in range(len(array)):
        # replace np.nan to no data values
        lcl_row_sum = np.sum((np.isnan(array[i])) * 1)
        if lcl_row_sum > 0:
            for j in range(len(array[i])):
                if np.isnan(array[i][j]):
                    array[i][j] = int(ndv)
        str_join = ' ' + ' '.join(np.array(array[i], dtype='str')) + '\n'
        exp_lst.append(str_join)
    flenm = folder + '/' + filename + '.asc'
    fle = open(flenm, 'w+')
    fle.writelines(exp_lst)
    fle.close()
    return flenm

Passo 5: abrir no QGIS o arquivo .asc

Você pode abrir o arquivo diretamente no QGIS. Para retornar ao velho geotiff, é só converter, usando o Translate. Existe um bom motivo para voltar ao geotiff: o arquivo asc é tão simples que não sabe qual é o sistema de referência das coordenadas. Isso você precisa saber previamente. Portanto, recomendo aqui manter esses arquivos como temporários e que não sejam distribuídos para não causar confusões por usuários que não tem como saber qual é o sistema de referência.