Bonjour, Dans le blog précédent, nous avons expliqué comment effectuer une analyse raster à l'aide des index h3 et de postgresql pour un raster à bande unique. Dans ce blog, nous expliquerons comment traiter les rasters multicanaux et créer facilement des indices. Nous utiliserons l'image sentinelle-2, créerons un NDVI à partir des cellules h3 traitées et visualiserons les résultats

Télécharger les données de Sentinel 2

Nous téléchargeons les données Sentinel 2 depuis https://apps.sentinel-hub.com/eo-browser/ à Pokhara, dans la région du Népal, juste pour nous assurer que le lac est dans la grille d'images afin qu'il soit facile pour nous pour valider le résultat NDVI

Pour télécharger l'image sentinelle avec toutes les bandes :

• Vous devez créer un compte
• Trouvez l'image dans votre zone, sélectionnez la grille qui couvre votre zone d'intérêt
• Zoomez sur la grille, et cliquez sur l'icône dans la barre verticale droite
• Après cela, allez dans l'onglet analytique et sélectionnez toutes les bandes avec un format d'image tiff 32 bits, haute résolution, format wgs1984 et toutes les bandes cochées

Vous pouvez également télécharger des indices prégénérés tels que NDVI, False Color Tiff uniquement ou des bandes spécifiques selon ce qui correspond le mieux à vos besoins. Nous téléchargeons toutes les bandes car nous voulons faire le traitement nous-mêmes

• Cliquez sur Télécharger

Prétraitement

Nous obtenons tous les groupes sous forme de tiff distinct de la sentinelle au fur et à mesure que nous téléchargeons le format brut

• créons une image composite :

Cela peut être fait via des outils SIG ou gdal

1. Utilisation de gdal_merge :

Nous devons renommer le fichier téléchargé en band1,band2 comme ceci pour éviter les barres obliques dans le nom de fichier
Traitons jusqu'à la bande 9 pour cet exercice, vous pouvez choisir la bande selon vos besoins

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. Utiliser QGIS :

• Charger toutes les bandes individuelles dans QGIS
• Allez dans Raster > Divers > Fusionner

• Lors de la fusion, vous devez vous assurer de cocher « placer chaque fichier d'entrée dans la bande septembre »

• Exportez maintenant votre tiff fusionné vers un géotiff brut en tant que composite

Ménage

• Assurez-vous que votre image est en WGS1984 dans notre cas, l'image est déjà dans ws1984 donc pas besoin de conversion
• Assurez-vous de ne pas avoir de nodata si oui remplissez-les avec 0

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

• Enfin, assurez-vous que votre image de sortie est en COG

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

J'utilise le script bash fourni dans le dépôt cog2h3 pour les automatiser

sudo bash pre.sh sentinel2_composite.tif

Processus et création de cellules h3

Maintenant, enfin, après avoir effectué le script de prétraitement, passons au calcul des cellules h3 pour chaque bande de l'image de rouage composite

• Installer cog2h3

  pip install cog2h3

• Exportez vos informations d'identification de base de données

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

• Courir

Nous utilisons la résolution 10 pour cette image sentinelle, mais vous verrez également dans le script lui-même qui imprimera la résolution optimale pour votre raster qui rend la cellule h3 plus petite que votre plus petit pixel dans le raster.

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

Il nous a fallu une minute pour calculer et stocker le résultat dans postgresql

Journaux :

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

Analyser

Puisque maintenant nous avons nos données dans postgresql, faisons quelques analyses

• Vérifiez que nous avons tous les groupes que nous avons traités (n'oubliez pas que nous avons traité du groupe 1 à 9)

select *
from sentinel

• Calculer le ndvi pour chaque cellule

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

Plan de requête :

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                                                                                       |

Comme vous pouvez le voir ici pour toutes les lignes de cette zone, le calcul est instantané. Cela est vrai pour tous les autres indices et vous pouvez calculer la jointure d'index complexes avec d'autres tables à l'aide de la clé primaire h3_ix et en tirer des résultats significatifs sans vous inquiéter car postgresql est capable de gérer des requêtes complexes et des jointures de tables.

Visualisation et vérification

Visualisons et vérifions si les indices calculés sont vrais

• Créer un tableau (pour visualiser dans QGIS)

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

• Ajoutons de la géométrie pour visualiser les cellules h3 Ceci n'est nécessaire que pour visualiser dans QGIS. Si vous créez vous-même une API minimale, vous n'en avez pas besoin car vous pouvez construire la géométrie directement à partir d'une requête.

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

• Créer un index sur la géométrie

create index on ndvi_sentinel(geometry);

• Connectez votre base de données dans QGIS et visualisez le tableau sur la base de la valeur ndvi Obtenons la zone près du lac ou du nuage Fewa

Comme nous le savons, une valeur comprise entre -1,0 et 0,1 devrait représenter des eaux profondes ou des nuages ​​denses
voyons si c'est vrai (rendre la première catégorie aussi transparente pour voir l'image sous-jacente)

• Vérifier les nuages :

• Vérifier le lac

Comme il y avait des nuages ​​autour du lac, les champs voisins sont donc couverts par des nuages, ce qui est logique

Merci d'avoir lu ! Rendez-vous sur le prochain blog