處理多波段柵格（Sentinel-with hndex 並建立索引

嗨，在先前的部落格中，我們討論如何使用 h3 索引和 postgresql 對單波段柵格進行柵格分析。在本部落格中，我們將討論如何處理多波段柵格並輕鬆建立索引。我們將使用 Sentinel-2 影像並從處理後的 h3 細胞創建 NDVI 並視覺化結果

下載哨兵2數據

我們正在從尼泊爾博卡拉地區的https://apps.sentinel-hub.com/eo-browser/下載sentinel 2數據，只是為了確保湖泊在圖像網格中，以便於查看我們驗證NDVI 結果

要下載所有頻段的哨兵影像：

• 您需要建立一個帳戶
• 尋找您所在區域的影像，選擇覆蓋您感興趣區域的網格
• 縮放到網格，然後點選右側垂直條上的圖示
• 之後進入分析標籤並選擇影像格式為 tiff 32 位元、高解析度、wgs1984 格式的所有波段並選取所有波段

您也可以下載預先產生的指數，例如 NDVI、僅假色 tiff 或最適合您需求的特定波段。我們正在下載所有樂隊，因為我們想自己進行處理

• 點選下載

預處理

當我們下載原始格式時，我們將所有樂團作為與哨兵分開的 tiff

• 讓我們建立一個合成影像：

這可以透過GIS工具或gdal來完成

1. 使用 gdal_merge:

我們需要將下載的檔案重新命名為 band1,band2 以避免檔案名稱中出現斜線
本練習最多處理頻段 9，您可以依需求選擇頻段

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 :

• 將所有單獨的波段載入到 QGIS
• 到光柵>雜項>合併

• 合併時，您需要確保勾選「將每個輸入檔案放入 sep band」

• 現在將合併的 tiff 作為複合材料匯出到原始 geotiff

家政

• 確保您的影像採用 WGS1984 在我們的例子中，圖像已經是 ws1984，所以不需要轉換
• 確保您沒有任何 nodata，如果有則用 0 填充

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

• 最後確保你的輸出影像是COG格式

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

我正在使用 cog2h3 儲存庫中提供的 bash 腳本來自動化這些

sudo bash pre.sh sentinel2_composite.tif

h3細胞的處理與創建

現在終於完成了預處理腳本，讓我們繼續計算複合齒輪影像中每個波段的 h3 單元

• 安裝cog2h3

  pip install cog2h3

• 匯出您的資料庫憑證

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

• 奔跑

我們對此哨兵影像使用解析度 10，但您也會在腳本本身中看到，它將列印柵格的最佳分辨率，使 h3 單元小於柵格中的最小像素。

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

我們花了一分鐘來計算結果並將結果儲存在 postgresql

日誌：

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

分析

現在我們的資料已經在postgresql了，讓我們來做一些分析

• 驗證我們是否擁有處理過的所有頻段（記住我們處理的是從頻段 1 到頻段 9）

select *
from sentinel

• 計算每個單元格的 ndvi

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

查詢計畫：

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                                                                                       |

As you can see here for all the rows in that area the calculation is instant . This is true for all other indices and you can compute complex indices join with other tables using the h3_ix primary key and derive meaningful result out of it without worrying as postgresql is capable of handling complex queries and table join.

Visualize and verification

Lets visualize and verify if the computed indices are true

• Create table ( for visualizing in QGIS )

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

• Lets add geometry to visualize the h3 cells This is only necessary to visualize in QGIS , if you build an minimal API by yourself you don't need this as you can construct geometry directly from query

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

• Create index on geometry

create index on ndvi_sentinel(geometry);

• Connect your database in QGIS and visualize the table on the basis of ndvi value Lets get the area near Fewa lake or cloud

As we know value between -1.0 to 0.1 should represent Deep water or dense clouds
lets see if thats true ( making first category as transparent to see the underlying image )

• Check clouds :

• Check Lake

As there were clouds around the lake hence nearby fields are covered by cloud which makes sense

Thank you for reading ! See you in next blog

以上是處理多波段柵格（Sentinel-with hndex 並建立索引的詳細內容。更多資訊請關注PHP中文網其他相關文章！