[postgis-users] Incorrect Results for ST_ValueCount with large polygons

[Michael Treglia]

Hi again all,

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

Simpler/quicker:

  SELECT  bbl, (pvc).VALUE, SUM((pvc).COUNT) AS tot_pix
 FROM base_rasters.landcover6in
  INNER join  results_scratch.polys_test
  ON ST_Intersects(landcover6in.rast, polys_test.geom_2263),
    ST_ValueCount(ST_Clip(landcover6in.rast,polys_test.geom_2263),1, TRUE,
ARRAY[1, 2, 3, 4, 5, 6, 7]) AS pvc
  GROUP BY (pvc).VALUE, bbl
 ORDER BY bbl, (pvc).VALUE;



A longer, less efficient version that gives the same results (but unions
the clipped raster tiles):


 SELECT
   bbl, (value_count).value, SUM((value_count).count) AS count
  FROM
    (
    select
    bbl,
      rid,
      ST_ValueCount(
        ST_Union(ST_Clip(rast, geom_2263, TRUE)), 1, TRUE, ARRAY[1, 2, 3,
4, 5, 6, 7]
      ) value_count
    FROM
      (SELECT bbl, geom_2263 FROM results_scratch.polys_test) v,
      (SELECT rid, rast FROM base_rasters.landcover6in) r
    WHERE ST_Intersects(rast, geom_2263)
    GROUP BY bbl, rid, geom_2263
    ) i
  GROUP BY bbl, value
  ORDER BY bbl, value

On Thu, Jun 21, 2018 at 4:23 PM, Michael Treglia <mtreglia at gmail.com> wrote:

> Hi all,
>
> I'm basically trying to compute # of pixels of land cover type by polygon,
> for lots of polygons.
>
> I'm working with a query a query like this:
> SELECT  bbl, (pvc).VALUE, SUM((pvc).COUNT) AS tot_pix
>  FROM base_rasters.landcover6in
>   INNER join  results_scratch.polys_test
>   ON ST_Intersects(landcover6in.rast, polys_test.geom_2263),
>     ST_ValueCount(ST_Clip(landcover6in.rast,polys_test.geom_2263),1) AS
> pvc
>   GROUP BY (pvc).VALUE, bbl
>  ORDER BY (pvc).VALUE;
>
> 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)
>
> 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?
>
> Thanks in advance for any suggestions!
> Mike
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20180622/d0f9de3b/attachment.html>