Intersecting a raster with a polygon using PostGIS - artefact error

Question

I'm using PostGIS2.0 to do some raster/polygon intersections. I'm having difficulty understanding which operation I should use, and what the quickest way of performing this is. My problem is as follows:

• I have a polygon and a raster
• I want to find all the pixels that fall within the polygon, and get the sum of the pixel value
• And (updated problem): I'm getting massive values for some pixels that do not exist in the original raster when I perform the query

I'm having difficulty understanding whether I should use ST_Intersects() or ST_Intersection(). I also don't know what the best approach for summing my pixels is. Here is the first approach that I have tried (#1):

SELECT 
    r.rast 
FROM
    raster as r, 
    polygon as p
WHERE
    ST_Intersects(r.rast, p.geom)

This returns a list of rast values, which I'm not sure what to do with. I tried calculating the summary statistics using ST_SummaryStats() but I'm not certain if this is the weighted-sum of all the pixels that lie within the polygon.

SELECT  
        (result).count,
        (result).sum    
FROM (  
         SELECT 
            ST_SummaryStats(r.rast) As result
         FROM
            raster As r, 
            polygon As p
         WHERE
            ST_Intersects(r.rast, p.geom)    
    ) As tmp

The other approach which I have tried (#2) uses ST_Intersection():

 SELECT
        (gv).geom,         
        (gv).val
 FROM 
 (
    SELECT 
        ST_Intersection(r.rast, p.geom) AS gv
    FROM 
        raster as r, 
        polygon as p           
    WHERE 
        ST_Intersects(r.rast, p.geom)

      ) as foo;

This returns a list of geometries that I analyze further, but I assume this is less efficient.

I'm unclear on which is the quickest order of operation also. Should I always choose raster, polygon or polygon, raster, or convert the polygon into a raster so that it is raster, raster ?

EDIT: I updated approach #2 with some details from R.K.'s link.

Using approach #2, I have noticed the following error in the results which is part of the reason why I did not understand the output. Here is the image of my original raster, and an outline of the polygon that is being used to intersect it, overlaid on top:

And here is the result of the intersection using PostGIS:

The problem with the result is that there are values of 21474836 being returned, which are not in the original raster. I've no idea why this is occurring. I suspect it is related to small numbers somewhere (dividing by almost 0), but it returns the wrong result.

Answer 1

I found a tutorial on intersecting vector polygons with a large raster coverage using PostGIS WKT Raster. Might have the answers you're looking for.

The tutorial used two datasets, a point shape file that was buffered to produce polygons and a series of 13 SRTM elevation rasters. There was a lot of steps in between but the query used to intersect the raster and vector looked like this:

 CREATE TABLE caribou_srtm_inter AS
 SELECT id, 
        (gv).geom AS the_geom, 
        (gv).val
 FROM (SELECT id, 
              ST_Intersection(rast, the_geom) AS gv
       FROM srtm_tiled,
            cariboupoint_buffers_wgs
       WHERE ST_Intersects(rast, the_geom)
      ) foo;

The values were then summarized using the following:

 CREATE TABLE result01 AS
 SELECT id, 
        sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 
        sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
 FROM caribou_srtm_inter
 GROUP BY id
 ORDER BY id;

I don't really know enough PostGIS to explain this but it sure sounds like what you were trying to accomplish. The tutorial should shed light on the intermediate steps. Good luck :)

Answer 2

With regard to point 2 in the original question - several of the postgis 2.0 development releases used a version of the GDAL library that cast floats to ints. If your raster has float values in it, and you were using a version of GDAL lower than 1.9.0, or a version of PostGIS 2.0 prerelease that didn't properly call GDALFPolygonize(), then you might be encountering this bug. Tickets in the PostGIS and GDAL bug trackers were filed and closed. This bug was active around the time of the original question.

In terms of performance, you'll find that using ST_Intersects(raster, geom) is much much faster than using ST_Intersects(geom, raster). The first version rasterizes the geometry and does a raster-space intersection. The second version vectorizes the geometry and does a vector-space intersection, which can be far more expensive process.

Answer 3

I was also having strange issues using ST_SummaryStats with ST_Clip. Querying the data differently told me the min value of my raster was 32, and then max 300, yet ST_SummaryStats was returning -32700 for the pixel values in my target polygon.

I ended up hacking around the issue as so:

WITH first AS (
   SELECT id, (ST_Intersection(geom, rast)).val
   FROM raster_table
   INNER JOIN vector_table ON ST_Intersects(rast, geom)
)
SELECT id, COUNT(val), SUM(val), AVG(val), stddev(val), MIN(val), MAX(val)
FROM first
GROUP BY id