I am looking for a PostGIS equivalent to the ArcGIS Zonal Statistics by Table function.

I have two rasters. I would like to group the data in the second raster by the values in the first raster, and report the summary statistics (count, sum, etc.) as a table. There are questions for grouping raster data by polygons, but I don't see anything relevant to grouping by rasters.

The desired result is produced by this query, but I am working with a global 30-arc second resolution raster, and would like to avoid dumping all the values. Note that the two rasters are aligned and the rids match for intersecting tiles. (I assume that unnest preserves the order of arrays - is this correct?)

SELECT id, sum(val), count(val) from (
   SELECT unnest(ST_DumpValues(r1.rast,1)) id,
          unnest(ST_DumpValues(r2.rast,1)) val 
   FROM rast1 r1, rast2 r2 WHERE r1.rid=r2.rid) t1

I would dump the raster areas sharing a same value with ST_DumpAsPolygons() and perform an intersection between them and the raster coverage:

WITH areas AS (
  SELECT (gv).val id, (gv).geom geom
  FROM (SELECT ST_DumpAsPolygons(rast) gv
        FROM rast1) foo
) ss AS (
  SELECT id, ST_SummaryStats(ST_Clip(rast, geom)) ss
  FROM areas, rast2
  -- WHERE ST_Intersects(rast, geom)
  -- GROUP BY id
SELECT id, (ss).count cnt, (ss).sum sum;

That would be much faster if:

  • rast2 was tiled (20x20?) at load time

  • rast2 was indexed (and the WHERE clause uncommented)

  • you would use ST_SummaryStatsAgg() instead of ST_SummaryStats() (and uncomment the GROUP BY clause)

  • Thanks for the idea. I get several syntax errors with this query. For example, ST_Intersection(rast, geom) returns a setof geomval, which doesn't seem to be a valid argument type for ST_SummaryStats. But more generally, my query in the question works - I was just hoping for a solution that does not involve dumping out all the values. And I'm not clear why converting the raster to polygons would be any faster. – amball Aug 13 '17 at 19:03
  • Right I changed ST_Intersection() for ST_Clip(). I would give it a try to check for performance. – Pierre Racine Aug 22 '17 at 14:36

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.