2

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
  • What SRID are you using. I am guessing it is not 32648, which is possibly why you are not getting anything for rast_alignto. – John Powell Jan 25 '16 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? – Cushen Jan 26 '16 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. – tilt Jan 30 '16 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. – Cushen Feb 1 '16 at 3:51
0

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

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.