Saturday, 14 November 2015

Zonal statistics of polygon overlays in PostGIS


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.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...