2

I am trying to ultimately calculate proportion of land cover types by polygon in PostGIS, so starting using st_valuecount to get counts of pixels with different land cover types. The data are pretty large (6" resolution land cover data and about 800,000 polygons, varying in size quite a bit).

The code below works as expected issue when applied to a subset of polygons, but when applied to the entire dataset, throws the following error:

SQL Error [XX000]: ERROR: invalid memory alloc request size 1144274944

Here's the code (adapted from PostGIS in Action, 2nd Ed.):

with
    cte as (
        select bbl,
            st_valuecount(
                st_clip(st_union(p.rast), geom_2263),1, true, ARRAY[0,1,2,3,4,5,6,7]
                ) as pv
            FROM
                base_rasters.landcover6in as p
                inner join
                results_scratch.polys
            on st_intersects(p.rast, geom_2263)
            group by bbl, geom_2263
        )
SELECT bbl, 
    (pv).value, sum((pv).count)
from cte
group by bbl, (pv).value
ORDER by bbl, (pv).value;

Of note, bringing the exact same datasets into QGIS 3.2.2 and running the 'Zonal Histogram' tool works without issue.

Any suggestions on what might be going on, or how to troubleshoot? (As much as anything, looking to have reproducible SQL code for doing similar tasks in the future.

Here's the software information (running on Win 7 x64 Professional): PostGIS Version

POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

PostgreSQL Version

PostgreSQL 10.3, compiled by Visual C++ build 1800, 64-bit

Lastly, I'll note I've previously posted about this on the PostGIS listserv but had not received responses, so posting here. I'll share any useful answers reciprocally.

Edit 1:

As noted in the comments, doing this operation without the union results in too many pixels being counted for some polygons (some of the pixels are of classes that do not even lie within the polygon). Thus, doing this without the union would likely work, but I'm observing issues in the clipping of the raster to the polygons, as follows:

Using the following code produces expected result for clipping to a single polygon:

select st_clip(st_union(p.rast), geom_2263)
 FROM base_rasters.landcover6in as p, 
      results_scratch.polys
      where st_intersects(p.rast, geom_2263) and bbl=5011720090
      group by geom_2263;

Here's the result - the yellow outline is the a focal polygon (bbl 5011720090); blue lines are borders of the raster tiles that intersect that polygon; the raster seen is the result of the clip, and the darker green filled polygons are the suite of polygons I'm pulling from. enter image description here

Here's the result when I remove the st_union:

select st_clip(p.rast, geom_2263)
 FROM base_rasters.landcover6in as p, 
      results_scratch.polys
        where st_intersects(p.rast, geom_2263) and bbl=5011720090;

Same color scheme as above; the black outside of the yellow/blue boundaries is 0 for nodata. You can see there are areas that are outside of the focal area (yellow boundaries) but within an associated tile that are maintained in this clip.

enter image description here

Any suggestions on what's going on here?

Edit 2

It appears the issue with st_clip across un-unioned raster tiles that I'm dealing with is a documented issue - https://trac.osgeo.org/postgis/ticket/3457 and documented in this stackexchange post.

Of course insights are welcome if I might be missing an alternative way of doing this in PostGIS.

Edit 3

Looks like there's now a bug fix, pending for the next PostGIS release, and with the ticket is now a way to redefine the function to fix for the interim.

  • 1
    You should probably give the OS version too, as I believe that could be relevant. – John Powell Aug 30 '18 at 14:06
  • 1
    What if you ST_Clip() before ST_Union()? St_Union must be running out of memory for some polygons. If you clip first, raster pieces to union will be smaller. – Pierre Racine Aug 31 '18 at 13:16
  • 1
    Also: Why do you have to ST_Union at all? Why not just aggregating the result of each clipped raster using 'bbl' and the polygon in the GROUP BY? Union is necessary only if you want the mean. For sum and count, it is not necessary. – Pierre Racine Aug 31 '18 at 13:17
  • Thanks @PierreRacine 1) If I re-org so ST_Union is outside the ST_Clip (st_valuecount(st_union(st_clip(p.rast, geom_2263))...), the result for larger polygons has extra pixels. (Though the sum isn't a multiple of the tile size I set on importing the raster, so its still clipping, but something's going awry). 2) If I remove ST_Union, I receive an error - 'p.rast needs to be in the group by or used in an aggregate function'. Putting p.rast into the 1st 'group by' I receive an error, that 'Some of the datatypes only support hashing, while others only support sorting.' Suggestions? – mtreg Aug 31 '18 at 14:23
  • 1) My two comments implied that there was no overlaps between the raster tiles in your raster coverage. If there are some, then some rasters parts, after the clipping, will overlap and produce extra pixels. 2) You have to remove the GROUP BY in the first subquery. You aggregate only in the last subquery. Understanding problems occuring in this type of query is very hard without having a visual look at the results. I generally look at them in OpenJump. Very easy... – Pierre Racine Sep 4 '18 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.