Sunday 26 November 2017

Error Quantifying Land Cover Classes by Polygon in PostGIS


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.



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?




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.



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.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...