Agora vamos falar sobre extração de dados de imagem, etapa do processo de PDI que poderá ser a principal, em especial quando a gente está querendo segmentar uma imagem e classificar os pixels dessas. Basicamente todo modelo que envolvem Machine Learning em imagem terá uma etapa parecida com essa, onde será necessário realizar a extração dos valores dos pixels de uma determinada classe apresentada.
Para começar os nossos testes vamos primeiro importar as bibliotecas principais, mas, dessa vez vamos também importar a biblioteca pandas, que vai ajudar a gente a visualizar os dados quando convertermos em data frame.
# importar bibliotecas
import arcpy
from arcpy import env
from arcpy.sa import *
import pandas as pdTambém iremos definir todos os diretórios, coordenadas da área de interesse e tudo mais que será preciso de configuração do projeto.
workspace = r'D:\\Projeto\\extract.gdb'
basepath = r'D:\\Projeto\\dados'
imagespath = f'{basepath}\\imagens' #Diretório das imagens
cbers = f'{imagespath}\\CBERS' #Diretório do CBERS
nimagem = 'CBERS_4A_WPM_20210609_202_136_L4_BAND'
composited_img = f"{workspace}\\CBERS_4A_WPM_RGB"
roi = f'{basepath}\\ROI.shp'
coords = "467321.538879395 7933874.64569092 481460.825317383 7949490.0112915"Um procedimento interessante a se fazer, é criar um laço de repetição para gerar uma lista com o diretório de cada banda da cena, então ficará mais fácil selecionar a banda pretendida somente usando o índice dessa.
bandas = []
for i in list(range(0, 5)):
#Diretório banda original
banda = f'{cbers}\\{nimagem}{i}.tif'
#Diretório nova imagem recortada
cliped = f'{workspace}\\rgb_cliped_BAND{i}'
#Recort da imagem
arcpy.Clip_management(banda, coords, cliped, roi)
bandas.append(f'{workspace}\\rgb_cliped_BAND{i}')Então para compor uma imagem, é somente criar uma lista com as bandas na sequência desejada e realizar o Composite Bands.
#Listando composição RGB
list_rgb = [bandas[3], bandas[2], bandas[1]]
arcpy.management.CompositeBands(list_rgb, composited_img)
composited = Raster(composited_img)
compositedAgora vamos fazer outro procedimento, neste caso uma aritmética de bandas criando um índice espectral entre as bandas red e nir, chamado de NDVI (Índice de Vegetação Normalizada).
#Crindo NDVI
red = Raster(bandas[3])
nir = Raster(bandas[4])
NDVI = ((nir - red)/ (nir + red))Com isso, podemos usar a primeira função do módulo extract, que é o Extract By Atributes, esse você pode informar um intervalo de valores, como nesse exemplo estou usando o valor de intervalo entre 0.33 - 0.66 que indica que a vegetação está moderadamente saudável.
# Extrair pelos valores da imagem
calculate = "VALUE >= 0.33 AND VALUE <= 0.66"
attExtract = ExtractByAttributes(NDVI, calculate)
attExtract
Outras funções que podemos fazer nesse caso é criar máscaras, mas imagens, sem precisar realizar um recorte por exemplo. Podemos extrair uma parte de uma imagem, a partir de um raio e um ponto, usando a função Extract By Circle, pode ser bastante útil quando queremos, por exemplo, extrair os pivôs centrais de uma fazenda, por exemplo.
# Extrair por um circulo
calculate = "474506.74 7942693.72"
raio = 4800
circleExtract= ExtractByCircle(composited, calculate, raio)
circleExtractOutro que é importante, é a extração por máscara, para o pessoal que trabalha com visão computacional e usa o OpenCV está acostumado com isso, e o Arcpy não é diferente, basicamente precisamos passar uma máscara que pode ser um shapefile, por exemplo, e a imagem e a função vai retirar a parte deseja. Por padrão a função usa o INSIDE, mas você pode passar OUTSIDE e destacar somente o que está fora a imagem.
#Extratir por uma máscara
mask = f'{workspace}\\buffer_polygon'
maskExtract = ExtractByMask(composited, mask)
maskExtractSeguindo esses passos, podemos também extrair parte de uma imagem a partir de polígonos gerados por pontos, mas nesse caso não ficou um resultado interessante, pois os pontos não seguiam uma origem.
# Extrair por um polígono
extPolygonOut = ExtractByPolygon(composited, pontos)
extPolygonOut
Você também pode extrair a partir de um retângulo, passando somente os pontos dos vértices e se ele será OUTSIDE ou INSIDE, como na máscara.
# Extrair por um retângulo
rectExtract = ExtractByRectangle(composited, Extent(469064, 7949080, 479480, 7935744),"OUTSIDE")
rectExtract
Agora vamos entrar de fato nas funções do módulo que trabalham com pontos, esse pode ser o conjunto de funções mais importantes para extrair valores de uma imagem. O primeiro é o Extract by Points. Nessa função eu tinha pontos selecionados nas áreas de coletas, realizei então o addXY para automaticamente incorporar os valores das coordenadas na tabela de atributos, e então criei um laço de repetição transformando na classe ponto.
# Definir diretório de pontos
points = f'{workspace}\\pontos_amostra
# Adicionar coordenadas
arcpy.AddXY_management(points)
df = pd.DataFrame.spatial.from_featureclass(points)
#Listar pontos
pontos = []
for x in df.POINT_X:
for y in df.POINT_Y:
pontos.append(arcpy.Point(x, y))Então realizei a extração dos valores dos pontos somente da banda verde. Ponto importante é que essa função gera um arquivo raster, então tive que realizar a conversão para polígono realizando a leitura desse podemos ver então o valor dos pixels interceptados pelos pontos.
pointEx = f'{basepath}\\pontos_ex.shp'
pointExtract = ExtractByPoints(bandas[2], pontos)
#Convertendo por ponto
arcpy.RasterToPolygon_conversion(pointExtract, pointEx)
#Ler como DataFrame
df2 = pd.DataFrame.spatial.from_featureclass(pointEx)
df2.head()
Extract values to Points segue a mesma lógica, necessita passar um ponto, a imagem e output, nesse caso ele gera um arquivo vetorial, outro ponto é que podemos também escolher quais valores que iremos extrair.
out_values = f'{workspace}\\pontos_values'
# Extrair Valor por pontos
ExtractValuesToPoints(points , composited, out_values, "INTERPOLATE", "VALUE_ONLY")
#Importar como DataFrame
df3 = pd.DataFrame.spatial.from_featureclass(out_values)
df3.head()
Por último e mais importante é o sample, esse gera diretamente o arquivo de forma tabular e se a gente passar uma composição já separa o valor de cada um em colunas específicas. Todas essas funções, quando não encontram um valor no pixel ou quando o ponto está entro dois ou mais pixel, ele interpola e gera um valor de interpolação, por isso é de extrema importância ter cuidado na seleção dos pontos.
ob_points = f'{workspace}\\observations' # Pode ser um polígono
# Amostrar por ponto
Sample([composited], points, ob_points)
#Ler como data frame
df4 = pd.DataFrame.spatial.from_table(ob_points)
df4.head()