使用 Uber hndexes 和 PostgreSQL 进行栅格分析

2024-08-13 0 573

嗨,在这篇博客中,我们将讨论如何使用 h3 索引轻松进行栅格分析。

客观的

为了学习,我们将计算出由 esri 土地覆盖确定的聚居区有多少建筑物。让我们针对矢量和栅格的国家级数据进行目标。

我们先找到数据

下载栅格数据

我已经从 esri land cover 下载了定居点区域。

  • HTTPs://livingatlas.arcgis.com/landcover/

让我们下载2023年,大小约362mb

使用 Uber hndexes 和 PostgreSQL 进行栅格分析

下载尼泊尔 osm 建筑

来源:http://download.geofabrik.de/asia/nepal.html

2882​​41523188

预处理数据

让我们在实际的 h3 单元计算之前对数据进行一些预处理
我们将在这一步中使用 gdal 命令行程序。在你的机器上安装 gdal

转换为云优化的 geotiff

如果您不知道 cog ,请在此处查看:HTTPS://www.cogeo.org/

  • 检查gdal_translate是否可用

1

gdal_translate --version

它应该打印您正在使用的 gdal 版本

  • 将光栅重新投影为 4326

您的栅格可能有不同的源 srs ,相应地更改它

1

gdalwarp esri-settlement-area-kathmandu-grid.tif esri-landcover-4326.tif -s_srs epsg:32645 -t_srs epsg:4326

  • 现在让我们将 tif 转换为云优化的 geotif

1

gdal_translate -of cog esri-landcover-4326.tif esri-landcover-cog.tif

将重新投影的 tiff 转换为 geotiff 大约需要一分钟

将osm数据插入postgresql表

我们正在使用 osm2pgsql 将 osm 数据插入到我们的表中

1

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

osm2pgsql 总共花费了 274 秒(4m 34​​ 秒)。

如果您有任何使用 ogr2ogr 的文件,也可以使用 geoJSon 文件

1

ogr2ogr -f postgresql  pg:"dbname=postgres user=postgres password=postgres" buildings_polygons_geoJSON.geojson -nln buildings

ogro2gr 对驱动程序有广泛的支持,因此您可以非常灵活地选择输入内容。输出是 postgresql 表

计算h3指数

postgresql

安装

1

2

pip install pgxnclient cmake

pgxn install h3

在您的数据库中创建扩展

1

2

create extension h3;

create extension h3_postgis cascade;

现在让我们创建建筑物表

1

2

3

4

5

6

create Table buildings (

  id serial primary key,

  osm_id bigint,

  building varchar,

  geometry geometry(polygon, 4326)

);

向表中插入数据

1

2

3

4

insert into buildings (osm_id, building, geometry)

select osm_id, building, way

from planet_osm_polygon pop

Where building is not null;

日志和计时

1

2

3

4

5

6

7

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 分辨率。在这里了解有关决议的更多信息

1

alter table buildings add column h3_index h3index generated always as (h3_lat_lng_to_cell(st_centroid(geometry), 8)) stored;

光栅操作

安装

1

pIP install h3 h3ronpy rasterio asyncio asyncpg aiohttp

确保重新投影的齿轮处于静态/

1

mv esri-landcover-cog.tif ./static/

运行 repo 中提供的脚本以从栅格创建 h3 像元。我正在按模式方法重新采样:这取决于您拥有的数据类型。对于分类数据模式更适合。在这里了解有关重采样方法的更多信息

1

Python cog2h3.py --cog esri-landcover-cog.tif --table esri_landcover --res 8 --sample_by mode

日志:

1

2

3

4

5

6

7

8

9

10

11

12

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

1

2

3

4

5

6

7

8

9

10

11

12

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;

让我们获取我们感兴趣的区域中所有被确定为建筑面积的建筑物

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

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;

查询计划:

使用 Uber hndexes 和 PostgreSQL 进行栅格分析

如果在建筑物的 h3_ix 列上添加索引,这可以进一步增强

1

create index on buildings(h3_ix);

拍摄计数时:我所在的地区有 24416 座建筑,其建筑等级属于 esri

确认

让我们验证输出是否为真:让我们以 geojson 形式获取建筑物

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

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细胞

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

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

使用 Uber hndexes 和 PostgreSQL 进行栅格分析

增加 h3 分辨率后可以提高准确性,并且还取决于输入和重采样技术

清理

删除我们不需要的桌子

1

2

3

4

5

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
  • 导出凭证

1

export database_url=postgresql://postgres:postgres@localhost:5432/postgres

  • 授予我们表选择权限

1

2

grant select on buildings to postgres;

grant select on esri_landcover to postgres;

  • 让我们在 h3 索引上创建几何图形以进行可视化(如果您从 st_asmvt 手动构建图块,则可以直接从查询中执行此操作)

1

2

3

alter table esri_landcover

add column geometry geometry(polygon, 4326)

generated always as (h3_cell_to_boundary_geometry(h3_ix)) stored;

  • 在该 h3 几何图形上创建索引以有效地可视化它

1

2

3

create index idx_esri_landcover_geometry

on esri_landcover

using gist (geometry);

  • 奔跑

1

./pg_tileserv

使用 Uber hndexes 和 PostgreSQL 进行栅格分析

  • 现在您可以根据图块值或任何您想要的方式可视化这些 mvt 图块,例如:maplibre!您还可以可视化处理结果或与其他数据集结合。 这是我根据 qgis 中的 cell_value 对 h3 索引进行的可视化 使用 Uber hndexes 和 PostgreSQL 进行栅格分析

代码库: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/
收藏 (0) 打赏

感谢您的支持,我会继续努力的!

打开微信/支付宝扫一扫,即可进行扫码打赏哦,分享从这里开始,精彩与您同在
点赞 (0)

免责声明
1. 本站所有资源来源于用户上传和网络等,如有侵权请邮件联系本站整改team@lcwl.fun!
2. 分享目的仅供大家学习和交流,您必须在下载后24小时内删除!
3. 不得使用于非法商业用途,不得违反国家法律。否则后果自负!
4. 本站提供的源码、模板、插件等等其他资源,都不包含技术服务请大家谅解!
5. 如有链接无法下载、失效或广告,请联系本站工作人员处理!
6. 本站资源售价或VIP只是赞助,收取费用仅维持本站的日常运营所需!
7. 如遇到加密压缩包,请使用WINRAR解压,如遇到无法解压的请联系管理员!
8. 因人力时间成本问题,部分源码未能详细测试(解密),不能分辨部分源码是病毒还是误报,所以没有进行任何修改,大家使用前请进行甄别!
9.本站所有源码资源都是经过本站工作人员人工亲测可搭建的,保证每个源码都可以正常搭建,但不保证源码内功能都完全可用,源码属于可复制的产品,无任何理由退款!

网站搭建学习网 Python 使用 Uber hndexes 和 PostgreSQL 进行栅格分析 https://www.xuezuoweb.com/12815.html

常见问题
  • 本站所有的源码都是经过平台人工部署搭建测试过可用的
查看详情
  • 购买源码资源时购买了带主机的套餐是指可以享受源码和所选套餐型号的主机两个产品,在本站套餐里开通主机可享优惠,最高免费使用主机
查看详情

相关文章

发表评论
暂无评论
官方客服团队

为您解决烦忧 - 24小时在线 专业服务

Fa快捷助手
手机编程软件开发

在手机上用手点一点就能轻松做软件

去做软件
链未云主机
免备案香港云主机

开通主机就送域名的免备案香港云主机

去使用
链未云服务器
免备案香港云服务器

支持售后、超低价、稳定的免备案香港云服务器

去使用