嗨，在這篇部落格中，我們將討論如何使用 h3 索引輕鬆進行柵格分析。
為了學習，我們將計算出由 ESRI 土地覆蓋確定的聚居區有多少建築物。讓我們針對向量和柵格的國家級資料進行目標。
我已從 Esri Land Cover 下載了定居點區域。
讓我們下載2023年，大小約362MB
資料來源：http://download.geofabrik.de/asia/nepal.html
wget http://download.geofabrik.de/asia/nepal-latest.osm.pbf
讓我們在實際的 h3 單元計算之前對資料進行一些預處理
我們將在這一步驟中使用 gdal 命令列程式。在你的機器上安裝 gdal
如果您不知道 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
gdal_translate -of COG esri-landcover-4326.tif esri-landcover-cog.tif
將重新投影的 tiff 轉換為 geotiff 大約需要一分鐘
我們正在使用 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 表
安裝
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 快速建立向量圖塊
export DATABASE_URL=postgresql://postgres:postgres@localhost:5432/postgres
GRANT SELECT ON buildings to postgres; GRANT SELECT ON esri_landcover to postgres;
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);
./pg_tileserv
原始碼庫：https://github.com/kshitijrajsharma/raster-analysis-using-h3
以上是使用 Uber hndexes 和 PostgreSQL 進行柵格分析的詳細內容。更多資訊請關注PHP中文網其他相關文章！
