7

I am trying to do a seemingly simple task -- extract mean pixel values from a raster based on a polygon overlay in PostGIS. I consulted several sources, including this blog, this blog, and this post). I loaded a Landsat tile into my database (in-db) and a polygon overlay consisting of three polygons (ID 1,2,3). The Landsat tile is a pansharpened image with a 15m spatial resolution and 4 bands.

Based on help files and the aforementioned posts/questions, I arrived at the following PostGIS query:

SELECT id, (SUM((ST_SummaryStats(a.rast, 1, true)).sum)/SUM((ST_SummaryStats(a.rast, 1, true)).count)) as mean FROM imagery.l8_2015_09_12 AS a, analysis_results.zonal_stats_test as b WHERE ST_Intersects(b.geom,a.rast) GROUP BY id;

The results are: 

1: 9902.49
2: 10079.68
3: 12355.90

However, when I test the results against ArcGIS's "zonal statistics to table" tool, I get the following results:

1: 8089.61
2: 7527.62
3: 12290.05

Close, but not identical. I am almost there but I was wondering if someone could help me figure out the source of the difference between both techniques. I suspect that there is something wrong with my ST_summaryStats method.

Improve this question
2
  • 1
    Does Zonal statistics of arcgis work with raster? If this answer is yes I belive that this is your answer. Sometimes, when you convert vector to raster you will lost informations of you data.
    – Diogo Caribé
    Dec 5, 2015 at 18:18
  • So perhaps a vector approach might be something to try. In that case, I would use a point overlay with st_value and then summarize within the boundaries of the polygon.
    – Cliff
    Dec 7, 2015 at 13:51

2 Answers 2

Reset to default
6

1) I'm surprised you do not "cut" the raster with ST_Clip(rast, geom) or ST_Intersection(rast, geom) before computing their stats. What you want are the stats on the pixels values intersecting the polygons no? That could be one reason why you end up with a different result but that would only be if only small portions of the rasters are cut, since your result is quite close to the ArcGIS one.

2) Even if you "cut" the raster properly before computing stats, the way PostGIS clip or intersect the raster might be different from ArcGIS: How intersecting pixels are chosen? By their centroid? By the total area intersecting? How nodata values are treated? How polygons are converted to raster before clipping/intersecting if it's the case? How rasters are converted to polygons if its the case? All this depend on the tools you select and sometimes the order in which you do the operations.

I suggest you work with a subset of your data and visualize/analyse the results step by step (in OpenJump for example) to understand the difference with ArcGIS.

See also this tutorial.

Improve this answer
5
  • Thanks for your comments. I used an intersect instead of a clip. In a previous version of this query I used a clip but it clipped to the 200x200 tile, not the pixel. So when I computed the mean I got the mean for every tile that intersected the geom. This query with a clip produces the exact same results:
    – Cliff
    Dec 7, 2015 at 15:06
  • --clip area of interest out of landsat tile WITH aoi AS ( SELECT a.*, st_clip(a.rast,b.geom) FROM imagery.l8_2015_09_12 AS a, analysis_results.zonal_stats_test as b WHERE ST_Intersects(b.geom,a.rast)) --sum and group statistics SELECT id, (SUM((ST_SummaryStats(a.rast, 1, true)).sum)/SUM((ST_SummaryStats(a.rast, 1, true)).count)) as mean --SELECT id, ST_SummaryStatsAgg(ST_Clip(a.rast, b.geom, true)) FROM aoi AS a, analysis_results.zonal_stats_test as b WHERE ST_Intersects(b.geom,a.rast) GROUP BY id;
    – Cliff
    Dec 7, 2015 at 15:06
  • 1
    You understand that using ST_Intersects(rast, geom) does not actually cut the raster right? For this you must use ST_Intersection() in the SELECT part...
    – Pierre Racine
    Dec 7, 2015 at 15:53
  • I have not tried st_intersection(). So I run st_intersection() first (instead of clip) to clip the raster, then run the original query on this clipped raster?
    – Cliff
    Dec 7, 2015 at 17:29
  • Yep! I would also investigate why ST_Clip() did not, apparently clip the raster...
    – Pierre Racine
    Dec 8, 2015 at 20:34
4

That could be the case of when you imported the Landsat raster to your database, you hopefully set a parameter -t 300x300, which splits the image in tiles of 300x300 pixels. So, imagine that you have your polygons in a shapefile that you imported to your database.

It is possible that you need more than one raster tile to cover a polygon feature, so, to prevent this issue you should add the St_union function to your query, just like this:

SELECT 
id,
(St_SummaryStats(St_Union(ST_Clip(rast,1,geom, true)))).*
FROM  imagery.l8_2015_09_12, analysis_results.zonal_stats_test 
WHERE st_intersects(rast,geom)
GROUP BY id;

Also, by using .* you get sum, count, max, min, stddev as bonus. If you want just the mean, change * for mean, that line would end up as this:

(St_SummaryStats(St_Union(ST_Clip(rast,1,geom, true)))).mean
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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