嗨，在这篇博客中，我们将讨论如何使用 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中文网其他相关文章！
