嗨，在這篇部落格中，我們將討論如何使用 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 --version
登入後複製
它應該會列印您正在使用的 gdal 版本
您的柵格可能有不同的來源 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;
登入後複製
CREATE INDEX idx_esri_landcover_geometry
ON esri_landcover
USING GIST (geometry);
登入後複製
- 現在您可以根據圖塊值或任何您想要的方式來視覺化這些 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中文網其他相關文章！
