Hola, En el blog anterior hablamos sobre cómo podemos hacer análisis de ráster usando índices h3 y postgresql para un ráster de banda única. En este blog hablaremos sobre cómo podemos procesar ráster multibanda y crear índices con facilidad. Usaremos la imagen Sentinel-2 y crearemos NDVI a partir de las celdas h3 procesadas y visualizaremos los resultados

Descargar datos de Sentinel 2

Estamos descargando los datos de Sentinel 2 de https://apps.sentinel-hub.com/eo-browser/ en Pokhara, área de Nepal, solo para asegurarnos de que el lago esté en la cuadrícula de imágenes para que sea fácil para nosotros para validar el resultado NDVI

Para descargar la imagen centinela con todas las bandas:

• Necesitas crear una cuenta
• Encuentra la imagen en tu área selecciona la grilla que cubre tu Área de interés
• Acerque la cuadrícula y haga clic en el icono en la barra vertical derecha
• Después de eso, vaya a la pestaña analítica y seleccione todas las bandas con formato de imagen como tiff de 32 bits, alta resolución, formato wgs1984 y todas las bandas marcadas

También puede descargar índices pregenerados como NDVI, solo tiff de color falso o bandas específicas, lo que mejor se adapte a sus necesidades. Estamos descargando todas las bandas porque queremos hacer el procesamiento nosotros mismos

• Haga clic en descargar

Preproceso

Obtenemos todas las bandas como tiff separadas del centinela a medida que descargamos el formato sin formato

• vamos a crear una imagen compuesta:

Esto se puede hacer a través de herramientas SIG o gdal

1. Usando gdal_merge:

Necesitamos cambiar el nombre del archivo descargado a band1,band2 así para evitar barras en el nombre del archivo
Procesemos hasta la banda 9 para este ejercicio, puede elegir la banda según sus requisitos

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 :

• Cargar todas las bandas individuales en QGIS
• Vaya a Ráster > Varios > Fusionar

• Al fusionar, debes asegurarte de marcar 'colocar cada archivo de entrada en la banda separada'

• Ahora exporta tu tiff fusionado a geotiff sin formato como compuesto

Gestión interna

• Asegúrate de que tu imagen esté en WGS1984 en nuestro caso, la imagen ya está en ws1984, por lo que no es necesaria la conversión
• Asegúrate de no tener ningún dato. Si es así, complétalo con 0.

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

• Finalmente asegúrate de que tu imagen de salida esté en COG

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

Estoy usando el script bash proporcionado en el repositorio cog2h3 para automatizarlos

sudo bash pre.sh sentinel2_composite.tif

Proceso y creación de células h3.

Ahora, finalmente, como hemos realizado el script de preprocesamiento, avancemos para calcular las celdas h3 para cada banda en la imagen compuesta del engranaje

• Instalar cog2h3

  pip install cog2h3

• Exporta las credenciales de tu base de datos

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

• Correr

Estamos usando la resolución 10 para esta imagen centinela, sin embargo, también verá en el script que imprimirá la resolución óptima para su ráster que hace que la celda h3 sea más pequeña que el píxel más pequeño en el ráster.

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

Nos tomó un minuto calcular y almacenar el resultado en 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

Analizar

Ya que ahora tenemos nuestros datos en postgresql, hagamos un análisis

• Verifique que tengamos todas las bandas que procesamos (recuerde que procesamos de la banda 1 a la 9)

select *
from sentinel

• Calcular ndvi para cada celda

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

Plan 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 puede ver aquí para todas las filas en esa área, el cálculo es instantáneo. Esto es cierto para todos los demás índices y puede calcular índices complejos unidos con otras tablas usando la clave principal h3_ix y obtener resultados significativos sin preocuparse, ya que postgresql es capaz de manejar consultas complejas y uniones de tablas.

Visualizar y verificar

Visualicemos y verifiquemos si los índices calculados son verdaderos

• Crear tabla (para visualizar en QGIS)

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

• Agreguemos geometría para visualizar las celdas h3 Esto solo es necesario para visualizar en QGIS, si construye una API mínima usted mismo, no la necesita, ya que puede construir geometría directamente desde la consulta.

ALTER TABLE ndvi_sentinel  
ADD COLUMN geometry geometry(Polygon, 4326) 
GENERATED ALWAYS AS (h3_cell_to_boundary_geometry(h3_ix)) STORED;

• Crear índice en geometría

create index on ndvi_sentinel(geometry);

• Conecte su base de datos en QGIS y visualice la tabla según el valor de ndvi Consigamos el área cerca del lago o la nube Fewa

Como sabemos, un valor entre -1,0 y 0,1 debería representar aguas profundas o nubes densas
veamos si eso es cierto (haciendo que la primera categoría sea transparente para ver la imagen subyacente)

• Comprueba las nubes:

• Verificar lago

Como había nubes alrededor del lago, los campos cercanos están cubiertos por nubes, lo cual tiene sentido

¡Gracias por leer! Nos vemos en el próximo blog