<div dir="ltr"><div>Thanks, All,</div><div><br></div><div>The below solution ended up being the most efficient, by far, provided by Florian off-thread. It runs in about 0.5 seconds.  Regina - the solution you provided took about 30 seconds. Thanks to you both for offering solutions! (Regina - happy to share the test polygons if useful, just let me know).<br></div><div><br></div><div>I have a related thought/question that might be good for this list, though if it seems like this isn't the appropriate venue, please let me know so I keep in mind for the future (and I welcome a better alternative):<br></div><div></div><div>Something that's in mind is this type of operation is really a common need in work that I and colleagues engage in (e.g., area or % of landcover for various units is a common need), and the performance of PostGIS can be unparalleled in my experience when dealing with relatively enormous datasets. And because there is some work needed to figure out the logic in terms of what will work, it can be tough to get others on board and really see the utility (the learning curve can be steep).  <br></div><div><br></div><div>This brings up a question for me - Would it be feasible/reasonable to ultimately set up a function that is more generalizable for this type of work, so its more straightforward, especially for new users (perhaps more akin to the polygon/polygon solution I was working with, or the query I initially tried for the raster/polygon overlay).  To be clear, I know it would take time, effort, and maybe some contracted work to make it happen, but would that be a feasible thing to have done, at least as an add-in function (sensu postgisaddons [<a href="https://github.com/pedrogit/postgisaddons]">https://github.com/pedrogit/postgisaddons]</a>), if not something that's ultimately baked in?  (I've thought of a couple of similar examples, but figured this challenge might be a good starting point to explore)</div><div><br></div><div>Thanks again! Per usual, I'm impressed by the performance in the end, and the awesome community here. </div><div>Best,<br></div><div>Mike<br></div><div><br></div><div></div><div><br></div><div>-- Took about 0.5 seconds<br></div><div>with step_1 as (<br>    select rast<br>    from base_rasters.landcover6in_2017 a<br>             join scratch.smallsample b<br>                  on st_intersects(a.rast, b.geom_2263)),<br>     step_2 as<br>         (<br>             select (geomval).geom, (geomval).val<br>             from (<br>                      select distinct (st_dumpaspolygons(rast)) geomval<br>                      from step_1) foo<br>         ),<br>     step_3 as<br>         (select st_intersection(step_2.geom, b.geom_2263) geom, val<br>          from step_2<br>                   inner join scratch.smallsample b on st_intersects(step_2.geom, b.geom_2263))<br>select sum(st_area(geom)), val<br>from step_3<br>group by val;</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Feb 1, 2022 at 4:32 PM Regina Obe <<a href="mailto:lr@pcorp.us">lr@pcorp.us</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div lang="EN-US"><div class="gmail-m_4055448553123884488WordSection1"><p class="MsoNormal">Try<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal">select<br>(foo.geomval).val,<br>sum(st_area((geomval).geom)) as area_ft2<br>from (<br>    select (st_intersection( ST_Clip(a.rast, ST_Envelope(b.geom_2263))  , b.geom_2263)) as geomval<br>   from base_rasters.landcover6in_2017 a<br>   join scratch.smallsample b<br>   on st_intersects(a.rast, b.geom_2263)) as foo<br>group by val;<u></u><u></u></p><p class="MsoNormal"><span style="font-size:11pt;font-family:"Calibri",sans-serif;color:rgb(31,73,125)"><u></u> <u></u></span></p><p class="MsoNormal"><span style="font-size:11pt;font-family:"Calibri",sans-serif;color:rgb(31,73,125)"><u></u> <u></u></span></p><div style="border-color:currentcolor currentcolor currentcolor blue;border-style:none none none solid;border-width:medium medium medium 1.5pt;padding:0in 0in 0in 4pt"><div><div style="border-color:rgb(225,225,225) currentcolor currentcolor;border-style:solid none none;border-width:1pt medium medium;padding:3pt 0in 0in"><p class="MsoNormal"><b><span style="font-size:11pt;font-family:"Calibri",sans-serif">From:</span></b><span style="font-size:11pt;font-family:"Calibri",sans-serif"> postgis-users [mailto:<a href="mailto:postgis-users-bounces@lists.osgeo.org" target="_blank">postgis-users-bounces@lists.osgeo.org</a>] <b>On Behalf Of </b>Mike Treglia<br><b>Sent:</b> Monday, January 31, 2022 10:57 AM<br><b>To:</b> Florian Nadler <<a href="mailto:florian.nadler@cybertec.at" target="_blank">florian.nadler@cybertec.at</a>><br><b>Cc:</b> PostGIS Users Discussion <<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>><br><b>Subject:</b> Re: [postgis-users] ST_Intersection - Slow with raster, fast with polygonized raster?<u></u><u></u></span></p></div></div><p class="MsoNormal"><u></u> <u></u></p><div><div><p class="MsoNormal">Thanks, Florian,<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Both of those alternatives take ~1 min 20 sec.<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Below are the execution plans for the queries I previously shared. Let me know if there's anything else that would be helpful.<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Best,<u></u><u></u></p></div><div><p class="MsoNormal">Mike<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">------Raster/Vector Overlay<u></u><u></u></p></div><div><p class="MsoNormal">explain select<br>(foo.geomval).val,<br>sum(st_area((geomval).geom)) as area_ft2<br>from (<br>    select (st_intersection(a.rast, b.geom_2263)) as geomval<br>   from base_rasters.landcover6in_2017 a<br>   join scratch.smallsample b<br>   on st_intersects(a.rast, b.geom_2263)) as foo<br>group by val;<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">HashAggregate  (cost=3904551.46..3904553.46 rows=200 width=16)<br>  Group Key: (foo.geomval).val<br>  ->  Subquery Scan on foo  (cost=0.41..456301.46 rows=26525000 width=40)<br>        ->  ProjectSet  (cost=0.41..191051.46 rows=26525000 width=32)<br>              ->  Nested Loop  (cost=0.41..51662.59 rows=26525 width=803)<br>                    ->  Seq Scan on smallsample b  (cost=0.00..1.28 rows=28 width=32)<br>                    ->  Index Scan using landcover6in_2017_st_convexhull_idx on landcover6in_2017 a  (cost=0.41..1844.10 rows=95 width=771)<br>                          Index Cond: ((rast)::geometry && b.geom_2263)<br>                          Filter: _st_intersects(b.geom_2263, rast, NULL::integer)<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">--Polygonized Raster/Vector Overlay<u></u><u></u></p></div><div><p class="MsoNormal">explain select lc_2017_value,<br>sum(st_area(st_intersection(a.geom_2263, b.geom_2263))) as area_ft2<br>from base_rasters.landcover6in_2017_polygons a<br>join scratch.smallsample b<br>on st_intersects(a.geom_2263, b.geom_2263)<br>group by lc_2017_value;<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">HashAggregate  (cost=12518962.15..12518962.23 rows=8 width=16)<br>  Group Key: a.lc_2017_value<br>  ->  Nested Loop  (cost=0.42..1295125.12 rows=446631 width=1388)<br>        ->  Seq Scan on smallsample b  (cost=0.00..1.28 rows=28 width=32)<br>        ->  Index Scan using landcover6in_2017_polygons_geom_idx on landcover6in_2017_polygons a  (cost=0.42..46238.47 rows=1595 width=1356)<br>              Index Cond: (geom_2263 && b.geom_2263)<br>              Filter: st_intersects(geom_2263, b.geom_2263)<u></u><u></u></p></div></div><p class="MsoNormal"><u></u> <u></u></p><div><div><p class="MsoNormal">On Mon, Jan 31, 2022 at 2:42 AM Florian Nadler <<a href="mailto:florian.nadler@cybertec.at" target="_blank">florian.nadler@cybertec.at</a>> wrote:<u></u><u></u></p></div><blockquote style="border-color:currentcolor currentcolor currentcolor rgb(204,204,204);border-style:none none none solid;border-width:medium medium medium 1pt;padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in"><div><p>Hi,<u></u><u></u></p><p>can you provide execution plans for both queries?<u></u><u></u></p><p>In the meantime, please try the following variants:<u></u><u></u></p><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">select<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(8,8,8)">(</span><span style="color:black">foo</span><span style="color:rgb(8,8,8)">.</span><span style="color:black">geomval</span><span style="color:rgb(8,8,8)">).val,<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><i><span style="color:rgb(8,8,8)">sum</span></i><span style="color:rgb(8,8,8)">(</span><i><span style="color:rgb(0,98,122)">st_area</span></i><span style="color:rgb(8,8,8)">((</span><span style="color:black">geomval</span><span style="color:rgb(8,8,8)">).geom)) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">area_ft2<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">from </span><span style="color:rgb(8,8,8)">(<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(8,8,8)">    </span><span style="color:rgb(0,51,179)">select </span><span style="color:rgb(8,8,8)">(</span><i><span style="color:rgb(0,98,122)">st_intersection</span></i><span style="color:rgb(8,8,8)">(</span><span style="color:black">a</span><span style="color:rgb(8,8,8)">.<i>rast</i>, </span><span style="color:black">b</span><span style="color:rgb(8,8,8)">.<i>geom_2263</i>)) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">geomval<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">from </span><span style="color:rgb(8,8,8)">base_rasters.landcover6in_2017 </span><span style="color:black">a<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">join </span><span style="color:rgb(8,8,8)">scratch.smallsample </span><span style="color:black">b<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">on </span><i><span style="color:rgb(0,98,122)">st_intersects</span></i><span style="color:rgb(8,8,8)">(</span><span style="color:black">a</span><span style="color:rgb(8,8,8)">.<i>rast</i>::</span><span style="color:rgb(0,51,179)">geometry</span><span style="color:rgb(8,8,8)">, </span><span style="color:black">b</span><span style="color:rgb(8,8,8)">.<i>geom_2263</i>)) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">foo<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">group by </span><span style="color:rgb(8,8,8)">val;<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">select<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(8,8,8)">(</span><span style="color:black">foo</span><span style="color:rgb(8,8,8)">.</span><span style="color:black">geomval</span><span style="color:rgb(8,8,8)">).val,<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><i><span style="color:rgb(8,8,8)">sum</span></i><span style="color:rgb(8,8,8)">(</span><i><span style="color:rgb(0,98,122)">st_area</span></i><span style="color:rgb(8,8,8)">((</span><span style="color:black">geomval</span><span style="color:rgb(8,8,8)">).geom)) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">area_ft2<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">from </span><span style="color:rgb(8,8,8)">(<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(8,8,8)">    </span><span style="color:rgb(0,51,179)">select </span><span style="color:rgb(8,8,8)">(</span><i><span style="color:rgb(0,98,122)">st_intersection</span></i><span style="color:rgb(8,8,8)">(</span><span style="color:black">a</span><span style="color:rgb(8,8,8)">.<i>rast</i>, </span><span style="color:black">b</span><span style="color:rgb(8,8,8)">.<i>geom_2263</i>)) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">geomval<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">from </span><span style="color:rgb(8,8,8)">base_rasters.landcover6in_2017 </span><span style="color:black">a<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">join </span><span style="color:rgb(8,8,8)">scratch.smallsample </span><span style="color:black">b<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:black">   </span><span style="color:rgb(0,51,179)">on </span><i><span style="color:rgb(0,98,122)">st_intersects</span></i><span style="color:rgb(8,8,8)">(</span><span style="color:black">a</span><span style="color:rgb(8,8,8)">.<i>rast</i>, </span><span style="color:black">b</span><span style="color:rgb(8,8,8)">.<i>geom_2263</i>) </span><span style="color:rgb(0,51,179)">where </span><span style="color:black">a</span><span style="color:rgb(8,8,8)">.<i>rast </i></span><span style="color:black">&& b</span><span style="color:rgb(8,8,8)">.<i>geom_2263</i>) </span><span style="color:rgb(0,51,179)">as </span><span style="color:black">foo<u></u><u></u></span></pre><pre style="background:white none repeat scroll 0% 0%"><span style="color:rgb(0,51,179)">group by </span><span style="color:rgb(8,8,8)">val<u></u><u></u></span></pre><div><p class="MsoNormal">Am 31.01.2022 um 04:32 schrieb Mike Treglia:<u></u><u></u></p></div><blockquote style="margin-top:5pt;margin-bottom:5pt"><div><div><p class="MsoNormal">Hi all,<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">I was exploring raster/vector overlays, basically to calculate the area of each of several landcover classes (from raster data) within an area or sets of areas (delineated in polygon/multipolygon data). In the past, I had ended up polygonizing the raster, making sure the polygonized layer had a spatial index, and developing a query from there. <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">But I wanted to come back to the problem and see if I could do it in what seemed like the more elegant way, using the raster directly without an intermediate step of polygonizing it.  While I'm able to do it, at least on a small sample area, it's at a huge performance cost. <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">I'm wondering if I'm just missing something I need to do to make the raster-based version operate more efficiently. (For what it's worth, rasters like this, are also slow to load into QGIS and such [though once loaded, are fairly quick to render], so I wonder if that's related at all and fixable.)<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">My test case is in New York City - I'm working with a six-inch resolution landcover dataset for the entirety of NYC (<a href="https://data.cityofnewyork.us/Environment/Land-Cover-Raster-Data-2017-6in-Resolution/he6d-2qns" target="_blank">https://data.cityofnewyork.us/Environment/Land-Cover-Raster-Data-2017-6in-Resolution/he6d-2qns</a>), which I'm overlaying with a small subset of parcels in a small geographic area. In this example, below, I'm just adding up the area for each landcover class across the focal parcels for simplicity.  When using the polygonized landcover data, it takes about 170ms. When using the raster landcover data, it takes about 1 min 20 seconds. <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Below are the queries, as well as the information about the raster table in PostGIS, and the raster2pgsql code I used to load the layer into the database. Curious to hear any thoughts from folks on the list. Of course if there's anything else I can share that would be helpful, let me know<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Thanks so much! Best,<u></u><u></u></p></div><div><p class="MsoNormal">Mike<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Here are the queries I'm running:<u></u><u></u></p></div><div><p class="MsoNormal">1) Original Raster over Polygons- ~1 min 20 sec<u></u><u></u></p></div><div><p class="MsoNormal">---------<u></u><u></u></p></div><div><p class="MsoNormal">select <br>(foo.geomval).val, <br>sum(st_area((geomval).geom)) as area_ft2 <br>from (<br>    select (st_intersection(a.rast, b.geom_2263)) as geomval<br>   from base_rasters.landcover6in_2017 a <br>   join scratch.smallsample b<br>   on st_intersects(a.rast, b.geom_2263)) as foo<br>group by val;<br><br>2) Polygonized Raster over Polygons ~0.170 ms<u></u><u></u></p></div><div><p class="MsoNormal">---------<u></u><u></u></p></div><div><p class="MsoNormal">select lc_2017_value, <br>sum(st_area(st_intersection(a.geom_2263, b.geom_2263))) as area_ft2<br>from base_rasters.landcover6in_2017_polygons a<br>join scratch.smallsample b<br>on st_intersects(a.geom_2263, b.geom_2263)<br>group by lc_2017_value;<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Here's the information for the raster table. (I had added a second index in hopes it might help but no, it didn't).<u></u><u></u></p></div><div><p class="MsoNormal">---------<u></u><u></u></p></div><div><p class="MsoNormal">                                Table "base_rasters.landcover6in_2017"<br> Column |  Type   | Collation | Nullable |                           Default<br>--------+---------+-----------+----------+-------------------------------------------------------------<br> rid    | integer |           | not null | nextval('base_rasters.landcover6in_2017_rid_seq'::regclass)<br> rast   | raster  |           |          |<br>Indexes:<br>    "landcover6in_2017_pkey" PRIMARY KEY, btree (rid)<br>    "landcover6in_2017_envidx" gist (st_envelope(rast))<br>    "landcover6in_2017_st_convexhull_idx" gist (st_convexhull(rast))<br>Check constraints:<br>    "enforce_height_rast" CHECK (st_height(rast) = ANY (ARRAY[128, 46]))<br>    "enforce_max_extent_rast" CHECK (st_envelope(rast) @ '0103000020D70800000100000005000000C3F528DC3DD72B41CCCCCCCCB46BFC40C3F528DC3DD72B413333333349B31041E27A14EE9E4930413333333349B31041E2<br>7A14EE9E493041CCCCCCCCB46BFC40C3F528DC3DD72B41CCCCCCCCB46BFC40'::geometry) NOT VALID<br>    "enforce_nodata_values_rast" CHECK (_raster_constraint_nodata_values(rast) = '{0.0000000000}'::numeric[])<br>    "enforce_num_bands_rast" CHECK (st_numbands(rast) = 1)<br>    "enforce_out_db_rast" CHECK (_raster_constraint_out_db(rast) = '{f}'::boolean[])<br>    "enforce_pixel_types_rast" CHECK (_raster_constraint_pixel_types(rast) = '{8BUI}'::text[])<br>    "enforce_same_alignment_rast" CHECK (st_samealignment(rast, '0100000000000000000000E03F000000000000E0BFC3F528DCBDB92E413333333349B3104100000000000000000000000000000000D708000001000100'::<br>raster))<br>    "enforce_scalex_rast" CHECK (round(st_scalex(rast)::numeric, 10) = round(0.5, 10))<br>    "enforce_scaley_rast" CHECK (round(st_scaley(rast)::numeric, 10) = round(- 0.5, 10))<br>    "enforce_srid_rast" CHECK (st_srid(rast) = 2263)<br>    "enforce_width_rast" CHECK (st_width(rast) = 128)<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Here's the raster2pgsql code I used to load the raster in the first place:<u></u><u></u></p></div><div><p class="MsoNormal">-------<u></u><u></u></p></div><div><p class="MsoNormal">raster2pgsql -s 2263 -d -C -t 128x128 -M -I -l 4,16 "NYC_2017_LiDAR_LandCover.tif" base_rasters.landcover6in_2017 | psql -h localhost -U postgres -d nycdb -v ON_ERROR_STOP=1<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Lastly - here's my PostGIS info (Running on Windows 10, 64 bit): <u></u><u></u></p></div><div><p class="MsoNormal">POSTGIS="3.1.1 3.1.1" [EXTENSION] PGSQL="130" GEOS="3.9.0-CAPI-1.14.1" SFCGAL="1.3.8" PROJ="7.1.1" GDAL="GDAL 3.2.1, released 2020/12/29" LIBXML="2.9.9" LIBJSON="0.12" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div></div><p class="MsoNormal"><u></u> <u></u></p><pre>_______________________________________________<u></u><u></u></pre><pre>postgis-users mailing list<u></u><u></u></pre><pre><a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><u></u><u></u></pre><pre><a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><u></u><u></u></pre></blockquote><pre>-- <u></u><u></u></pre><pre>CYBERTEC PostgreSQL International GmbH<u></u><u></u></pre><pre>Römerstraße 19, A-2752 Wöllersdorf<u></u><u></u></pre><pre>Web: <a href="https://www.cybertec-postgresql.com" target="_blank">https://www.cybertec-postgresql.com</a><u></u><u></u></pre></div></blockquote></div></div></div></div>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">https://lists.osgeo.org/mailman/listinfo/postgis-users</a><br>
</blockquote></div>