<div dir="ltr">Hello,<div><br></div><div>I'm following up on this blog post to see if there is a fast way to join points to rasters.</div><div><a href="http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html">http://blog.cleverelephant.ca/2018/09/postgis-external-storage.html</a><br></div><div><br></div><div>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. </div><div><br></div><div><br></div><div>CREATE INDEX mn_smokers_gist_4326  ON mn_smokers  USING GIST (geom);<br>CREATE INDEX mn_smokers_gist_26915  ON mn_smokers  USING GIST (ST_Transform(geom,26915));<br>CREATE INDEX mn_dasymetric_gist_4326 ON mn_dasymetric  USING GIST (ST_Transform(ST_Envelope(rast),4326));<br>CREATE INDEX mn_dasymetric_gist_26915  ON mn_dasymetric  USING GIST (ST_Transform(ST_Envelope(rast),26915));<br></div><div><br></div><div>Type: Nested Loop (Inner); ; Cost: 6.22 - 1594516.11<br>Type: Limit; ; Cost: 5.11 - 6.22<br>Type: Nested Loop (Inner); ; Cost: 5.11 - 471.95<br>Type: Seq Scan; Rel: minnesota_counties ; Cost: 0.00 - 34.09<br>Type: Bitmap Heap Scan; Rel: mn_smokers ; Cost: 5.11 - 437.50<br>Type: Bitmap Index Scan; Rel: mn_smokers_gist_4326 ; Cost: 0.00 - 5.10<br>Type: CTE Scan; Rel: points ; Cost: 0.00 - 4.00<br>Type: Materialize; ; Cost: 0.00 - 59.57<br>Type: Seq Scan; Rel: mn_dasymetric ; Cost: 0.00 - 47.38<br></div><div><br></div><div><br></div><div>with points as<br>(<br>select sp_id as point_id, smoker as point_value, ms.geom<br>from mn_smokers ms <br>inner join minnesota_counties mc on ST_Intersects(ms.geom, mc.geom)<br>)<br>, point_raster_agg as<br>(<br>select geom, count(point_id) as pixel_point_value<br>from points <br>group by geom<br>)<br>, output_rast as<br>(<br>select r.rid, pixel_point_value,  ST_SetValue(rast, ST_Transform(p.geom, 26915), pixel_point_value ) as setrast<br>,rast, geom<br>from point_raster_agg p<br>inner join mn_dasymetric r on ST_intersects(r.rast, ST_Transform(p.geom, 26915))<br></div></div>