Operações e condicionais com Arcpy

---

Irei começar agora uma série de artigos mais voltados para a parte de visão computacional, na prática com imagens de satélite. Para esse primeiro aqui estou usando a biblioteca Arcpy (Bibliotecas do ArcGIS). Para iniciar qualquer projeto de PDI a primeira coisa é definir os nossos diretórios de trabalho, nesse caso usei um geodatabase como meu workspace, coloquei uma pasta base onde colocaria os dados externos e defini dentro dela uma pasta que contém as imagens de satélite.

# Importar bibliotecas
import arcpy
from arcpy.sa import *

# Configurar ambiente
workspace = r'D:\\Projeto\\Spatial Analyst.gdb'
basepath = r'D:\\Projeto\\dados'
imagespath = f'{basepath}\\imagens'

Neste caso vou trabalhar com as imagens do CBERS 4A, baixei esses dados diretamente do site do INPE e descompactei todas as imagens na pasta.

Pegando o nome da nossa imagem, criando uma lista e colocando o número de bandas que o satélite possui, que nesse caso são 5, podemos criar uma lista com o diretório de cada uma das imagens, facilitando a chamada das bandas separadamente.

# Listar bandas na cena
cbers = f'{imagespath}\\CBERS'

bandas = []
nbands = 5
nimagem = 'CBERS_4A_WPM_20210609_202_136_L4_BAND'

for i in list(range(0, nbands)):
    bandas.append(f'{cbers}\\{nimagem}{i}.tif')

Dentro do arcgis podemos trabalhar diretamente com cada uma das bandas se quiséssemos, com o simples comando abaixo, neste caso estou pegando a banda 1 que para esse satélite é a banda da cor azul.

# Importando como Raster
band1 = Raster(bands[1])
band

Compondo Bandas e Clipando imagem

Para compor as bandas, ou seja, colocar cada banda dentro de um canal de cor é simples, necessita somente do comando Composite Bands e passar a ordem exata das bandas, exemplo: queremos uma imagem RGB e temos o CBERS 4A que tem a banda 0 a pancromática e as banda 1, 2, 3, 4 são respectivamente sensores que captam as bandas Blue, Green, Red e Infravermelho.

Sabendo isso é só criar uma lista com o diretório de cada uma das bandas necessárias.

# Diretorio de imagem composta
composited = f"{workspace}\\CBERS_4A_WPM_20210609_RGB"

# Lista de sequência de composição
list_rgb = [bandas[3], bandas[2], bandas[1]]

# Comport imagem
arcpy.management.CompositeBands(list_rgb, composited)

composited = Raster(composited)
composited

Para melhorar o desempenho e realizar as operações aqui de forma mais rápida, vou realizar o clip da imagem por uma área de interesse, passamos então uma coordenada base, a área de interesse e pronto, temos uma imagem recortada.

# Arquivo de área de interesse
roi = f'{basepath}\\ROI.shp'

# Coordenadas padrão
coords = "467321.538879395 7933874.64569092 481460.825317383 7949490.0112915"

# Recortar imagem
rgb_clipped = f'{workspace}\\clipped'
arcpy.Clip_management(composited, coords, clipped, roi)

rgb_clipped = Raster(rgb_clipped)
rgb_clipped

Realizar operações simples fica fácil quando importamos como Raster, você pode fazer qualquer operação como somar um valor constante nas bandas, inverter a escala de cores, negativar o valor da escala de cores ou até mesmo elevar exponencial por um valor.

Dá mesma forma fazer uma operação de normalização para um índice também fica muito fácil, sendo preciso colocar exatamente como a equação pede e pronto, como mostro no exemplo do NDVI.

# NDVI
(band4 - band3) / (band4 + band3)

Condicional

Agora vamos para as funções condicionais, basicamente é a criação de um if/else, você passa a imagem que vai trabalhar, o valor se a condição for verdadeira, valor da falsa e a condição. É uma função muito interessante, pois você pode cria mascaras para um modelo de deep learning, por exemplo, selecionando uma região de cor da imagem como uma classe.

Nesse exemplo passei a imagem RGB recortada, o valor que colocarei na imagem caso a condição seja verdadeira 250 e o valor caso 0, e passando a condição de uma faixa de valores entre 160 e 200 para cada um dos canais.

# Realizar um condicional
cond1 = Con(rgb_clipped, 250, 0, "Value >= 160 and Value <=200")
cond1

Mas você não precisa passar uma faixa de valor, mas um threshold, por exemplo, colocando que todo valor maior que 70 é 0 e menor 255. E com isso basicamente conseguindo delimitar quase toda a área do corpo de água que apresenta na imagem.

# Realizar um condicional
cond2 = Con(rgb_cliped, 0, 255, "Value <= 70")
1/cond2  # Inverso estabelece máscara binária

Set Null

Outra função muito legal é o set null, essa é uma função bem simples que lhe serve setar qual faixa de valor ou um valor da imagem terá o valor nulo, você pode usar, por exemplo, para retirar aqueles pixels outlier de sua imagem que possuem valor, ou até mesmo criar uma masca, como no exemplo abaixo:

Basicamente coloquei que todos pixels com valor maior ou igual a 200, passando um valor constante ou raster contantes para caso não exista tal valor.

# Definir valor nulo
nulo = SetNull(rgb_clipped, 0, "VALUE >= 200")
1 / nulo # Inverso estabelece máscara binária

Pick value

O último é o pick value que é um condicional um pouco mais complexo, ele pega o valor de posição de um raster e passa uma lista de imagens e o tipo de operação se é banda única ou banda múltipla.

Aqui basicamente estamos fazer uma "substituição" entre a máscara que fiz anteriormente com a imagem composta, pegando o valor de uma imagem e colocando na outra.

# Realizar pick entre a máscara e a imagem
pick = Pick(cond2, [rgb_clipped], "MULTI_BAND")
pick

An overview of the Conditional toolset—ArcGIS Pro | Documentation