3

I have a table in postgis with a tiled raster and a table with polygons overlaying the raster. I want to get the number of pixels for each category (pixel value) that fall inside a specific polygon.

If I use the following query the result I am getting is wrong because, from what I understand, it hasn't unioned the tiles of the raster:

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

With the following query the result I am getting is the right one, but it takes longer than it should, since it unions all the tiles before clipping

WITH
  unionedraster AS (SELECT ST_Union(rast) As rast FROM rastertable)   
SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Clip(myraster.rast, myvector.geom)) As pvc
  FROM unionedraster As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

I thought the following would work, but the results are the same with the first query (like no union took part).

SELECT (pvc).value, SUM((pvc).count) As total
FROM (
  SELECT ST_ValueCount(ST_Union(ST_Clip(myraster.rast, myvector.geom))) As pvc
  FROM rastertable As myraster, vectortable As myvector
  WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
) As f
GROUP BY (pvc).value
ORDER BY (pvc).value;

What is the correct way to do this?


Upon further investigation, it seems that the problem lies in the clipping of the raster tiles. I managed to export the result of the clipping as an image using this query:

COPY (
WITH clipraster As (
    SELECT ST_Clip(myraster.rast, myvector.geom) As pvc
    FROM rastertable As myraster, vectortable As myvector
    WHERE ST_Intersects(myraster.rast, myvector.geom) AND myvector.gid = 6
)
SELECT encode(ST_AsPNG(
    clipraster.pvc
), 'hex') AS png FROM clipraster
) TO 'clip.hex';

and the resulting image was the following

enter image description here

The blue line is the polygon to which the raster tiles should be clipped.

However, it is obvious that some of the tiles where not clipped.

Any idea to what is the problem?

  • I would think that in your 2nd version, if you put a WHERE ST_Intersects(....) in your initial CTE query, so you only union those rasters that intersect your polygon and remove the ST_Intersects form the ST_ValueCount(ST_Clip(... part you will get a significant speed up. Currently, as you have unioned the entire raster, then it all passes the ST_Intersects in the inner select, which means you are doing excessive work in the clipping. – John Powell Mar 3 '18 at 9:42
  • Another possible speed up would be to clip first, and then union, in the initial query, see this post. However, this would seem to run against what you found in your 3rd query. – John Powell Mar 3 '18 at 10:01
  • John Powell, I have seen the post you mention and this is what I tried to do in my third query, which didn't work for some reason. As for your first comment, I understand that I should union only where the raster and vector intersect, but I can't find the correct syntax (which works) to do this. – gijo Mar 5 '18 at 7:40
  • OK, it was too long for a comment, so I put it as an answer instead. I believe that is what you are looking for, based on a test on my own data. – John Powell Mar 5 '18 at 11:19

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.