こんにちは。前のブログでは、単一バンド ラスターに対して h3 インデックスと postgresql を使用してラスター解析を行う方法について説明しました。このブログでは、マルチバンド ラスターを処理し、インデックスを簡単に作成する方法について説明します。 Sentinel-2 画像を使用し、処理された h3 細胞から NDVI を作成し、結果を視覚化します

センチネル 2 データをダウンロードする

ネパール地域ポカラの https://apps.sentinel-hub.com/eo-browser/ からセンチネル 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 バンドに配置する」にチェックを入れる必要があります

• マージされた 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 に保存するのに 1 分かかりました

ログ:

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                                                                                       |

ここでわかるように、その領域のすべての行の計算は瞬時に行われます。これは他のすべてのインデックスに当てはまり、postgresql は複雑なクエリとテーブル結合を処理できるため、h3_ix 主キーを使用して他のテーブルと結合する複雑なインデックスを計算し、心配することなくそこから意味のある結果を導き出すことができます。

見える化と検証

計算されたインデックスが正しいかどうかを視覚化して検証してみましょう

• テーブルの作成 (QGIS で視覚化するため)

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

• h3 セルを視覚化するためにジオメトリを追加しましょう これは QGIS で視覚化する場合にのみ必要です。最小限の API を自分で構築する場合は、クエリから直接ジオメトリを構築できるため、これは必要ありません。

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

• ジオメトリにインデックスを作成

create index on ndvi_sentinel(geometry);

• QGIS でデータベースを接続し、ndvi 値に基づいてテーブルを視覚化します ペワ湖または雲の近くのエリアを取得しましょう

ご存知のとおり、-1.0 ～ 0.1 の値は深海または濃い雲を表すはずです
それが本当かどうか見てみましょう（最初のカテゴリを透明にして、下にある画像が見えるようにします）

• 雲をチェック :

• チェックレイク

湖の周りに雲があったので、近くの畑が雲に覆われているのは当然です

読んでいただきありがとうございます！次のブログでお会いしましょう