[postgis-users] ST_Intersection - Slow with raster, fast with polygonized raster?

Mike Treglia mtreglia at gmail.com
Sun Jan 30 19:32:00 PST 2022

Hi all,

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.

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.

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.)

My test case is in New York City - I'm working with a six-inch resolution
landcover dataset for the entirety of NYC (
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.

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

Thanks so much! Best,

Here are the queries I'm running:
1) Original Raster over Polygons- ~1 min 20 sec
sum(st_area((geomval).geom)) as area_ft2
from (
    select (st_intersection(a.rast, b.geom_2263)) as geomval
   from base_rasters.landcover6in_2017 a
   join scratch.smallsample b
   on st_intersects(a.rast, b.geom_2263)) as foo
group by val;

2) Polygonized Raster over Polygons ~0.170 ms
select lc_2017_value,
sum(st_area(st_intersection(a.geom_2263, b.geom_2263))) as area_ft2
from base_rasters.landcover6in_2017_polygons a
join scratch.smallsample b
on st_intersects(a.geom_2263, b.geom_2263)
group by lc_2017_value;

Here's the information for the raster table. (I had added a second index in
hopes it might help but no, it didn't).
                                Table "base_rasters.landcover6in_2017"
 Column |  Type   | Collation | Nullable |                           Default
 rid    | integer |           | not null |
 rast   | raster  |           |          |
    "landcover6in_2017_pkey" PRIMARY KEY, btree (rid)
    "landcover6in_2017_envidx" gist (st_envelope(rast))
    "landcover6in_2017_st_convexhull_idx" gist (st_convexhull(rast))
Check constraints:
    "enforce_height_rast" CHECK (st_height(rast) = ANY (ARRAY[128, 46]))
    "enforce_max_extent_rast" CHECK (st_envelope(rast) @
    "enforce_nodata_values_rast" CHECK
(_raster_constraint_nodata_values(rast) = '{0.0000000000}'::numeric[])
    "enforce_num_bands_rast" CHECK (st_numbands(rast) = 1)
    "enforce_out_db_rast" CHECK (_raster_constraint_out_db(rast) =
    "enforce_pixel_types_rast" CHECK (_raster_constraint_pixel_types(rast)
= '{8BUI}'::text[])
    "enforce_same_alignment_rast" CHECK (st_samealignment(rast,
    "enforce_scalex_rast" CHECK (round(st_scalex(rast)::numeric, 10) =
round(0.5, 10))
    "enforce_scaley_rast" CHECK (round(st_scaley(rast)::numeric, 10) =
round(- 0.5, 10))
    "enforce_srid_rast" CHECK (st_srid(rast) = 2263)
    "enforce_width_rast" CHECK (st_width(rast) = 128)

Here's the raster2pgsql code I used to load the raster in the first place:
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

Lastly - here's my PostGIS info (Running on Windows 10, 64 bit):
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)"
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20220130/1da8d987/attachment.html>

More information about the postgis-users mailing list