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

我有一个巨大的表（gps_points），其中有一个存储2D点的几何列。我试图做到的是运行一个查询输出类似使用PostGIS从点数据生成热量/密度图

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

其中“ID”是我的总边框和“频率”里面的小矩形的唯一标识符是点的数量落在特定的矩形内。

所以我定义了一个PostGIS的表：

create table sub_rects (
id int, 
geom geometry)

我那么外部运行的脚本，我在那里生成1000×1000这样的矩形和创建它们的多边形，所以我得到这样一个万行：

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))'));

当然，除了每个多边形都有一套新的坐标来匹配它在1000x1000网格中的实际位置和我的GPS数据的边界框坐标之外，并且ID会针对每个元组进行更新。

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

最后，我可以用

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

这给我，我试图输出运行这个表，我原来的数据表（gps_points）。问题是，它需要永久加载所有的小多边形，并且每次我想要生成具有不同数量的矩形的地图或运行具有不同底层坐标的数据集时，我必须删除sub_rects并生成重新加载它。

有没有更好的方法来做到这一点？我不需要图形输出。我只需要生成数据。不必从外部生成支持表（sub_rects）将会非常好，我怀疑实现相同目标的计算方法更少。我更愿意不必安装任何额外的软件。

ETA：按在评论请求，这里是查询计划（我家的机器，所以更小的数据集和其他表名，但同样的计划上）：

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)

ETA 2：

按在回答我修改了代码的建议建议有此：

（选择（选择 st_geomfromtext ROW_NUMBER（）OVER（由GEOM顺序）ID，GEOM（CONCAT（ '多边形（（' ，X || ''|| y，'，'，x + xstep ||' '|| y，'，'，x + xstep ||''|| y + ystep，'，'，x ||''|| （选择x，y FROM（选择 generate_series（xmin，xmin + xdelta，xstep）x）x，（（x，y）请选择 generate_series（ymin，ymin + ydelta，ystep）y）y）foo）bar）;

其中XMIN，YMIN，xdelta，ydelta，XSTEP和ystep全部由外部脚本计算，但可能只是如果包裹上面的函数调用中也被计算为一个Postgres功能的一部分。根据这个生成一个临时表并根据这个查询运行查询比我最初做的要快两个数量级。

来源

2017-07-24 Ged Mourrain

您能否发布您的查询的解释计划 –

两件事。首先在sql级别创建表（例如，从pg_admin开始）。

create table polygons as 
select st_geomfromtext(concat('Polygon((',x||' '||y,',',x||' 
'||y+0.2,',',x+0.4||' '||y+0.2,',',x+0.4||' '||y,',',x||' '||y,'))')) geom 
    FROM (select generate_series(0,199.9,0.2) x) x, 
     (select generate_series(0,199.9,0.4) y) y

创建索引

创建使用要旨（的geom）多边形索引;

然后使用您的查询或这一个。检查哪一个在你的情况下会更快

group by id;

来源

2017-07-24 21:52:49

谢谢！您创建的片段正是我所需要的 - 我修改了它以生成索引列，以便在后续查询中参考，但它快速且美观地解决了实际问题。至于使用st_dwithin进行匹配，它确实加快了速度，但仅以损失非常大部分匹配为代价，所以我将坚持在边界框上进行加工，这对于矩形和点应该足够了。 –

这很奇怪，用ST_DWithin（r.geom，p.geom，0）你会失去匹配。它应该采取它内部的每一点。我认为问题可能出现在两个多边形之间的边界上的点，因为它们会被计数两次。如果你可以显示你认为丢失的案例很少，我可以检查st_dwithin和&& –

之间有什么区别，但是我没有地方可以上传我运行过的数据集。太糟糕了，因为我希望对发生这种差异的原因有第二意见。我也对非Postgis版本进行了类似的处理，其中x，y为点，只是数字边界比较而不是Postgis矩形，虽然速度慢了3倍，但结果与边界框匹配相匹配，并且根本不匹配Dwithin 。 –

下面是从边界框生成网格的例子：

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;

来源

2017-07-25 19:49:46

谢谢！非常有趣的链接，尽管我认为Grzegorz下面提供的代码足以满足我的目的。要将其作为临时表运行，这是一个很好的建议，我肯定会实施它。 –

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

回答

相关问题