4

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.

Improve this question
4
  • What SRID are you using. I am guessing it is not 32648, which is possibly why you are not getting anything for rast_alignto. Commented Jan 25, 2016 at 21:13
  • Apologies, I am using 4326, I have changed the question to show that. I'm not sure that would affect the query either anyway? Commented Jan 26, 2016 at 21:50
  • I would be surprised if your raster is in 4326. If so, you might be better of transforming it to a cartesian projection first, along with your geometry. What is ST_Extent(rast) telling you? Your second query obviously doesn't work because the last FROM refers to ras69 instead of ref. Commented Jan 30, 2016 at 19:29
  • Thank you for your comments, I have discovered a solution and added this as an answer below. This didn't use the linked answer in the question so I have removed this part from the question for clarity sake. Commented Feb 1, 2016 at 3:51

1 Answer 1

Reset to default
1

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
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.