嗨，在这篇博客中，我们将讨论如何使用 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/

运行存储库中提供的脚本以从 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

分析

让我们创建一个函数来获取多边形中的 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 geom 上创建索引以有效地可视化它

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/