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_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
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
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.