6.3. Scripts para suporte à análise de imagens
Para analisar as imagens relacionadas com as glebas do Sicor, podemos construir scripts em Python e combinar algumas bibliotecas disponíveis, que nos auxiliam muito na manipulação dos dados. Tais bibliotecas servem de suporte para abrir uma imagem, para obter a matriz de pixels, reprojetar, combinar as bandas e visualizar as imagens. As bibliotecas GDAL e Rasterio são muito úteis para facilitar o acesso às imagens nos formatos tradicionais, como por exemplo GeoTIFF.
Os scripts que são criados basicamente manipulam números em matrizes digitais, e a biblioteca NumPy possui muitas funções prontas para facilitar essa manipulação. Além disso, a visualização das imagens e gráficos como histogramas, pode ser feita com ajuda da biblioteca Matplotlib.
Como exemplo, uma imagem do sensor WFI (Wide Field Imager, ou Imageador de Amplo Campo de Visada) do satélite AMAZONIA-1, quando baixada do catálogo oficial de imagens do INPE, possui um conjunto de arquivos na seguinte nomenclatura:
AMAZONIA_1_WFI_20240721_036_019_L4_BAND1.tif
AMAZONIA_1_WFI_20240721_036_019_L4_BAND2.tif
AMAZONIA_1_WFI_20240721_036_019_L4_BAND3.tif
AMAZONIA_1_WFI_20240721_036_019_L4_BAND4.tif
Veja que nesse caso o nome dos arquivos já contém informações importantes sobre a imagem:
AMAZONIA_1 |
nome do satélite |
WFI |
nome do sensor |
20240721 |
data da realização do imageamento, no formato YYYYMMDD |
036_019 |
localização da imagem na grade de imageamento, path e row |
L4 |
nível de processamento da imagem, nesse caso L4 é o melhor produto disponível, sendo 0 o nível de dados brutos, 1 quando foi aplicada correção radiométrica sobre o nível anterior, 2 quando foi aplicada correção de sistema sobre o nível anterior, 3 quando foi aplicada correção geométrica utilizando pontos de controle sobre o nível anterior, e por fim 4 quando a imagem do nível anterior foi ortorretificada (observação: cada satélite possui uma quantidade diferente de níveis de processamento, e formas variadas de representação desses níveis) |
BAND1, … BAND4 |
indica a banda específica, no caso do AMAZONIA-1, BAND1 corresponde ao canal azul, BAND2 ao canal verde, BAND3 ao canal vermelho, e BAND4 ao canal infravermelho (observação: nem sempre a banda número 1 de um sensor vai corresponder ao canal verde, e nem sempre a BAND4 ao canal infravermelho; isso sempre depende da quantidade de bandas que cada sensor possui. |
Nota
Alguns satélites possuem mais de um sensor, e alguns deles podem capturar as mesmas bandas, que serão representadas em números diferentes. Por exemplo, o CBERS-4 possui a BAND5 correspondendo ao canal verde, e a BAND13 também correspondendo ao canal verde, porém cada uma obtida por um sensor diferente (MUX e WFI).
Quando trabalhamos com arquivos separados para cada banda, é preciso carregar cada banda de forma individual, utilizando a biblioteca Rasterio, assumindo que os metadados serão equivalentes. Para iniciar os experimentos de análise de imagens, apresentamos alguns trechos de scripts em Python.
Script 1. Carregar a imagem e construir uma matriz com os pixels da banda carregada.
Solução:
with rio.open('AMAZONIA_1_WFI_20220819_035_020_L4_BAND1.tif') as banda1:
matriz_azul = banda1.read(1)
# as linhas a seguir são para apresentar algumas informações sobre a imagem
resolucao_espacial = banda1.res[0]
linhas = banda1.height
colunas = banda1.width
print(f'as imagens têm resolução espacial de {resolucao_espacial} m, número de linhas: {linhas}, número de colunas: {colunas}')
print(f'área total imageada: {linhas * colunas * resolucao_espacial * resolucao_espacial / 1000000} km²')
with rio.open('AMAZONIA_1_WFI_20220819_035_020_L4_BAND2.tif') as banda2:
matriz_verde = banda2.read(1)
with rio.open('AMAZONIA_1_WFI_20220819_035_020_L4_BAND3.tif') as banda3:
matriz_vermelho = banda3.read(1)
with rio.open('AMAZONIA_1_WFI_20220819_035_020_L4_BAND4.tif') as banda4:
matriz_iv = banda4.read(1)
Cada par de linhas do código acima solicita a abertura do arquivo de imagem indicado pelo nome, armazenando-o na variável correspondente, banda1, banda2, etc. Em seguida, a função read com o parâmetro 1 indica que a primeira banda contida no arquivo, que nesse caso só possui uma banda, será carregada e a matriz de pixels será armazenada nas variáveis matriz_verde, matriz_azul, matriz_vermelho, matriz_iv.
Script 2. Realizar a visualização das imagens, utilizando as funções disponíveis na biblioteca matplotlib, módulo pyplot.
Solução:
plt.figure(figsize=(19, 3))
plt.subplot(141)
plt.imshow(matriz_azul, cmap = 'gray')
plt.title('Canal Azul')
plt.subplot(142)
plt.imshow(matriz_verde, cmap = 'gray')
plt.title('Canal Verde')
plt.subplot(143)
plt.imshow(matriz_vermelho, cmap = 'gray')
plt.title('Canal Vermelho')
plt.subplot(144)
plt.imshow(matriz_iv, cmap = 'gray')
plt.title('Canal Infravermelho')
plt.colorbar()
plt.show()
A Figura 6.9 abaixo apresenta o resultado do Script 2.
Figura 6.9 - Visualização das 4 bandas carregadas, em tons de cinza.
A visualização dessas imagens está prejudicada por alguns fatores. Primeiramente, a quantidade de pixels de cada banda é muito grande, e o espaço de visualização é muito pequeno. A área representada por esta imagem, considerando a resolução espacial de 64 metros, é de aproximadamente 950 mil km². No entanto, uma gleba armazenada nos dados Sicor, tendo sido escolhida a de maior área no ano de 2022, possui uma área de aproximadamente 15 km². Assim, para focar a visualização dos pixels da imagem somente na área de interesse de uma gleba, podemos recortar a imagem em uma região no entorno da geometria da gleba, considerando um buffer de alguns metros para observação da vizinhança.
A Figura 6.10 abaixo apresenta um recorte das 4 bandas carregadas, em tons de cinza, e com visualizações independentes.
Figura 6.10 - Visualização de um recorte das 4 bandas carregadas, em tons de cinza.
A barra de cores ao lado de cada figura nos indica que essa visualização está com um esquema de cores independente para cada imagem. Esse esquema de cores é baseado nos valores mínimos e máximos encontrados no recorte, e isso automaticamente transforma a visualização numa rampa de cores entre preto e branco. Isso pode nos trazer uma falsa impressão de que o canal infravermelho, por exemplo, possui a resposta espectral semelhante às demais imagens. No entanto podemos perceber pela barra de cores, que os menores valores do canal infravermelho, representados em pixels escuros, correspondem ao número 100, e os maiores valores estão em torno de 500, valor tal muito maior do que nos canais visíveis, como azul, verde e vermelho.
Para observar a distribuição de valores de cada banda, de forma independente, podemos utilizar ferramentas estatísticas. Uma maneira muito utilizada para complementar a visualização das imagens é através do histograma. Ele apresenta, de forma gráfica, a distribuição dos pixels na imagem.
A Figura 6.11 abaixo apresenta 4 histogramas integrados das bandas carregadas. Nesta figura fica bem claro que os pixels do canal infravermelho são muito mais altos do que os demais. Nota-se também a ocorrência dos menores valores no canal vermelho.
Figura 6.11 - Visualização do histograma composto pelas 4 bandas carregadas.
É possível também combinar a visualização de 2 bandas em um gráfico de dispersão, também chamado de scatterplot. Com este gráfico, é possível verificar a diferença, ou a semelhança, entre as bandas. Note como os canais verde e azul possuem uma grande semelhança (o que também se mostrou pela grande sobreposição no histograma).
A Figura 6.12 abaixo apresenta 2 gráficos de dispersão entre pares de bandas, das mesmas 4 bandas carregadas.
Figura 6.12 - Visualização de gráficos de dispersão entre pares de bandas, das mesmas 4 bandas carregadas.
Os intervalos de valores apresentados nos gráficos de dispersão são interessantes pois nos ajudarão posteriormente a identificar formas de melhorar a visualização da imagem numa composição colorida.
6.3.1. Visualização de imagens Sentinel-2 do Brazil Data Cube (BDC)
Podemos utilizar também as imagens do satélite Sentinel-2, disponíveis no projeto Brazil Data Cube (BDC). O BDC é um projeto de pesquisa, desenvolvimento e inovação tecnológica do Instituto Nacional de Pesquisas Espaciais (INPE), que está produzindo dados a partir de grandes volumes de imagens de sensoriamento remoto de média resolução espacial para todo o território nacional. Os dados produzidos no BDC incluem coleções de dados prontos para análise (conhecidos como ARD - Analysis-Ready Data) [76], cubos de dados multidimensionais e mosaicos a partir de imagens dos satélites CBERS-4/4A, Sentinel-2 e Landsat-8.
Um dos produtos disponíveis no acervo do BDC é o resultado da função de composição temporal de imagens Sentinel-2, a qual combina pixels de diversas imagens utilizando técnicas como média, mediana ou best pixel, usando a abordagem Least Cloud-cover First (LCF). A Figura 6.13 apresenta as técnicas utilizadas para gerar esse tipo de produto
Figura 6.13 - Abordagem Least Cloud-cover First (LCF, primeiro os pixels com menor cobertura de nuvens) para obtenção dos melhores pixels.
Script 3. Carregar a geometria de uma gleba, disponível em arquivo Shapefile, e buscar imagens Sentinel-2 no BDC que incluam a área da gleba, em um período de agosto de 2022 até julho de 2023.
Solução:
catalogo_bdc = pystac_client.Client.open("https://data.inpe.br/bdc/stac/v1")
items = catalogo_bdc.search(
collections=['S2-16D-2'],
datetime='2022-08-01/2023-07-31',
intersects=shapely.geometry.mapping(gleba_4326.geometry)
)
total_rasters = items.matched()
print('total de elementos encontrados', total_rasters)
i = 1
for item in items.items():
print(f"{i}: imagem da data {item.properties['datetime']}, link B04 {item.assets['B04'].href}")
i += 1
Para demonstrar os scripts de análise de imagens apresentados nesta seção, veja o seguinte Jupyter Notebook: