Friday, 7 April 2017

Raster statistics for polygon areas PostGIS


I'm trying to get zonal raster statistics for each polygon area. I have used the following method from this site:


SELECT bufid, (ST_SummaryStatsAgg(ST_Clip(rast, geom, true))).*
FROM rasttable, buftable
WHERE ST_Intersects(rast, geom)
GROUP BY bufid

This works great for a test area with 11 polygons. However, when I try to run it on the same raster with a larger table of polygons (count of 190,000) I keep getting the following error:




ERROR: rt_raster_from_two_rasters: The two rasters provided do not have the same alignment CONTEXT: PL/pgSQL function st_clip(raster,integer[],geometry,double precision[],boolean) line 8 at RETURN



How might I go about resolving this?


I'm using PostGIS 2.2.0 on PostgreSQL 9.4.4, Windows 7 64-bit.



Answer



I've come up with a solution which works:


  SELECT x.bufid, (ST_SummaryStatsAgg(x.intersectx,1,true)).* 
FROM
(SELECT bufid, ST_Intersection(rast,1,ST_AsRaster(geom, rast),1) as intersectx
FROM rasttable, buftable

WHERE ST_Intersects(geom, rast)) as x
WHERE x.intersectx IS NOT NULL
GROUP BY bufid

I think the issue was with the geometry of the vector data (buftable) whereby there was no intersection occurring between it and the raster data for some rows (due to invalid geometry?).


The following part takes care of that issue:



WHERE x.intersectx IS NOT NULL

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