Carlos Silva

23 de fev. de 20226 min de leitura

Detecção de anomalias - Identificando áreas florestais colhidas com PyCaret e a API Python do GEE

Atualizado: 16 de mar. de 2022

Introdução

Operações de colheita florestal podem ser bem complexas, a operação fica isolada e nem sempre os dados são transportados da melhor maneira, as tecnologias embarcadas apesar se promissoras tem conexão com a internet limitadas no ambiente rural e uma única companhia pode ter várias equipes de campo diferentes, sendo próprias ou terceirizadas. Todo esse cenário torna a vida dos responsáveis pelas informações florestais um terror, gerenciar caixa de e-mail, ligações e mensagens de texto para conseguir atualizar a informação milhares de talhões florestais, garantir uma cadeia de produção leve e ágil pode ser um tremendo desafio nesse caso. Por isso é comum explorar tecnologias que estabeleçam fluxos de dados concisos, atualizados e padronizados, que garantam uma cadeia de produção robusta.

Corroborando com essa dificuldade, minha proposta hoje é prototipar um roadmap em Python que explore dados orbitais armazenados em nuvem com a API Python do Google Earth Engine (GEE) e criar um modelo de identificação de anomalias, que nesse caso servirá para identificar áreas que foram colhidas.

Preparando o ambiente

Essas foram as bibliotecas usadas na construção desse protótipo. Em especial a Geemap e a PyCaret foram fundamentais, recomendo dar uma olhada nelas. O Geemap entrega uma ambiente interativo e diverso melhorando e otimizando os recursos disponíveis da API Python oficial do GEE e o PyCaret tem uma proposta de ser uma biblioteca "Low-Code" para machine learning.

#Para para carregamento e geoprocessamento
import ee
from ipygee import*
import geemap
import geopandas as gp

#Para Análise e Exploração dos Dados
import os
import datetime
import pandas as pd
import matplotlib.pyplot as plt

#Para modelagem
from pycaret.anomaly import setup , create_model , assign_model , plot_model , predict_model , save_model , load_model

Para usar os dados do GEE é preciso estabelecer uma conexão com a API da seguinte forma.

ee.Authenticate()
ee.Initialize()

Será gerado um link de autenticação de usuário Google e a criação de um Token de acesso.

Dados do Projeto

Os dados desse projeto são um conjunto centroides de áreas de floresta plantada da região da cidade de Três Lagoas do estado do Mato Grosso do Sul, Brasil. Processei eles me baseando em dados do projeto MapBiomas Coleção 5 de 2019. O link para o código completo está aqui , nesse post, vou pontuar principais partes do código.

#Pasta do projeto
FolderProj = '/Anomaly Detection - Identifying harvested areas with sentinel imagery with the GEE/'

PointFiler = FolderProj+'/pointsPlantedForest.shp'

pointsPlantedForest = gp.read_file(PointFiler)

#Faço uma cópia do arquivo de pontos. Dissolvo (ou resumo) os vários pontos em um só.
FileExtend = pointsPlantedForest.copy().dissolve() 

#Crio um polígono do tamanho das proporções máximas do polígono dissolvido anteriormente aumento em 5 km de borda. Essa será  área de estudo
FileExtend['geometry'] = FileExtend['geometry'].envelope.buffer(5000)

#Mostro a tabela de pontos
pointsPlantedForest

Agora é só carregar pra dentro do GEE com a ajuda do geemap.

Points = geemap.shp_to_ee(PointFiler) 
AOI = geemap.geopandas_to_ee(FileExtend)

O que é NDVI?

O índice de vegetação de diferença normalizada (NDVI) é um indicador simples que pode avaliar se uma região tem vegetação verde viva através de sensores remotos de observação da terra em satélites. Objetivamente valores de 0 a 1 indicam alguma atividade vegetativa, o comportamento desse índice é fortemente influenciado pelo clima das estações do ano de como cada vegetação reage a isso. Veremos isso na prática no decorrer desse post.

Extraindo série de dados NDVI

Com a conexão ao GEE ativa podemos acessar seu amplo catalogo de sensores remotos, para esse protótipo escolhi usar o catalogo de imagens NDVI MODIS cujo os detalhes você pode encontrar aqui. O processo de carregar essa série de imagens foi conforme o seguinte código.

#Função para transformar os valores dos pixeis para o intervalo de -1 a 1 em tipo float e insere a propriedade da data das imagens
def NDVIf(image):
  return image.expression('float(b("NDVI")/10000)').copyProperties(image, ["system:time_start"])

#Definindo o intervalo de data da serie de imagens.
startDate = ee.Date.fromYMD(2009,1,1);
endDate = ee.Date.fromYMD(2019,12,31);

#Carrego o dataset de imagens NDVI MODIS
modis = ee.ImageCollection("MODIS/006/MOD13Q1")
#Aplico os filtros na coleção de imagens NDVI MODIS, de área, datas e a função de representação de valores do NDVI
NDVI_Modis = modis.filterBounds(AOI).filterDate(startDate,endDate).map(NDVIf)

Através de uma função, abaixo obtemos conseguimos agregar os dados em uma série temporal de NDVI e armazenar em um conjunto de dados do pandas (pandas.DataFrame). Essa função extraí a média dos valores NDVI das imagens disponíveis para o conjunto de pontos de floresta plantada nas datas de obtenção de cada imagem ou conjunto de imagens.

#Função do ipygee para obtenção de gráficos derivados de coleção de imagens do GEE
point_ndvi_Modis = chart.Image.series(**{'imageCollection': NDVI_Modis,
                                   'region': Points,
                                   'reducer': ee.Reducer.mean(),
                                   'bands' : 'NDVI',
                                   'scale': 25,
                                   'xProperty': 'system:time_start'})
point_ndvi_Modis.renderWidget(width='40%')
#Crio um dataframe pandas com os dados da serie temporal
p1_dataframeModis=point_ndvi_Modis.dataframe
p1_dataframeModis.describe()

Análise e Exploração dos Dados

Foi selecionado um intervalo entre 2009 e 2019 para comportar nossos dados iniciais, observamos abaixo que há um certo padrão sazonal, que é esperado para uma série temporal de NDVI, que como disse é muito influenciada pelas condições climáticas das estações do ano. Além disso é observado duas ocorrências negativas bem marcadas.

Analisando esse dado em um gráfico box-plot fica destacado dois pontos bem distantes dos limites inferiores, que junto com a lógica do dado não faz sentido além de serem Outliers já que não faz sentido uma floresta decair para uma valor de NDVI tão baixo e se recuperar tão rápido logo em seguida, considerando que o índice reflete o comportamento natural de uma vegetação. Isso é bom, porque essa queda rápida é justamente o que queremos identificar.

Modelando no PyCaret

Como comentado anteriormente uma série temporal NDVI é fortemente influenciada pelas condições climáticas sazonais advindas das estações do ano, neste caso a sazonalidade é sendo representada como tendo a baixas no meio do segundo semestre dos anos e a alta no meio do primeiro semestre. Então decidi adicionar as colunas de 'dia do ano' e 'semana no ano' aos dados para que esse tipo de comportamento seja aprendido com mais facilidade no nosso modelo de anomalia. Também separei os dados em treino e validação na proporção de 95% e 5% respectivamente.

dataset = p1_dataframeModis
#Criando novas colunas
dataset['dia_do_ano'] = [i.dayofyear for i in dataset.index]
dataset['semana_do_ano'] = [i.weekofyear for i in dataset.index]
#Separando treino e teste
data = dataset.sample(frac=0.95, random_state=786)
data_unseen = dataset.drop(data.index)
data

No PyCaret existe um processo de analise e limpeza de dados que ocorre na utilização da função setup, ele é necessário antes da produção de qualquer modelo no PyCaret, essa função cuida de uma série de questões como valores vazios, categorização dos tipos de coluna, normalização dos dados e uma série de outros tratamentos necessários para cada tipo de modelo. Ao executar essa função o PyCaret também elabora uma espécie de perfil dos dados que estão sendo carregados.

exp_ano101 = setup(data, normalize = True, session_id = 123)

Após a inicialização dos dados com a função setup, é só chamar a função para criar o modelo(create_model) e comunicar qual modelo se quer criar, nesse protótipo o modelo escolhido foi o 'iforest' ou 'Isolation Forest' que tem um comportamento análogo a de uma árvore de decisão, porém foca em encontrar os valores que são facilmente isolados ao se recortar nossa base de dados . O parâmetro fration definido representa a porcentagem de outliers esperados, no nosso caso baseado nos nossos dados de treino temos aproximadamente 1% (2/240).

#Cria o modelo.
iforest = create_model('iforest', fraction = 0.01)
#Atribui rótulo de dados aos dados de treino do modelo.
iforest_results = assign_model(iforest)
#Mostrando quais itens dentro dos dados de treino foram considerados anomalias.
iforest_results[iforest_results['Anomaly'] == 1].head()

#Mostrando o modelo graficamente
plot_model(iforest,feature="NDVI")

data_unseen['NDVI'][0] = 0.5 #Inseri um dado anômalo na base de teste.
unseen_predictions = predict_model(iforest, data=data_unseen)
unseen_predictions[unseen_predictions['Anomaly'] == 1].head()

Após isso o modelo já está pronto, e pode ser salvo para ser usado em processos.

#Salvando o modelo.
save_model(iforest,'Anomaly_Detection_NDVI_IForest_Model')

Pondo em prática

Agora temos em mãos a seguinte série de dados NDVI MODIS de um ponto de floresta plantada da mesma região conforme abaixo. Considere os mesmos métodos de obtenção de dados só que agora para um ponto em especifico na Cidade de Três Lagoas-MS, Brasil (-51.8186,-20.4644).

dataset = p1_dataframeModis2
#Adicionando as outras colunas auxiliares
dataset['dia_do_ano'] = [i.dayofyear for i in dataset.index]
dataset['semana_do_ano'] = [i.weekofyear for i in dataset.index]
new_prediction = predict_model(Anomaly_Detection_Model, data=dataset)
#Excluindo anomalias positivas, já que não é nosso foco
anomaly = new_prediction[new_prediction['Anomaly'] == 1][new_prediction['NDVI']<0.6]

Analisando mais de perto as últimas anomalias, elas iniciaram no fim de 2017, isso pode ser confirmado analisando a série de imagens do ponto que usamos de teste.

Conclusão

É valido que o protótipo pode ser funcional, as técnicas aplicadas aqui em conjunto com base de dados consolidadas de um cadastro florestal e imagens de maior resolução poderiam ser de grande valia na consolidação de cadeia de dados de produção florestal robusta. Para continuar esse trabalho novas rodadas de desenvolvimento seriam necessárias a partir daqui, principalmente focando nos seguintes itens.

Busca de imagens de sensores remotos de maior resolução para monitorar áreas menores que 250 m². As amostras pontuais desse protótipo eram de áreas maiores que essa restrição implicadas pelas imagens MODIS.
Usar uma máscara de nuvens para ignorar ou classificar anomalias causadas pela a ocorrência das mesmas.
Aumentar a quantidade de amostras, diversificar e acurar os dados. Como o dados base desse protótipo foram baseados nos dados do MapBiomas, eu não pude contar com informações como datas de colheita e plantio, espécie e outra série de dados que poderiam auxiliar na construção de um modelo melhor.