<div dir="ltr"><div dir="ltr"><div dir="ltr"><div>Hi All,<br></div><div><br></div><div>As one more follow-up on this thread... it sounds like the reason for the memory error I received is likely the large size of unioned rasters. For st_valuecount, it sounds like the st_union operation is not actually necessary, but without it, the st_clip operation seems to falter when a polygon overlaps multiple raster tiles (which is why I previously was seeing incorrect pixel counts for focal polygons without the st_union).</div><div><br></div><div>Some of this discovery on my part is documented in a StackExchange thread here - <a href="https://gis.stackexchange.com/questions/294482/error-quantifying-land-cover-classes-by-polygon-in-postgis">https://gis.stackexchange.com/questions/294482/error-quantifying-land-cover-classes-by-polygon-in-postgis</a> (and thanks to Pierre Racine for the helpful comments that led me to figure out what was going on).</div><div><br></div><div>And I believe the issue with st_clip is already documented here: <a href="https://trac.osgeo.org/postgis/ticket/3457">https://trac.osgeo.org/postgis/ticket/3457</a> (If I'm missing something, or if there's a way around this in PostGIS as-is, suggestions are welcome).</div><div><br></div><div>Cheers,</div><div>Mike<br></div><div><br></div></div></div></div><br><div class="gmail_quote"><div dir="ltr">On Thu, Jul 19, 2018 at 5:04 PM Michael Treglia <<a href="mailto:mtreglia@gmail.com">mtreglia@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi All,</div><div><br></div><div>As partial resolution to this, below are a couple of working examples of SQL, thanks to an example from Regina Obe's book. <br></div><div><br></div><div>These work on smaller subsets of the data I've used for testing, though unfortunately in running on the entire datasets, after running for hours, I receive an error: Invalid memory allocation request size 1109458944. Thoughts/suggestions welcome, but otherwise I hope these examples are useful to tothers.</div><div><br></div><div>Mike<br></div><div><br></div><div><br></div><div>--SQL CODE<br></div><div>with <br> cte as (<br> select bbl,<br> st_valuecount(<br> st_clip(st_union(p.rast), geom_2263),1, true, ARRAY[0,1,2,3,4,5,6,7]<br> ) as pv<br> FROM<br> base_rasters.landcover6in as p<br> inner join<br> results_scratch.polys_test<br> on st_intersects(p.rast, geom_2263)<br> group by bbl, geom_2263<br> )<br>SELECT bbl, <br> (pv).value, sum((pv).count)<br>from cte<br>group by bbl, (pv).value<br>ORDER by bbl, (pv).value;</div><div>--End SQL<br></div><div><br></div><div>And to get the results with land cover by polygon in a 'wide' format (e.g., amounts of different land cover in columns):</div><div><br></div><div>
--SQL CODE
<br></div><div>with <br> cte as (<br> select bbl,geom_2263,<br> st_valuecount(<br> st_clip(st_union(p.rast), geom_2263),1, true, ARRAY[0,1,2,3,4,5,6,7]<br> ) as pv<br> FROM<br> base_rasters.landcover6in as p<br> inner join<br> results_scratch.polys_test<br> on st_intersects(p.rast, geom_2263)<br> group by bbl, geom_2263<br> )<br>SELECT bbl, geom_2263, <br> (sum((pv).count) filter (where (pv).value = 0)::numeric/sum((pv).count))::numeric(5,4) as prop_na,<br> (sum((pv).count) filter (where (pv).value = 1)::numeric/sum((pv).count))::numeric(5,4) as prop_canopy,<br> (sum((pv).count) filter (where (pv).value = 2)::numeric/sum((pv).count))::numeric(5,4) as prop_grass,<br> (sum((pv).count) filter (where (pv).value = 3)::numeric/sum((pv).count))::numeric(5,4) as prop_bare,<br> (sum((pv).count) filter (where (pv).value = 4)::numeric/sum((pv).count))::numeric(5,4) as prop_water,<br> (sum((pv).count) filter (where (pv).value = 5)::numeric/sum((pv).count))::numeric(5,4) as prop_bldgs,<br> (sum((pv).count) filter (where (pv).value = 6)::numeric/sum((pv).count))::numeric(5,4) as prop_roads,<br> (sum((pv).count) filter (where (pv).value = 7)::numeric/sum((pv).count))::numeric(5,4) as prop_paved,<br> ((sum((pv).count) filter (where (pv).value = 0)) * 0.25)::numeric(12,2) as area_na,<br> ((sum((pv).count) filter (where (pv).value = 1)) * 0.25)::numeric(12,2) as area_canopy,<br> ((sum((pv).count) filter (where (pv).value = 2)) * 0.25)::numeric(12,2) as area_grass,<br> ((sum((pv).count) filter (where (pv).value = 3)) * 0.25)::numeric(12,2) as area_bare,<br> ((sum((pv).count) filter (where (pv).value = 4)) * 0.25)::numeric(12,2) as area_water,<br> ((sum((pv).count) filter (where (pv).value = 5)) * 0.25)::numeric(12,2) as area_bldgs,<br> ((sum((pv).count) filter (where (pv).value = 6)) * 0.25)::numeric(12,2) as area_roads,<br> ((sum((pv).count) filter (where (pv).value = 7)) * 0.25)::numeric(12,2) as area_paved<br>from cte<br>group by bbl, geom_2263<br>ORDER by bbl;</div><div>--End SQL<br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Jul 11, 2018 at 4:59 PM, Michael Treglia <span dir="ltr"><<a href="mailto:mtreglia@gmail.com" target="_blank">mtreglia@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi everyone,</div><div><br></div><div>Sorry for multiple emails on this, but wanted to ping this list again to see if there is anything obvious I'm missing or should check regarding this issue. I've reiterated context/the problem/sample sql code below, and if there's other info that would be helpful, please let me know.</div><div><br></div><div>Thanks! Sincerely,</div><div>Mike<br></div><div><br></div><div># To reiterate: <br></div><div>## The goal/issue:</div><div>I'm trying to calculate # of pixels from a categorical raster, by value, within distinct polygons of a vector layer. I have code, below, that seems to work accurately for small polygons, though I'm getting inaccuracies in large polygons - the addition of some pixels, including some from a value that aren't actually contained within a focal polygon. (I'm comparing to available calculations for the same overlay, along with results from QGIS Zonal Histogram tool). I'd like to get this working for this specific overlay, but then be able to apply the SQL code to comparable situations.</div><div><br></div><div>##Data:</div><div>The polygons are parcels, ranging quite widely in size;</div><div>The raster dataset is a 6" land cover dataset loaded as an in-db raster via raster2pgsql;<br></div><div>
<span class="m_-5658772553576828442m_-2378123914380777227gmail-pl-c">raster2pgsql -s 2263 -d -C -t 128x128 -M -I -l 4,16 "landcover_2010_nyc_05ft.img" base_rasters.landcover6in | psql -h localhost -U postgres -d dbname -v ON_ERROR_STOP=1</span>
<br></div><div><br></div><div>## System Info: <br></div><div>I'm running PostGIS 2.4.4 with PostgreSQL v 10.3 on Windows 7 x64</div><div><br></div><div>## Two sets of sample code that give the same result:</div><div>1) <br></div><span><div><span class="m_-5658772553576828442m_-2378123914380777227gmail-im">SELECT bbl, (pvc).VALUE, SUM((pvc).COUNT) AS tot_pix<br> FROM base_rasters.landcover6in<br> INNER join results_scratch.polys_test<br> ON ST_Intersects(landcover6in.rast, polys_test.geom_2263), <br></span> ST_ValueCount(ST_Clip(landcover6in.rast,polys_test.geom_2263),1, TRUE, ARRAY[1, 2, 3, 4, 5, 6, 7]) AS pvc<br> GROUP BY (pvc).VALUE, bbl<br> ORDER BY bbl, (pvc).VALUE;
<br></div><div><br></div></span><div>2) <br></div><span><div>
SELECT<br> bbl, (value_count).value, SUM((value_count).count) AS count<br> FROM<br> (<br> select<br> bbl,<br> rid,<br> ST_ValueCount(<br> ST_Union(ST_Clip(rast, geom_2263, TRUE)), 1, TRUE, ARRAY[1, 2, 3, 4, 5, 6, 7]<br> ) value_count<br> FROM<br> (SELECT bbl, geom_2263 FROM results_scratch.polys_test) v,<br> (SELECT rid, rast FROM base_rasters.landcover6in) r<br> WHERE ST_Intersects(rast, geom_2263)<br> GROUP BY bbl, rid, geom_2263<br> ) i<br> GROUP BY bbl, value<br> ORDER BY bbl, value
<br></div><div><br></div><div><br></div></span></div><div class="m_-5658772553576828442HOEnZb"><div class="m_-5658772553576828442h5"><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jun 22, 2018 at 5:12 PM, Michael Treglia <span dir="ltr"><<a href="mailto:mtreglia@gmail.com" target="_blank">mtreglia@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi again all,<br></div><div><br></div><div>Sorry for the multiple messages on this - realizing the query I posted was somewhat off from what I intended - below are examples that are getting at what I'm going for but exhibiting the issue I described (sometimes including pixels in counts aren't actually be within polygons [i.e., the categorical value isn't actually within the polygon based on visual inspection & this same analysis with other tools]).</div><div><br></div><div>Simpler/quicker:</div><div><br></div><div><span> SELECT bbl, (pvc).VALUE, SUM((pvc).COUNT) AS tot_pix<br> FROM base_rasters.landcover6in<br> INNER join results_scratch.polys_test<br> ON ST_Intersects(landcover6in.rast, polys_test.geom_2263), <br></span> ST_ValueCount(ST_Clip(landcover6in.rast,polys_test.geom_2263),1, TRUE, ARRAY[1, 2, 3, 4, 5, 6, 7]) AS pvc<br> GROUP BY (pvc).VALUE, bbl<br> ORDER BY bbl, (pvc).VALUE;<br> </div><div><br></div><div><br></div><div>A longer, less efficient version that gives the same results (but unions the clipped raster tiles):</div><div><br></div><div><br></div><div> SELECT<br> bbl, (value_count).value, SUM((value_count).count) AS count<br> FROM<br> (<br> select<br> bbl,<br> rid,<br> ST_ValueCount(<br> ST_Union(ST_Clip(rast, geom_2263, TRUE)), 1, TRUE, ARRAY[1, 2, 3, 4, 5, 6, 7]<br> ) value_count<br> FROM<br> (SELECT bbl, geom_2263 FROM results_scratch.polys_test) v,<br> (SELECT rid, rast FROM base_rasters.landcover6in) r<br> WHERE ST_Intersects(rast, geom_2263)<br> GROUP BY bbl, rid, geom_2263<br> ) i<br> GROUP BY bbl, value<br> ORDER BY bbl, value<br></div></div><div class="m_-5658772553576828442m_-2378123914380777227HOEnZb"><div class="m_-5658772553576828442m_-2378123914380777227h5"><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Jun 21, 2018 at 4:23 PM, Michael Treglia <span dir="ltr"><<a href="mailto:mtreglia@gmail.com" target="_blank">mtreglia@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Hi all,</div><div><br></div><div>I'm basically trying to compute # of pixels of land cover type by polygon, for lots of polygons. <br></div><div><br></div><div>I'm working with a query a query like this:<br></div><div></div><div>SELECT bbl, (pvc).VALUE, SUM((pvc).COUNT) AS tot_pix<br> FROM base_rasters.landcover6in<br> INNER join results_scratch.polys_test<br> ON ST_Intersects(landcover6in.rast, polys_test.geom_2263), <br> ST_ValueCount(ST_Clip(landcover6in.rast,polys_test.geom_2263),1) AS pvc<br> GROUP BY (pvc).VALUE, bbl<br> ORDER BY (pvc).VALUE;</div><div><br></div><div>For smaller polygons, spot-checking, it looks like I'm generally getting the right results. However, for larger polygons, not always. For example, though a polygon might not have roads, I'm getting non-0 counts for road pixels within a polygon (polygons have unique identifiers as 'bbl'). The total area of pixels counted ends up being larger than the polygon (larger than would be attributable to differences between polygons edges & pixels)<br></div><div><br></div><div>I'm thinking it's an issue when there are multiple raster tiles included... but not entirely sure. Any suggestions on whether there's something off in this query that might be causing an issue like I describe?</div><div><br></div><div>Thanks in advance for any suggestions!</div><div>Mike<br></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</div></div></blockquote></div><br></div>
</blockquote></div>