使用 Uber hndexes 和 PostgreSQL 進行柵格分析

嗨，在這篇部落格中，我們將討論如何使用 h3 索引輕鬆進行柵格分析。

客觀的

為了學習，我們將計算出由 ESRI 土地覆蓋確定的聚居區有多少建築物。讓我們針對向量和柵格的國家級資料進行目標。

我們先找到數據

下載柵格數據

我已從 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 大約需要一分鐘

將osm資料插入postgresql表

我們正在使用 osm2pgsql 將 osm 資料插入我們的表中

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

osm2pgsql 總共花了 274 秒（4m 34​​ 秒）。

如果您有使用 ogr2ogr 的文件，也可以使用 geojson 文件

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

ogro2gr 對驅動程式有廣泛的支持，因此您可以非常靈活地選擇輸入內容。輸出是 Postgresql 表

計算h3指數

PostgreSQL

安裝

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/

執行儲存庫中提供的腳本以從柵格建立 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);

拍攝計數時：我所在的地區有 24416 棟建築，其建築等級屬於 ESRI

確認

讓我們驗證輸出是否為真：讓我們以 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 幾何圖形上建立索引以有效地視覺化它

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/

以上是使用 Uber hndexes 和 PostgreSQL 進行柵格分析的詳細內容。更多資訊請關注PHP中文網其他相關文章！