database - 使用 PostGIS 从点数据生成热图/密度图

Question

我有一个巨大的表 (gps_points)，其中包含存储二维点的几何列。我想要完成的是运行一个输出类似

的查询

 id | freq
-------------
  1 | 365
  2 | 1092
  3 | 97
...

其中“id”是我的总边界框内的一个小矩形的唯一标识符，“freq”是落在该特定矩形内的点数。

所以我定义了一个 PostGIS 表:

create table sub_rects (
id int,
geom geometry)

然后我在外部运行一个脚本，在那里我生成 1000x1000 个这样的矩形并创建它们的多边形，所以我得到了一百万行这样的代码:

insert into sub_rects values(1,ST_GeomFromText('POLYGON((1.1 1.2, 1.1 1.4, 1.5 1.4, 1.5 1.2, 1.1 1.2))'));

当然除了每个多边形都会获得一组新的坐标以匹配其在我的 gps 数据的边界框坐标上的 1000x1000 网格中的实际位置，并且每个元组的 ID 都会更新。

然后我在这个表上生成一个空间索引和一个主键索引。

最后我可以运行这个表和我的原始数据表 (gps_points)

select id, count(*) from sub_rects r join gps_points g on r.geom && g.geom group by id;

这给了我想要的输出。问题是加载所有的小多边形需要很长时间，每次我想生成一个具有不同数量矩形的 map 或运行具有不同底层坐标的数据集时，我必须删除 sub_rects 并生成和重新加载它。

有更好的方法吗？我不需要图形输出。我只需要生成数据。不必在外部生成支持表 (sub_rects) 会非常好，而且我怀疑有一些计算成本更低的方法可以完成同样的事情。我更希望不必安装任何额外的软件。

预计到达时间:根据评论中的要求，这里是查询计划(在我家的机器上，所以较小的数据集和其他表名，但计划相同):

gisdb=# explain analyse select g.id id, count(*) from gridrect g join broadcast b on g.geom && b.wkb_geometry group by g.id;

    QUERY PLAN
-------------------------------------------------------------------------------------------------------------------------------------------------------------------
 GroupAggregate  (cost=0.57..177993.58 rows=10101 width=12) (actual time=14.740..3528.600 rows=1962 loops=1)
   Group Key: g.id
   ->  Nested Loop  (cost=0.57..144786.36 rows=6621242 width=4) (actual time=13.948..3050.741 rows=1366376 loops=1)
         ->  Index Scan using gridrect_id_idx on gridrect g  (cost=0.29..485.30 rows=10201 width=124) (actual time=0.079..6.582 rows=10201 loops=1)
         ->  Index Scan using broadcast_wkb_geometry_geom_idx on broadcast b  (cost=0.29..12.78 rows=137 width=32) (actual time=0.011..0.217 rows=134 loops=10201)
               Index Cond: (g.geom && wkb_geometry)
 Planning time: 0.591 ms
 Execution time: 3529.320 ms
(8 rows)

预计到达时间 2:

根据答案中的建议，我修改了那里建议的代码:

(SELECT row_number() OVER (ORDER BY geom) id, geom
 FROM (SELECT st_geomfromtext(
                  concat('Polygon((', x || ' ' || y, ',', x + xstep || ' ' || y, ',', x + xstep || ' ' || y + ystep,
                         ',', x || ' ' || y + ystep, ',', x || ' ' || y, '))')) geom
       FROM (SELECT x, y
             FROM (SELECT generate_series(xmin, xmin + xdelta, xstep) x) x,
                  (SELECT generate_series(ymin, ymin + ydelta, ystep) y) y) foo) bar);

其中 xmin、ymin、xdelta、ydelta、xstep 和 ystep 都是由外部脚本计算的，但如果将上述内容包装在函数调用中，也可能作为 Postgres 函数的一部分进行计算。从中生成一个临时表并对其运行查询比我最初做的快两个数量级。

Answer 1

这是从边界框生成网格的示例:

https://gis.stackexchange.com/questions/16374/how-to-create-a-regular-polygon-grid-in-postgis

要生成密度数据，请先尝试创建一个包含所有数据的临时表，然后获取计数。根据我的经验，以下比将所有内容组合成一个查询要快一些:

create temp table rect_points as 
select r.id as rect_id, p.id as point_id 
from sub_rects r, gps_points p
where p.geom && r.geom;

create index idx on rect_points (rect_id);

select rect_id, count(*) from rect_points group by rect_id;