使用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從點數據生成熱量/密度圖

回答

相關問題