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 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.
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.
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.
Any suggestions on what's going on here?
Of course insights are welcome if I might be missing an alternative way of doing this in PostGIS.
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.