Olá, No blog anterior falamos sobre como podemos fazer análise raster usando índices h3 e postgresql para um raster de banda única. Neste blog falaremos sobre como podemos processar raster multibanda e criar índices com facilidade. Usaremos a imagem sentinela-2 e criaremos NDVI a partir das células h3 processadas e visualizaremos os resultados

Baixar dados do sentinela 2

Estamos baixando os dados do sentinela 2 de https://apps.sentinel-hub.com/eo-browser/ em Pokhara, região do Nepal, apenas para ter certeza de que o lago está na grade de imagens para que seja fácil para para validarmos o resultado do NDVI

Para baixar a imagem sentinela com todas as bandas:

• Você precisa criar uma conta
• Encontre a imagem na sua área selecione a grade que cobre sua área de interesse
• Amplie a grade e clique no ícone na barra vertical direita
• Depois disso vá na aba analítica e selecione todas as bandas com formato de imagem como tiff 32 bits, alta resolução, formato wgs1984 e todas as bandas marcadas

Você também pode baixar índices pré-gerados, como NDVI, apenas tiff de cor falsa ou bandas específicas, o que melhor atender às suas necessidades. Estamos baixando todas as bandas porque queremos fazer o processamento sozinhos

• Clique em download

Pré-processar

Obtemos todas as bandas como tiff separadas do sentinela conforme baixamos o formato bruto

• vamos criar uma imagem composta:

Isso pode ser feito através de ferramentas GIS ou gdal

1. Usando gdal_merge:

Precisamos renomear o arquivo baixado para band1,band2 assim para evitar barras no nome do arquivo
Vamos processar até a banda 9 para este exercício, você pode escolher a banda conforme sua necessidade

gdal_merge.py -separate -o sentinel2_composite.tif band1.tif band2.tif band3.tif band4.tif band5.tif band6.tif band7.tif band8.tif band9.tif 

1. Usando QGIS :

• Carregar todas as bandas individuais no QGIS
• Vá para Raster > Diversos > Mesclar

• Ao mesclar, você precisa marcar 'colocar cada arquivo de entrada na banda de setembro'

• Agora exporte seu tiff mesclado para geotiff bruto como composto

Limpeza

• Certifique-se de que sua imagem esteja em WGS1984 no nosso caso a imagem já está em ws1984 então não há necessidade de conversão
• Certifique-se de não ter nenhum nodata, se sim, preencha-os com 0

  gdalwarp -overwrite -dstnodata 0 "$input_file" "${output_file}_nodata.tif"

• Finalmente, certifique-se de que sua imagem de saída esteja em COG

  gdal_translate -of COG "$input_file" "$output_file"

Estou usando o script bash fornecido no repositório cog2h3 para automatizá-los

sudo bash pre.sh sentinel2_composite.tif

Processo e criação de células h3

Agora, finalmente, como fizemos o script de pré-processamento, vamos prosseguir para calcular células h3 para cada banda na imagem de engrenagem composta

• Instale cog2h3

  pip install cog2h3

• Exporte suas credenciais de banco de dados

  export DATABASE_URL="postgresql://user:password@host:port/database"

• Correr

Estamos usando a resolução 10 para esta imagem sentinela, no entanto, você também verá no próprio script que imprimirá a resolução ideal para o seu raster, o que torna a célula h3 menor que o menor pixel do raster.

  cog2h3 --cog sentinel2_composite_preprocessed.tif --table sentinel --multiband --res 10

Demoramos um minuto para calcular e armazenar o resultado no postgresql

Registros:

2024-08-24 08:39:43,233 - INFO - Starting processing
2024-08-24 08:39:43,234 - INFO - COG file already exists at sentinel2_composite_preprocessed.tif
2024-08-24 08:39:43,234 - INFO - Processing raster file: sentinel2_composite_preprocessed.tif
2024-08-24 08:39:43,864 - INFO - Determined Min fitting H3 resolution for band 1: 11
2024-08-24 08:39:43,865 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:44,037 - INFO - Resampling Done for band 1
2024-08-24 08:39:44,037 - INFO - New Native H3 resolution for band 1: 10
2024-08-24 08:39:44,738 - INFO - Calculation done for res:10 band:1
2024-08-24 08:39:44,749 - INFO - Determined Min fitting H3 resolution for band 2: 11
2024-08-24 08:39:44,749 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:44,757 - INFO - Resampling Done for band 2
2024-08-24 08:39:44,757 - INFO - New Native H3 resolution for band 2: 10
2024-08-24 08:39:45,359 - INFO - Calculation done for res:10 band:2
2024-08-24 08:39:45,366 - INFO - Determined Min fitting H3 resolution for band 3: 11
2024-08-24 08:39:45,366 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:45,374 - INFO - Resampling Done for band 3
2024-08-24 08:39:45,374 - INFO - New Native H3 resolution for band 3: 10
2024-08-24 08:39:45,986 - INFO - Calculation done for res:10 band:3
2024-08-24 08:39:45,994 - INFO - Determined Min fitting H3 resolution for band 4: 11
2024-08-24 08:39:45,994 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:46,003 - INFO - Resampling Done for band 4
2024-08-24 08:39:46,003 - INFO - New Native H3 resolution for band 4: 10
2024-08-24 08:39:46,605 - INFO - Calculation done for res:10 band:4
2024-08-24 08:39:46,612 - INFO - Determined Min fitting H3 resolution for band 5: 11
2024-08-24 08:39:46,612 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:46,619 - INFO - Resampling Done for band 5
2024-08-24 08:39:46,619 - INFO - New Native H3 resolution for band 5: 10
2024-08-24 08:39:47,223 - INFO - Calculation done for res:10 band:5
2024-08-24 08:39:47,230 - INFO - Determined Min fitting H3 resolution for band 6: 11
2024-08-24 08:39:47,230 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:47,239 - INFO - Resampling Done for band 6
2024-08-24 08:39:47,239 - INFO - New Native H3 resolution for band 6: 10
2024-08-24 08:39:47,829 - INFO - Calculation done for res:10 band:6
2024-08-24 08:39:47,837 - INFO - Determined Min fitting H3 resolution for band 7: 11
2024-08-24 08:39:47,837 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:47,845 - INFO - Resampling Done for band 7
2024-08-24 08:39:47,845 - INFO - New Native H3 resolution for band 7: 10
2024-08-24 08:39:48,445 - INFO - Calculation done for res:10 band:7
2024-08-24 08:39:48,453 - INFO - Determined Min fitting H3 resolution for band 8: 11
2024-08-24 08:39:48,453 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:48,461 - INFO - Resampling Done for band 8
2024-08-24 08:39:48,461 - INFO - New Native H3 resolution for band 8: 10
2024-08-24 08:39:49,046 - INFO - Calculation done for res:10 band:8
2024-08-24 08:39:49,054 - INFO - Determined Min fitting H3 resolution for band 9: 11
2024-08-24 08:39:49,054 - INFO - Resampling original raster to: 200.786148m
2024-08-24 08:39:49,062 - INFO - Resampling Done for band 9
2024-08-24 08:39:49,063 - INFO - New Native H3 resolution for band 9: 10
2024-08-24 08:39:49,647 - INFO - Calculation done for res:10 band:9
2024-08-24 08:39:51,435 - INFO - Converting H3 indices to hex strings
2024-08-24 08:39:51,906 - INFO - Overall raster calculation done in 8 seconds
2024-08-24 08:39:51,906 - INFO - Creating or replacing table sentinel in database
2024-08-24 08:40:03,153 - INFO - Table sentinel created or updated successfully in 11.25 seconds.
2024-08-24 08:40:03,360 - INFO - Processing completed

Analisar

Já que agora temos nossos dados no postgresql, vamos fazer algumas análises

• Verifique se temos todas as bandas que processamos (lembre-se de que processamos das bandas 1 a 9)

select *
from sentinel

• Calcular ndvi para cada célula

explain analyze 
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel

Plano de consulta:

QUERY PLAN                                                                                                       |
----------------------------------------------------------------------------------------------------------------- 
Seq Scan on sentinel  (cost=0.00..28475.41 rows=923509 width=16) (actual time=0.014..155.049 rows=923509 loops=1)|
Planning Time: 0.080 ms                                                                                          |
Execution Time: 183.764 ms                                                                                       |

Como você pode ver aqui, para todas as linhas nessa área, o cálculo é instantâneo. Isso é verdade para todos os outros índices e você pode calcular a junção de índices complexos com outras tabelas usando a chave primária h3_ix e obter resultados significativos sem se preocupar, pois o postgresql é capaz de lidar com consultas complexas e junção de tabelas.

Visualizar e verificar

Vamos visualizar e verificar se os índices calculados são verdadeiros

Criar tabela (para visualização no QGIS)

criar tabela ndvi_sentinel como( selecione h3_ix , (band8-band4)/(band8 band4) como ndvi de public.sentinel )

create table ndvi_sentinel
as(
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel )

Vamos adicionar geometria para visualizar as células h3 Isso só é necessário para visualizar no QGIS, se você construir uma API mínima sozinho, não precisará disso, pois poderá construir geometria diretamente da consulta

ALTER TABLE ndvi_sentinel ADICIONAR geometria geométrica da COLUNA (Polígono, 4326) GERADO SEMPRE COMO (h3_cell_to_boundary_geometry(h3_ix)) ARMAZENADO;

create table ndvi_sentinel
as(
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel )

Criar índice na geometria

criar índice em ndvi_sentinel(geometry);

create table ndvi_sentinel
as(
select h3_ix , (band8-band4)/(band8 band4) as ndvi
from public.sentinel )

Conecte seu banco de dados no QGIS e visualize a tabela com base no valor ndvi Vamos chegar à área perto do lago ou nuvem Fewa

Como sabemos, valores entre -1,0 e 0,1 devem representar águas profundas ou nuvens densas

vamos ver se isso é verdade (tornando a primeira categoria tão transparente para ver a imagem subjacente)

Verifique as nuvens:

Verificar Lago

Como havia nuvens ao redor do lago, os campos próximos estão cobertos por nuvens, o que faz sentido

Obrigado por ler! Vejo você no próximo blog