• João Ataide

Extração de dados em imagens com Arcpy

Atualizado: 2 de ago.


 

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 pd


També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)
composited




Agora 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)
circleExtract



Outro 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)
maskExtract





Seguindo 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()

An overview of the Extraction toolset—ArcGIS Pro | Documentation