<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hi,</p>
<p>can you provide execution plans for both queries?</p>
<p>In the meantime, please try the following variants:</p>
<pre style="background-color:#ffffff;color:#080808;font-family:'JetBrains Mono',monospace;font-size:9,0pt;"><span style="color:#0033b3;">select
</span>(<span style="color:#000000;">foo</span>.<span style="color:#000000;">geomval</span>).val,
<span style="font-style:italic;">sum</span>(<span style="color:#00627a;font-style:italic;">st_area</span>((<span style="color:#000000;">geomval</span>).geom)) <span style="color:#0033b3;">as </span><span style="color:#000000;">area_ft2
</span><span style="color:#0033b3;">from </span>(
<span style="color:#0033b3;">select </span>(<span style="color:#00627a;font-style:italic;">st_intersection</span>(<span style="color:#000000;">a</span>.<span style="font-style:italic;">rast</span>, <span style="color:#000000;">b</span>.<span style="font-style:italic;">geom_2263</span>)) <span style="color:#0033b3;">as </span><span style="color:#000000;">geomval
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">from </span>base_rasters.landcover6in_2017 <span style="color:#000000;">a
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">join </span>scratch.smallsample <span style="color:#000000;">b
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">on </span><span style="color:#00627a;font-style:italic;">st_intersects</span>(<span style="color:#000000;">a</span>.<span style="font-style:italic;">rast</span>::<span style="color:#0033b3;">geometry</span>, <span style="color:#000000;">b</span>.<span style="font-style:italic;">geom_2263</span>)) <span style="color:#0033b3;">as </span><span style="color:#000000;">foo
</span><span style="color:#0033b3;">group by </span>val;</pre>
<pre style="background-color:#ffffff;color:#080808;font-family:'JetBrains Mono',monospace;font-size:9,0pt;">
<span style="color:#0033b3;">select
</span>(<span style="color:#000000;">foo</span>.<span style="color:#000000;">geomval</span>).val,
<span style="font-style:italic;">sum</span>(<span style="color:#00627a;font-style:italic;">st_area</span>((<span style="color:#000000;">geomval</span>).geom)) <span style="color:#0033b3;">as </span><span style="color:#000000;">area_ft2
</span><span style="color:#0033b3;">from </span>(
<span style="color:#0033b3;">select </span>(<span style="color:#00627a;font-style:italic;">st_intersection</span>(<span style="color:#000000;">a</span>.<span style="font-style:italic;">rast</span>, <span style="color:#000000;">b</span>.<span style="font-style:italic;">geom_2263</span>)) <span style="color:#0033b3;">as </span><span style="color:#000000;">geomval
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">from </span>base_rasters.landcover6in_2017 <span style="color:#000000;">a
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">join </span>scratch.smallsample <span style="color:#000000;">b
</span><span style="color:#000000;"> </span><span style="color:#0033b3;">on </span><span style="color:#00627a;font-style:italic;">st_intersects</span>(<span style="color:#000000;">a</span>.<span style="font-style:italic;">rast</span>, <span style="color:#000000;">b</span>.<span style="font-style:italic;">geom_2263</span>) <span style="color:#0033b3;">where </span><span style="color:#000000;">a</span>.<span style="font-style:italic;">rast </span><span style="color:#000000;">&& b</span>.<span style="font-style:italic;">geom_2263</span>) <span style="color:#0033b3;">as </span><span style="color:#000000;">foo
</span><span style="color:#0033b3;">group by </span>val</pre>
<div class="moz-cite-prefix">Am 31.01.2022 um 04:32 schrieb Mike
Treglia:<br>
</div>
<blockquote type="cite"
cite="mid:CAPKp32v4NYoPk_Q4nWV+9UcA8LJQExVq7ZXROP6L709gzcFjZA@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<div dir="ltr">
<div>Hi all,</div>
<div><br>
</div>
<div>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. <br>
</div>
<div><br>
</div>
<div>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. <br>
</div>
<div><br>
</div>
<div>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.)<br>
</div>
<div><br>
</div>
<div>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"
moz-do-not-send="true" class="moz-txt-link-freetext">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. <br>
</div>
<div><br>
</div>
<div>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<br>
</div>
<div><br>
</div>
<div>Thanks so much! Best,</div>
<div>Mike<br>
</div>
<div><br>
</div>
<div>Here are the queries I'm running:</div>
<div>1) Original Raster over Polygons- ~1 min 20 sec<br>
</div>
<div>---------<br>
</div>
<div>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<br>
</div>
<div>---------<br>
</div>
<div>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;</div>
<div><br>
</div>
<div><br>
</div>
<div>Here's the information for the raster table. (I had added a
second index in hopes it might help but no, it didn't).</div>
<div>---------</div>
<div> 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)</div>
<div><br>
</div>
<div>Here's the raster2pgsql code I used to load the raster in
the first place:</div>
<div>-------<br>
</div>
<div>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</div>
<div><br>
</div>
<div>Lastly - here's my PostGIS info (Running on Windows 10, 64
bit): <br>
</div>
<div>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</div>
<div><br>
</div>
</div>
<br>
<fieldset class="moz-mime-attachment-header"></fieldset>
<pre class="moz-quote-pre" wrap="">_______________________________________________
postgis-users mailing list
<a class="moz-txt-link-abbreviated" href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a>
<a class="moz-txt-link-freetext" href="https://lists.osgeo.org/mailman/listinfo/postgis-users">https://lists.osgeo.org/mailman/listinfo/postgis-users</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
CYBERTEC PostgreSQL International GmbH
Römerstraße 19, A-2752 Wöllersdorf
Web: <a class="moz-txt-link-freetext" href="https://www.cybertec-postgresql.com">https://www.cybertec-postgresql.com</a></pre>
</body>
</html>