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):

    raster as r, 
    polygon as p
    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.

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

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

        ST_Intersection(r.rast, p.geom) AS gv
        raster as r, 
        polygon as p           
        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:

enter image description here

And here is the result of the intersection using PostGIS:

enter image description here

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.

  • Regarding point two, you want to get the sum of the values of the pixels that intersect the polygon? – R.K. Feb 6 '12 at 15:04
  • Yes, I have used ST_SummaryStats() for #1, but am not sure how to do it for #2. – djq Feb 6 '12 at 15:09
  • Posted a link to a reference. I hope that helps. – R.K. Feb 6 '12 at 15:28
  • 2
    The maximum value of the scale in your map is the maximum of a 32-bit signed integer. I do not know what that means in your case, but it could have to do with NA values. The range in your query might have null values that are not handled properly. en.wikipedia.org/wiki/2147483647#2147483647_in_computing – yellowcap Feb 7 '12 at 14:59
  • 6
    FWIW, 21474836 is not the maximum value of a 32 bit signed int. However, 2^31-1 = 2147483647 is the max, and notice that 21474836 = 2147483647 / 100 (integer division). This hints that internally some NA's are generated (perhaps along border cells), they are represented as 2^31-1, and then the code "forgets" these are NA and (perhaps in a resampling process?) it divides them by 100. – whuber Feb 7 '12 at 16:20

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, 
              ST_Intersection(rast, the_geom) AS gv
       FROM srtm_tiled,
       WHERE ST_Intersects(rast, the_geom)
      ) foo;

The values were then summarized using the following:

 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

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 :)

  • Thanks @R.K. I did read through that tutorial. I think my calculation is more basic, yet I'm still figuring out this basic step! – djq Feb 6 '12 at 16:23

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.


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

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.