[postgis-users] Fast spatial join point to raster

[David Haynes]

Hello,

I'm following up on this blog post to see if there is a fast way to join
points to rasters.
http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html

I'm attempting to do a dasymetric mapping process, which is assigning areal
unit values to raster pixels. For this example, I have created two tables.
Table 1, mn_smokers has are geom (x,y) coordinates and spatial indices on
4326 and 26915. The resulting raster will be stored in mn_dasymetric and it
also has indices at 4326 and 26915. When running the explain query it seems
that index isn't really being used.


CREATE INDEX mn_smokers_gist_4326  ON mn_smokers  USING GIST (geom);
CREATE INDEX mn_smokers_gist_26915  ON mn_smokers  USING GIST
(ST_Transform(geom,26915));
CREATE INDEX mn_dasymetric_gist_4326 ON mn_dasymetric  USING GIST
(ST_Transform(ST_Envelope(rast),4326));
CREATE INDEX mn_dasymetric_gist_26915  ON mn_dasymetric  USING GIST
(ST_Transform(ST_Envelope(rast),26915));

Type: Nested Loop (Inner); ; Cost: 6.22 - 1594516.11
Type: Limit; ; Cost: 5.11 - 6.22
Type: Nested Loop (Inner); ; Cost: 5.11 - 471.95
Type: Seq Scan; Rel: minnesota_counties ; Cost: 0.00 - 34.09
Type: Bitmap Heap Scan; Rel: mn_smokers ; Cost: 5.11 - 437.50
Type: Bitmap Index Scan; Rel: mn_smokers_gist_4326 ; Cost: 0.00 - 5.10
Type: CTE Scan; Rel: points ; Cost: 0.00 - 4.00
Type: Materialize; ; Cost: 0.00 - 59.57
Type: Seq Scan; Rel: mn_dasymetric ; Cost: 0.00 - 47.38


with points as
(
select sp_id as point_id, smoker as point_value, ms.geom
from mn_smokers ms
inner join minnesota_counties mc on ST_Intersects(ms.geom, mc.geom)
)
, point_raster_agg as
(
select geom, count(point_id) as pixel_point_value
from points
group by geom
)
, output_rast as
(
select r.rid, pixel_point_value,  ST_SetValue(rast, ST_Transform(p.geom,
26915), pixel_point_value ) as setrast
,rast, geom
from point_raster_agg p
inner join mn_dasymetric r on ST_intersects(r.rast, ST_Transform(p.geom,
26915))
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210418/e99c385a/attachment.html>