Hallo, im vorherigen Blog haben wir darüber gesprochen, wie wir eine Rasteranalyse mit H3-Indizes und Postgresql für ein einzelnes Bandraster durchführen können. In diesem Blog werden wir darüber sprechen, wie wir Multiband-Raster verarbeiten und ganz einfach Indizes erstellen können. Wir werden das Sentinel-2-Bild verwenden und NDVI aus den verarbeiteten h3-Zellen erstellen und die Ergebnisse visualisieren

Sentinel 2-Daten herunterladen

Wir laden die Sentinel 2-Daten von https://apps.sentinel-hub.com/eo-browser/ in der Region Pokhara, Nepal, herunter, nur um sicherzustellen, dass sich der See im Bildraster befindet, damit er leicht zu erkennen ist uns, um das NDVI-Ergebnis zu validieren

Um das Sentinel-Bild mit allen Bändern herunterzuladen:

• Sie müssen ein Konto erstellen
• Suchen Sie das Bild in Ihrem Bereich. Wählen Sie das Raster aus, das Ihren Interessenbereich abdeckt.
• Zoomen Sie auf das Raster und klicken Sie auf das Symbol in der rechten vertikalen Leiste
• Danach gehen Sie zur Registerkarte „Analyse“ und wählen Sie alle Bänder mit dem Bildformat TIFF 32 Bit, hohe Auflösung, WGS1984-Format aus und markieren Sie alle Bänder.

Sie können auch vorgenerierte Indizes wie NDVI, nur Falschfarben-TIFF oder bestimmte Bänder herunterladen, je nachdem, was Ihren Anforderungen am besten entspricht. Wir laden alle Bands herunter, da wir die Bearbeitung selbst durchführen wollen

• Klicken Sie auf „Herunterladen“

Vorprozess

Wir erhalten alle Bänder als separates TIFF vom Sentinel, während wir das Rohformat heruntergeladen haben

• Erstellen wir ein zusammengesetztes Bild:

Dies kann über GIS-Tools oder gdal erfolgen

1. Verwenden von gdal_merge:

Wir müssen die heruntergeladene Datei wie folgt in band1,band2 umbenennen, um Schrägstriche im Dateinamen zu vermeiden
Lassen Sie uns für diese Übung bis Band 9 vorgehen. Sie können das Band nach Ihren Wünschen auswählen

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

• Alle einzelnen Bänder in QGIS laden
• Gehen Sie zu Raster > Verschiedenes > Zusammenführen

• Beim Zusammenführen müssen Sie sicherstellen, dass Sie die Option „Jede Eingabedatei im Sep-Band platzieren“ aktivieren.

• Exportieren Sie nun Ihr zusammengeführtes TIFF als zusammengesetztes Roh-Geotiff

Hauswirtschaft

• Stellen Sie sicher, dass Ihr Bild in WGS1984 ist In unserem Fall befindet sich das Bild bereits in ws1984, sodass keine Konvertierung erforderlich ist
• Stellen Sie sicher, dass Sie keine Knotendaten haben. Wenn ja, füllen Sie diese mit 0

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

• Stellen Sie abschließend sicher, dass Ihr Ausgabebild im COG ist

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

Ich verwende das im cog2h3-Repo bereitgestellte Bash-Skript, um diese zu automatisieren

sudo bash pre.sh sentinel2_composite.tif

Prozess und Entstehung von h3-Zellen

Nachdem wir nun endlich das Vorverarbeitungsskript erstellt haben, können wir mit der Berechnung der h3-Zellen für jedes Band im zusammengesetzten Zahnradbild fortfahren.

• Installieren Sie cog2h3

  pip install cog2h3

• Exportieren Sie Ihre Datenbankanmeldeinformationen

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

• Laufen

Wir verwenden die Auflösung 10 für dieses Sentinel-Bild, Sie werden jedoch auch im Skript selbst sehen, dass die optimale Auflösung für Ihr Raster gedruckt wird, wodurch die h3-Zelle kleiner als Ihr kleinstes Pixel im Raster wird.

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

Es hat eine Minute gedauert, bis wir das Ergebnis berechnet und in Postgresql gespeichert haben

Protokolle:

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

Analysieren

Da wir jetzt unsere Daten in Postgresql haben, machen wir eine Analyse

• Stellen Sie sicher, dass wir über alle von uns verarbeiteten Bänder verfügen (denken Sie daran, dass wir von Band 1 bis 9 verarbeitet haben).

select *
from sentinel

• Berechnen Sie ndvi für jede Zelle

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

Abfrageplan:

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                                                                                       |

Wie Sie hier sehen können, erfolgt die Berechnung für alle Zeilen in diesem Bereich sofort. Dies gilt für alle anderen Indizes, und Sie können mithilfe des h3_ix-Primärschlüssels komplexe Indizes berechnen und mit anderen Tabellen verknüpfen und daraus aussagekräftige Ergebnisse ableiten, ohne sich Sorgen machen zu müssen, da Postgresql in der Lage ist, komplexe Abfragen und Tabellenverknüpfungen zu verarbeiten.

Visualisieren und überprüfen

Lasst uns visualisieren und überprüfen, ob die berechneten Indizes wahr sind

• Tabelle erstellen (zur Visualisierung in QGIS)

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

• Fügen wir Geometrie hinzu, um die h3-Zellen zu visualisieren Dies ist nur zur Visualisierung in QGIS erforderlich. Wenn Sie selbst eine minimale API erstellen, benötigen Sie dies nicht, da Sie Geometrie direkt aus der Abfrage erstellen können

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

• Index zur Geometrie erstellen

create index on ndvi_sentinel(geometry);

• Verbinden Sie Ihre Datenbank in QGIS und visualisieren Sie die Tabelle auf der Grundlage des ndvi-Werts Lass uns den Bereich in der Nähe des Fewa-Sees oder der Wolke erfassen

Wie wir wissen, sollte ein Wert zwischen -1,0 und 0,1 tiefes Wasser oder dichte Wolken darstellen
Mal sehen, ob das stimmt (die erste Kategorie transparent machen, um das darunter liegende Bild zu sehen)

• Wolken prüfen:

• See prüfen

Da es rund um den See Wolken gab, sind die umliegenden Felder von Wolken bedeckt, was Sinn macht

Vielen Dank fürs Lesen! Wir sehen uns im nächsten Blog