Uber hndexes と PostgreSQL を使用したラスター解析-Python チュートリアル-php.cn

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

客観的

学習のために、ESRI 土地被覆によって決定される居住エリアに何棟の建物があるかを計算します。ベクターとラスターの両方で国家レベルのデータを目指しましょう。

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

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

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

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

サイズ約 362MB の 2023 年をダウンロードしてみましょう

Raster Analysis Using Uber hndexes and PostgreSQL

ネパールの 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;

ログイン後にコピー

クエリプラン:

Raster Analysis Using Uber hndexes and PostgreSQL

建物の 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

ログイン後にコピー

Raster Analysis Using Uber hndexes and PostgreSQL

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

ログイン後にコピー

Raster Analysis Using Uber hndexes and PostgreSQL

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