Re: [postgis-users] Fast spatial join point to raster
Ok thanks I aligned my query to the functional indices and it is much faster. inner join mn_dasymetric r on ST_intersects(st_envelope(r.rast), ST_Transform(p.geom, 26915)) On Sun, Apr 18, 2021 at 1:21 PM Paul Ramsey wrote: > The indexing of rasters is (and remains a pain in the butt) (maybe > something I should look at shortly) because the index is bound to geometry > and almost always declraed functionally. You've built functional indexes, > > GIST (ST_Transform(ST_Envelope(rast),26915)) > > But you aren't calling them with queries that use the SAME functional > signature. The rule of functional indexes: > > * An index like this: CREATE INDEX foo_x ON (function(bar)) > * must be called like this: SELECT * FROM foo WHERE function(bar) = 'baz' > > Note that function(bar) appears in the CREATE and in the SELECT. > > P > > > On Apr 18, 2021, at 10:53 AM, David Haynes wrote: > > > > 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)) > > ___ > > postgis-users mailing list > > postgis-users@lists.osgeo.org > > https://lists.osgeo.org/mailman/listinfo/postgis-users > > ___ > postgis-users mailing list > postgis-users@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/postgis-users > ___ postgis-users mailing list postgis-users@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/postgis-users
Re: [postgis-users] Fast spatial join point to raster
The indexing of rasters is (and remains a pain in the butt) (maybe something I should look at shortly) because the index is bound to geometry and almost always declraed functionally. You've built functional indexes, GIST (ST_Transform(ST_Envelope(rast),26915)) But you aren't calling them with queries that use the SAME functional signature. The rule of functional indexes: * An index like this: CREATE INDEX foo_x ON (function(bar)) * must be called like this: SELECT * FROM foo WHERE function(bar) = 'baz' Note that function(bar) appears in the CREATE and in the SELECT. P > On Apr 18, 2021, at 10:53 AM, David Haynes wrote: > > 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)) > ___ > postgis-users mailing list > postgis-users@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/postgis-users ___ postgis-users mailing list postgis-users@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/postgis-users
[postgis-users] Fast spatial join point to raster
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)) ___ postgis-users mailing list postgis-users@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/postgis-users
