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
0

This thread is a little old. I thought it would be good to share my method for solving the issue presented. I have GEOS 3.8.0 installed so the ST_Union is pretty snappy. I haven't tried testing it on large datasets.

I use a function to determine the bounding tiles, merge those into one raster, and then apply a mask with a buffer returning the masked raster in the raster's original format. The code is specific to my current application and needs to be made more general, e.g. pulling rasters from a different table (not just ned_13), taking the raster table's SRID, and selecting the correct UTM zone's SRID based on the geometry's location.

CREATE OR REPLACE FUNCTION
ST_GeoMask ( geom GEOMETRY, buf DOUBLE PRECISION DEFAULT 60, 
                    srid INTEGER DEFAULT 6345 )
RETURNS RASTER AS
$$
    DECLARE
        mask GEOMETRY := ST_Transform (
            ST_Buffer (
                ST_Transform( geom, srid ), 
                buf -- use rounded edge buffer the default
            ),
            4269 -- mask in NAD83 coordinates EPSG::4269
        );
        rasterout RASTER;
    BEGIN
        SELECT
            ST_Clip ( ST_Union ( rast ), mask, TRUE ) INTO rasterout
        FROM public.ned_13
        WHERE ST_Intersects( rast, geom );
        RETURN rasterout;
    END;  
$$
LANGUAGE plpgsql RETURNS NULL ON NULL INPUT;
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.