Laboratório 4#

Mosaicos e classificação supervisionada de imagens no GEE

Objetivos:

Remoção de nuvens
Composições e mosaicos de imagens
Classificação supervisionada no GEE

[ ]:

import ee
import geemap

[ ]:

# Ref: https://developers.google.com/earth-engine/apidocs/ee-authenticate
# Para inicializar a sessão para execução insira o id do projeto em ee.Initialize().
ee.Authenticate()
ee.Initialize(project='id_projeto')

Preparar dados#

Filtrar coleção de imagens (nuvens, datas e região)#

Dessa vez vamos utilizar uma coleção de imagens ao invés de uma imagem única, isso garantirá maior robustez ao nosso modelo de classificação.

Adicionalmente neste laboratório, diferentemente do anterior, utilizar imagens sem nuvens serão essenciais. Para eliminar as nuvens, captura-se múltiplas imagens da mesma área, remove-se as nuvens e, usando a técnica do mosaico, integram-se todas em uma única imagem. Isso preenche lacunas deixadas pelas nuvens removidas, pois os pixels faltantes em uma imagem são substituídos por pixels válidos de outras. O GEE prioriza o pixel da imagem que está no topo da coleção para compor o mosaico.

Vamos ver como fazer isso.

[ ]:

# Essa é uma função padrão do
def mask_s2_clouds(image):
    """Masks clouds in a Sentinel-2 image using the QA band.

    Args:
        image (ee.Image): A Sentinel-2 image.

    Returns:
        ee.Image: A cloud-masked Sentinel-2 image.
    """
    qa = image.select("QA60")

    # Bits 10 and 11 are clouds and cirrus, respectively.
    cloud_bit_mask = 1 << 10
    cirrus_bit_mask = 1 << 11

    # Both flags should be set to zero, indicating clear conditions.
    mask = qa.bitwiseAnd(cloud_bit_mask).eq(0).And(qa.bitwiseAnd(cirrus_bit_mask).eq(0))

    return image.updateMask(mask).divide(10000)

Para nos ajudar com a análise, vamos filtrar somente as imagens que passam por São Paulo. Para isso definimos um ponto e depois utilizamos o método filterBounds() do GEE

[ ]:

sao_paulo = ee.Geometry.Point(-46.711, -23.641)
sao_paulo

Vamos definir os limites para cobertura de nuvem e filtrar a coleção de imagens. Agora vamos filtrar o período em que queremos trabalhar com as imagens.

Alguns comentários:

Selecionar um período muito longo pode resultar em uma coleção muito grande e pesada para trabalhar
Para fins de exemplos vamos cobrir o ano de 2020

💡 A dica aqui é trabalhar com o intervalo de datas e valor de cobertura por nuvens até se obter uma coleção de tamanho suficiente para a composição da imagem de forma apropriada

[ ]:

limite_cobertura_nuvem = 5  # máximo de 10% de cobertura de nuvem
data_inicial, data_final = "2018-01-01", "2022-01-30"

dataset = (
    ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
    .filterBounds(sao_paulo)
    .filterDate(data_inicial, data_final)
    .filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE", limite_cobertura_nuvem))
    .map(mask_s2_clouds)
)
dataset

Podemos verificar o número de imagens resultantes na coleção com o método size():

[ ]:

print(f"total de imagens na coleção: {dataset.size().getInfo()}")

Recortar a região de interesse (Region of Interest - ROI)#

Primeiramente vamos definir um polígino de interesse para nossa análise.

A fins de exemplo, vamos definir programaticamente um polígono retangular que cobre a cidade de São Paulo. Mas tenha em mente que você também pode selecionar no mapa usando a ferramenta de desenho.

[ ]:

# roi: região de interesse
roi = ee.Geometry.BBox(
    west=-46.88655, south=-23.6906, east=-46.611642, north=-23.590958
)
roi

Como no trabalho a ideia é utilizar uma coleção de imagens, o método clip() não pode ser aplicado diretamente.

Para efetuar o procedimento uma função personalizada deve ser criada. No caso, abaixo criamos a função clip_img_collection(), que recebe uma coleção de imagens e uma geometria (no caso, o nosso polígono de interesse) e retorna a coleção já recortada.

[ ]:

def clip_img_collection(dataset, geometria):
    """Recorta todas as imagens de uma coleção de imagens do Earth Engine para
    uma região específica definida por uma geometria.

    Essa função cria internamente uma função de recorte específica para a
    geometria fornecida e aplica essa função a cada imagem na coleção de imagens
    usando o método map().

    Parameters
    ----------
    dataset : ee.ImageCollection
        A coleção de imagens a ser recortada. Cada imagem da coleção será
        recortada para se adequar à geometria fornecida.

    geometria : ee.Geometry
        A geometria que define a região de interesse para onde as imagens serão
        recortadas. A geometria pode ser um ponto, linha, polígono, etc.

    Returns
    -------
    ee.ImageCollection
        Uma nova coleção de imagens com cada imagem recortada para a região
        definida pela geometria fornecida.
    """

    def clip_func(image):
        return image.clip(geometria)

    return dataset.map(clip_func)

[ ]:

# Faz o clip da região desejada em toda coleção.
clipped_dataset = clip_img_collection(dataset, roi)
clipped_dataset

Remoção de valores inválidos#

Após a remoção das nuvens, para melhorar ao conjunto de dados, vamos eliminar valores discrepantes.

Sabemos que a reflectância deve ser menor do que 1, então vamos criar uma máscara para eliminar valores maiores do que 1.

As máscaras são geradas a partir do comparador lt() (less than - menor do que) do GEE.

[ ]:

# Função para remoção de valores inválidos
def remove_valores_invalidos(collection, bandas):
    """Aplica uma máscara para remover valores inválidos em uma coleção de
    imagens.

    Os valores considerados inválidos são aqueles maiores que 1 para as bandas
    especificadas.

    Parameters
    ----------
    collection : ee.ImageCollection
        A coleção de imagens a ser processada.
    bandas : list of str
        Lista com os nomes das bandas que serão processadas.

    Returns
    -------
    ee.ImageCollection
        A coleção de imagens com valores inválidos removidos.
    """

    def mask_out_of_range_reflectance(imagem):
        """Mascara valores de reflectância acima de 1"""
        for banda in bandas:
            imagem = imagem.updateMask(imagem.select(banda).lt(1))
        return imagem

    return collection.map(mask_out_of_range_reflectance)

[ ]:

# Executa a função de remoção dos valores inválidos na coleção.
clipped_dataset = remove_valores_invalidos(clipped_dataset, ["B2", "B3", "B4", "B8"])
clipped_dataset

Composição e Mosaico#

Para se criar um mosaico a partir de uma coleção de imagens, basta invocar o método mosaic() da coleção de interesse.

Já para a criação de uma composição, um método do tipo reducer que faz um cálculo pixel a pixel deve ser chamado a partir do objeto que contém a coleção. Nesse caso, os pixels transparentes são ignorados durante os cálculos.

Os cálculos são feitos de forma separada para cada banda das imagens da coleção, resultando assim em uma imagem com as mesmas bandas da coleção, onde cada pixel de cada banda é preenchido por um valor calculado a partir da banda de mesmo nome na coleção.

[ ]:

# Criação de mosaico a partir da coleção.
mosaico = clipped_dataset.mosaic()
# TODO: não vamos usar o mosaico pra nada? Roteiro antigo não ajudou muito.

# Criação de composição a partir da coleção.
composition = clipped_dataset.mean()

Vamos conferir o resultado da composição através do método bandTypes() que retorna as bandas da imagem.

[ ]:

composition.bandTypes()
# TODO: for a strange reason, it is giving me numbers higher than 1.0, needs checking

Calcular NDVI#

Agora queremos calcular a banda de índice de vegetação (NDVI) para cada imagem da coleção. Vamos seguir os mesmos passos do laboratório anterior!

Alguns comentários:

B2 é a banda azul
B3 é a banda verde
B4 é a banda vermelha
B8 é a banda infravermelho próximo (NIR)
MSK_CLDPRB é a banda de probabilidade de nuvens
SCL é fruto de um processo de classificação pelo algoritmo de classificação de cena (Scene Classification)

💡 Sempre verifique a documentação oficial: https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR_HARMONIZED#bands

[ ]:

# calcular a banda NDVI conforme visto no laboratório anterior.

img_bd_nir = composition.select("B8")
img_bd_red = composition.select("B4")

## adicionar a banda NDVI à coleção de imagens
composition = composition.addBands(
    img_bd_nir.subtract(img_bd_red).divide(img_bd_nir.add(img_bd_red)).rename("NDVI"),
    ["NDVI"],
)
composition

Seleção de bandas para a classificação#

Vamos criar listas com nomes das bandas para utilizarmos ao longo deste laboratório.

[ ]:

# Lista com somente as bandas de interesse para a classificação
classification_bands = ["B2", "B3", "B4", "B8", "NDVI"]

[ ]:

# Descarta as bandas usadas para se eliminar efeitos de nuvens.
composition = composition.select(classification_bands)
composition

Classificação supervisionada#

Definindo amostras conhecidas#

Vamos começar a definir visualmente regiões conhecidas para treinamento do modelo de classificação. Em outras palavras, precisamos gerar amostras sobre a imagem para que o modelo aprenda a classificar os pixels.

Essa é a magia da classificação supervisionada… O ser humano aponta para alguns polígonos e diz para o computador: «olha, isso aqui é água, isso aqui é vegetação, isso aqui é solo exposto, etc». Em seguida, o computador aprende a classificar os pixels da imagem de acordo com o que foi definido pelo ser humano.

Essa é uma etapa tão especial que merece ser apresentada pelo próprio Dr. Qiusheng Wu, criador do geemap e Prof. da Universidade do Tennessee, nos EUA. Veja:

[ ]:

import geemap

# TODO: the process described in the video is temporally broken
# see: https://github.com/gee-community/geemap/discussions/1816
geemap.show_youtube("https://youtu.be/VWh5PxXPZw0")

Vamos definir um dicionário para armazenar os parâmetros de visualização básicos para o mapa.

[ ]:

visualization = {
    "min": 0.0,
    "max": 0.3,
    "bands": ["B4", "B3", "B2"],
}

Agora de fato começaremos a definir as regiões de treinamento. Devido a um erro no geemap, vamos coletar classe a classe e depois juntar tudo em um único FeatureCollection.

Você deve seguir os seguintes passos:

Execute a célula que começa com my_map = geemap.Map(
Clique no botão toolbar no canto superior direito do mapa
Clique no botão expand tool bar (sinal de +)
Clique no botão collect training samples (símbolo de um dedo indicador)
Altere o campo Required para classe
Altere o campo Integer para um número inteiro de sua escolha. Recomenda-se começar com 0 e ir aumentando de 1 em 1 para cada classe. É extremamente importante que cada classe receba um número diferente neste campo, por exemplo, 0=agua, 1=vegetacao, 2=solo exposto, etc.
Altere o campo Optional para label
Altere o campo String para o nome da classe, por exemplo, agua, vegetacao, solo exposto, etc.
Clique no botão Draw a polygon (símbolo de um hexágono) no canto esquerdo do mapa
Desenhe um polígono sobre a imagem que aparece no mapa.
Ao desenhar todos seus polígonos, execute a célula que vem logo em seguida a esta. Não volte a executar a mesma célula novamente, pois isso irá apagar todas as suas amostras de treinamento.

[ ]:

# executar esta célula e desenhar polígonos de AGUA
my_map = geemap.Map(lite_mode=False)
my_map.set_center(-46.711, -23.641, 12)
my_map.addLayer(composition, visualization, "mosaico")
my_map

[ ]:

# executar esta célula somente após desenhar todos os polígonos de AGUA
agua = my_map.user_rois
agua

[ ]:

# executar esta célula e desenhar polígonos de VEGETACAO RASTEIRA
my_map = geemap.Map(lite_mode=False)
my_map.set_center(-46.711, -23.641, 12)
my_map.addLayer(composition, visualization, "mosaico")
my_map

[ ]:

# executar esta célula somente após desenhar todos os polígonos de VEGETACAO RASTEIRA
vegeta_baixa = my_map.user_rois
vegeta_baixa

[ ]:

# executar esta célula e desenhar polígonos de VEGETACAO ALTA
my_map = geemap.Map(lite_mode=False)
my_map.set_center(-46.711, -23.641, 12)
my_map.addLayer(composition, visualization, "mosaico")
my_map

[ ]:

# executar esta célula somente após desenhar todos os polígonos de VEGETACAO ALTA
vegeta_alta = my_map.user_rois
vegeta_alta

[ ]:

# executar esta célula e desenhar polígonos de CONSTRUÇÕES
my_map = geemap.Map(lite_mode=False)
my_map.set_center(-46.711, -23.641, 12)
my_map.addLayer(composition, visualization, "mosaico")
my_map

[ ]:

# executar esta célula somente após desenhar todos os polígonos de CONSTRUÇÕES
construcoes = my_map.user_rois
construcoes

[ ]:

# executar esta célula e desenhar polígonos de SOLO EXPOSTO
my_map = geemap.Map(lite_mode=False)
my_map.set_center(-46.711, -23.641, 12)
my_map.addLayer(composition, visualization, "mosaico")
my_map

[ ]:

# executar esta célula somente após desenhar todos os polígonos de SOLO EXPOSTO
solo_exposto = my_map.user_rois
solo_exposto

Mesclando amostras em uma única feature collection#

Os objetos agua, vegeta_baixa, vegeta_alta, construcoes e solo_exposto são do tipo ee.FeatureCollection. Vamos mesclar todos em uma única feature collection para facilitar a manipulação dos dados. Utilizaremos o método merge() do GEE.

[ ]:

# Criação de um único objeto do tipo ee.FeatureCollection com todos os objetos.
collection = (
    agua.merge(vegeta_baixa).merge(vegeta_alta).merge(construcoes).merge(solo_exposto)
)
collection

Sampling#

O procedimento anterior apenas delimitou as regiões de onde se deseja coletar amostra, associando-as com suas respectivas classes, e não colheu nenhuma amostra propriamente.

Vamos fazer a coleta dos valores de pixels nas regiões.

Adicionalmente, vamos definir a escala que queremos trabalhar, nesse caso 10 metros pois é o menor valor (mais preciso) disponível para o Sentinel-2.

[ ]:

# Ajuste este valor para aumentar ou diminuir a área de amostragem
SCALE = 10  # meters

[ ]:

# Extrai o valor de todos os pixels nas regiões amostradas.
amostra_treinamento = composition.sampleRegions(
    collection=collection, properties=["classe"], scale=SCALE
)

[ ]:

amostra_treinamento

Treinamento#

Para treinar o modelo, basta utilizar um classificador com os parâmetros desejados e utilizar o método train().

Vamos utilizar o classificador Random Forest.

O primeiro parâmetro do train() é a amostra a ser utilizada para treinar o modelo, o segundo é a propriedade que contém o número que identifica a classe e o último, a lista com as bandas que serão utilizadas para classificação.

Na instanciação do classificador, smileRandomForest(), o parâmetro é a quantidade de árvores de decisão que devem ser criadas. Esse valor pode, e deve, ser ajustado de acordo com os resultados no conjunto de validação.

Outros parâmetros podem ser encontrados na documentação do GEE e a explicação da utilidade dos mesmos em literatura apropriada.

[ ]:

# Instancia um classificador na memória com os parâmetros dados e treinando no conjunto de treinamento.
classificador_treinado = ee.Classifier.smileRandomForest(50).train(
    amostra_treinamento, "classe", classification_bands
)
classificador_treinado

Um primeiro teste pode ser feito com a própria amostra de treinamento, com a criação da matriz de confusão a partir do método confusionMatrix().

A acurácia provavelmente será altíssima, visto que foi o próprio conjunto utilizado no treinamento que foi avaliado.

Caso o desempenho tenha sido ruim sobre o conjunto de treinamento, é um provável caso de underfitting.

[ ]:

# // Matriz de confusão do conjunto de treinamento e acurácia.
matriz_confusao = classificador_treinado.confusionMatrix()
matriz_confusao

[ ]:

print(
    f"Acurácia no conjunto de treinamento: {matriz_confusao.accuracy().getInfo()*100:.4f}%"
)

Classificação#

Para classificar a imagem a partir do treinamento, invoca-se o método classify() a partir da imagem e com o classificador treinado escolhido como parâmetro.

O resultado é um objeto do tipo imagem com uma única banda em que os pixels armazenam o valor relativo às classes que foram atribuídas no momento do desenho das regiões das amostras em tela.

[ ]:

# Classifica a imagem com o classificador treinado e com os parâmetros definidos.
imagem_classificada = composition.classify(classificador_treinado)
imagem_classificada

Visualização#

Para visualizar a imagem, primeiro uma paleta de cores deve ser definida.

Essa paleta deve ser feita como uma lista de strings em que cada cor é uma sequência de seis dígitos em hexadecimal.

A ordem em que as cores estão na string representa o número inteiro associado com a classe a partir de 0.

[ ]:

# Paleta de cores para adicionar a imagem classificada ao mapa.
paleta_cores_classes = [
    "#1f77b4",  # Água - um azul mais suave e claro.
    "#98df8a",  # Vegetação rasteira - um verde claro para diferenciar da vegetação alta.
    "#2ca02c",  # Vegetação alta - um verde mais vibrante e menos saturado.
    "#7f7f7f",  # Construção - um cinza médio que representa áreas construídas.
    "#ff7f0e",  # Solo exposto - um laranja mais vibrante e atraente.
]

Finalmente, podemos visualizar a imagem classificada.

[ ]:

final_map = geemap.Map(lite_mode=False)
final_map.centerObject(roi, 12)
vis_params = {
    "min": 0,
    "max": 4,
    "bands": ["classification"],
    "palette": paleta_cores_classes,
}
final_map.addLayer(
    imagem_classificada.clip(roi), vis_params, "Imagem Classificada", True
)
final_map

Atividade#

Este laboratório não requer uma atividade específica. Você pode aproveitar esse tempo para aplicar os conceitos no seu trabalho prático.

[ ]:

# Exportar a imagem classificada para um shapefile
geemap.ee_export_vector(
    ee_object=collection,
    filename="collection.shp",
    selectors=None,
    verbose=True,
    keep_zip=True,
    timeout=300,
    proxies=None,
)

[ ]:

# %pip install pycrs

[ ]:

geemap.shp_to_ee("collection.shp")

Laboratório 4

Conteúdo