こんにちは。このブログでは、h3 インデックスを使用してラスター解析を簡単に行う方法について説明します。

客観的

学習のために、ESRI Land Coverで定められた定住エリアに何棟の建物があるかを計算してみます。ベクターとラスターの両方で国家レベルのデータを目指しましょう。

まずはデータを探しましょう

ラスターデータのダウンロード

Esri Land Cover から決済エリアをダウンロードしました。

• https://livingatlas.arcgis.com/landcover/

2023 年をダウンロードしましょう。サイズは約 362MB です

ネパールの OSM 建物をダウンロード

出典 : http://download.geofabrik.de/asia/nepal.html

wget http://download.geofabrik.de/asia/nepal-latest.osm.pbf

データの前処理

実際の h3 セルの計算の前にデータに前処理を適用しましょう
このステップでは gdal コマンドライン プログラムを使用します。マシンに gdal をインストール

クラウドに最適化された Geotiff への変換

cog を知らない場合は、ここでチェックアウトしてください: https://www.cogeo.org/

• gdal_translate が利用可能かどうかを確認する

gdal_translate --version

使用している gdal のバージョンが出力されるはずです

• ラスターを 4326 に再投影

ラスターには異なるソース SRS がある可能性があります。それに応じて変更してください

gdalwarp esri-settlement-area-kathmandu-grid.tif esri-landcover-4326.tif -s_srs EPSG:32645 -t_srs EPSG:4326

• 次に、tif をクラウドに最適化された geotif に変換しましょう

gdal_translate -of COG esri-landcover-4326.tif esri-landcover-cog.tif

再投影された tiff を geotiff に変換するのに約 1 分かかりました

OSM データを postgresql テーブルに挿入する

osm2pgsql を使用して osm データをテーブルに挿入しています

osm2pgsql --create nepal-latest.osm.pbf -U postgres

osm2pgsql には全体で 274 秒 (4 分 34 秒) かかりました。

ogr2ogr
を使用している場合は、geojson ファイルも使用できます

ogr2ogr -f PostgreSQL  PG:"dbname=postgres user=postgres password=postgres" buildings_polygons_geojson.geojson -nln buildings

ogro2gr はドライバーを幅広くサポートしているため、入力内容についてかなり柔軟に対応できます。出力はPostgresqlテーブル

h3 インデックスを計算する

ポストグレSQL

インストール

pip install pgxnclient cmake
pgxn install h3

データベースに拡張機能を作成

create extension h3;
create extension h3_postgis CASCADE;

次に、建物テーブルを作成しましょう

CREATE TABLE buildings (
  id SERIAL PRIMARY KEY,
  osm_id BIGINT,
  building VARCHAR,
  geometry GEOMETRY(Polygon, 4326)
);

テーブルにデータを挿入

INSERT INTO buildings (osm_id, building, geometry)
SELECT osm_id, building, way
FROM planet_osm_polygon pop
WHERE building IS NOT NULL;

ログとタイミング:

Updated Rows    8048542
Query   INSERT INTO buildings (osm_id, building, geometry)
    SELECT osm_id, building, way
    FROM planet_osm_polygon pop
    WHERE building IS NOT NULL
Start time  Mon Aug 12 08:23:30 NPT 2024
Finish time Mon Aug 12 08:24:25 NPT 2024

次に、 centroid を使用してこれらの建物の h3 インデックスを計算しましょう。ここで 8 は私が取り組んでいる h3 解像度です。解像度の詳細については、こちらをご覧ください

ALTER TABLE buildings ADD COLUMN h3_index h3index GENERATED ALWAYS AS (h3_lat_lng_to_cell(ST_Centroid(geometry), 8)) STORED;

ラスター操作

インストール

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

再投影された歯車が静的であることを確認してください/

mv esri-landcover-cog.tif ./static/

リポジトリで提供されているスクリプトを実行して、 raster から h3 セルを作成します。私はモード方式でリサンプリングしています。これは、持っているデータのタイプによって異なります。カテゴリカル データの場合、モードの適合性が高くなります。リサンプリング方法の詳細については、こちらをご覧ください

python cog2h3.py --cog esri-landcover-cog.tif --table esri_landcover --res 8 --sample_by mode

ログ:

2024-08-12 08:55:27,163 - INFO - Starting processing
2024-08-12 08:55:27,164 - INFO - COG file already exists: static/esri-landcover-cog.tif
2024-08-12 08:55:27,164 - INFO - Processing raster file: static/esri-landcover-cog.tif
2024-08-12 08:55:41,664 - INFO - Determined Min fitting H3 resolution: 13
2024-08-12 08:55:41,664 - INFO - Resampling original raster to : 1406.475763m
2024-08-12 08:55:41,829 - INFO - Resampling Done
2024-08-12 08:55:41,831 - INFO - New Native H3 resolution: 8
2024-08-12 08:55:41,967 - INFO - Converting H3 indices to hex strings
2024-08-12 08:55:42,252 - INFO - Raster calculation done in 15 seconds
2024-08-12 08:55:42,252 - INFO - Creating or replacing table esri_landcover in database
2024-08-12 08:55:46,104 - INFO - Table esri_landcover created or updated successfully in 3.85 seconds.
2024-08-12 08:55:46,155 - INFO - Processing completed

分析

ポリゴン内の get_h3_indexes 関数を作成しましょう

CREATE OR REPLACE FUNCTION get_h3_indexes(shape geometry, res integer)
  RETURNS h3index[] AS $$
DECLARE
  h3_indexes h3index[];
BEGIN
  SELECT ARRAY(
    SELECT h3_polygon_to_cells(shape, res)
  ) INTO h3_indexes;

  RETURN h3_indexes;
END;
$$ LANGUAGE plpgsql IMMUTABLE;

対象エリア内で建築エリアとして識別されている建物をすべて取得しましょう

WITH t1 AS (
  SELECT *
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT *
FROM buildings bl
JOIN t1 ON bl.h3_ix = t1.h3_ix;

クエリプラン:

建物の h3_ix 列にインデックスを追加すると、これをさらに強化できます

create index on buildings(h3_ix);

撮影時のカウント: ESRI から分類された建築クラスを持つ建物が私のエリアに 24416 棟ありました

検証

出力が true かどうかを確認してみましょう : 建物を geojson として取得しましょう

WITH t1 AS (
  SELECT *
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT jsonb_build_object(
  'type', 'FeatureCollection',
  'features', jsonb_agg(ST_AsGeoJSON(bl.*)::jsonb)
)
FROM buildings bl
JOIN t1 ON bl.h3_ix = t1.h3_ix;

h3 セルも取得しましょう

with t1 as (
  SELECT *, h3_cell_to_boundary_geometry(h3_ix)
  FROM esri_landcover el
  WHERE h3_ix = ANY (
    get_h3_indexes(
      ST_GeomFromGeoJSON('{
        "coordinates": [
          [
            [83.72922006065477, 28.395029869336483],
            [83.72922006065477, 28.037312312532066],
            [84.2367635433626, 28.037312312532066],
            [84.2367635433626, 28.395029869336483],
            [83.72922006065477, 28.395029869336483]
          ]
        ],
        "type": "Polygon"
      }'), 8
    )
  ) AND cell_value = 7
)
SELECT jsonb_build_object(
  'type', 'FeatureCollection',
  'features', jsonb_agg(ST_AsGeoJSON(t1.*)::jsonb)
)
FROM t1

h3 解像度を上げた後に精度を上げることができますが、入力とリサンプリング手法にも依存します

掃除

不要なテーブルを削除

drop table planet_osm_line;
drop table planet_osm_point;
drop table planet_osm_polygon;
drop table planet_osm_roads;
drop table osm2pgsql_properties;

オプション - ベクター タイル

タイルを視覚化するには、pg_tileserv

• pg_tileserv をダウンロード 上記のリンクからダウンロードするか、docker
を使用してください
• 認証情報をエクスポートする

export DATABASE_URL=postgresql://postgres:postgres@localhost:5432/postgres

• テーブル選択権限を付与します

GRANT SELECT ON buildings to postgres;
GRANT SELECT ON esri_landcover to postgres;

• 視覚化のために h3 インデックスにジオメトリを作成しましょう (st_asmvt からタイルを手動で構築している場合は、クエリから直接これを実行できます)

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

• その h3 geom にインデックスを作成して効果的に視覚化します

CREATE INDEX idx_esri_landcover_geometry 
ON esri_landcover 
USING GIST (geometry);

• 走る

  ./pg_tileserv

• これで、タイルの値に基づいて、または任意の方法でこれらの MVT タイルを視覚化できるようになりました。たとえば、maplibre !処理結果を視覚化したり、他のデータセットと組み合わせたりすることもできます。 これは、QGIS の cell_value に基づいて h3 インデックスに対して行った視覚化です。

ソースリポジトリ: https://github.com/kshitijrajsharma/raster-analysis-using-h3

参考文献：

• https://blog.rustprooflabs.com/2022/04/postgis-h3-intro
• https://jsonsingh.com/blog/uber-h3/
• https://h3ronpy.readthedocs.io/en/latest/