Janelas móveis auto-reguladas em camadas matriciais

Um problema típico no processamento de camadas matriciais é a análise de vizinhança das células (píxel) vizinhas. Inclusive células não tão vizinhas assim, mas que localizem-se dentro de um raio de distância da célula de referência, como no caso do TPI, que analisa o desvio do valor de cota local em relação à média de uma vizinhança “alargada”.

Quando estamos programando em python, essa análise de vizinhança é feita pela técnica da “janela móvel” – que no caso é bidimensional no mapa. Assim, a matriz do mapa é varrida célula por célula com dois laços de iteração aninhados. Mas antes precisamos carregar uma biblioteca muito útil para matrizes, numpy:

# [1] -- importar bibliotecas
import numpy as np

Como exemplo, vou aqui usar uma matriz aleatória de 100 x 100 células. Na prática teremos que carregar a matriz a partir de uma camada raster usando GDAL ou Rasterio.

# [2] -- carregar matriz do mapa
n_size = 100
grd = np.random.random(size=(n_size, n_size))
# computar atributos
n_rows = len(grd)
n_cols = len(grd[0])

Uma janela de inspeção é facilmente criada pela técnica de fatiamento (slicing) da matriz numpy. Por exemplo, grd[10:21, 15:26] retorna uma fatia de 10×10 da matriz original grd, entre as linhas 10 e 20 e 15 e 25 (o último número é excludente).

# [3] -- exemplo de fatiamento
grd_view = grd[10:21, 15:26]
print(grd_view)
print(np.shape(grd_view))

Na ideia de uma janela móvel de raio r ao redor de uma célula de linha i e coluna j Isso equivale à seguinte sintaxe: grd[i - r:i + r + 1, i - r:i + r + 1] Mas se você deseja uma janela com um raio pré-definido, isso só funciona no miolo do mapa, longe das bordas em que o raio extrapola os limites do próprio mapa. Eventualmente algum dos lados da janela (inferior, superior, lateral esquerdo e lateral direito) irá exceder o limite da matriz, o que certamente retornará um erro.

Aqui, eu bolei uma solução para auto-regular os diferentes raios de busca (inferior, superior, lateral esquerdo e lateral direito). Para implementar a solução, é preciso declarar duas funções:

# [4] -- declarar funcoes auxiliares para auto-regular o raio de inspecao

def lower_radius(n_position, n_radius_max=3):
    """
    :param n_position: int cell position
    :param n_radius_max: int window radius
    :return: lower radius
    """
    if n_position < n_radius_max:
        return n_position
    else:
        return n_radius_max

def upper_radius(n_position, n_max, n_radius_max=3):
    """
    :param n_position: int cell position
    :param n_max: int cell upper bound   
    :param n_radius_max: int window radius
    :return: lower radius
    """
    if n_position <= n_max - 1 - n_radius_max:
        return n_radius_max
    else:
        return n_max - 1 - n_position

Essas funções são genéricas para linhas e para colunas, retornando o raio inferior e superior (ou lateral). Agora é preciso definir o raio da janela:

# [5] -- definir o raio (maximo) da janela movel 
n_window_radius = 2

E executar finalmente o laço duplo para varrer a matriz da camada raster:

# [6] -- varrer a matriz
for i in range(n_rows):
    # obter os limites de linha
    n_row_low = i - lower_radius(n_position=i, n_radius_max=n_window_radius)
    n_row_upp = i + 1 + upper_radius(n_position=i, n_max=n_rows, n_radius_max=n_window_radius)
    for j in range(n_cols):
        # obter os limites de coluna
        n_col_low = j - lower_radius(n_position=j, n_radius_max=n_window_radius)
        n_col_upp = j + 1 + upper_radius(n_position=j, n_max=n_cols, n_radius_max=n_window_radius)
        # janela movel auto-regulada:
        grd_window = grd[n_row_low:n_row_upp, n_col_low:n_col_upp]

A variável grd_window aqui recebe a matriz da janela para cada célula de análise. Aqui é o momento para customizar o código e executar as operações e processos de interesse.