Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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:

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.

share|improve this question
    
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
1  
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
2  
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

2 Answers 2

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

share|improve this answer
    
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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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