From f17a5176a0a266e82d10138552dbc1a67167e54a Mon Sep 17 00:00:00 2001
From: Regina Obe <lr@pcorp.us>
Date: Fri, 16 Mar 2012 12:25:42 +0000
Subject: [PATCH] Provide summarizing summary stats example using intersection
 of geometry. Pretty fast.

git-svn-id: http://svn.osgeo.org/postgis/trunk@9509 b70326c6-7e19-0410-871a-916f4a2858ee
---
 doc/reference_raster.xml | 37 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 36 insertions(+), 1 deletion(-)

diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml
index b4cb2ef44..e2118eb90 100644
--- a/doc/reference_raster.xml
+++ b/doc/reference_raster.xml
@@ -4601,6 +4601,41 @@ FROM (SELECT rid, band, ST_SummaryStats(rast, band) As stats
 				</programlisting>	
 			</refsection>
 			
+			<refsection>
+				<title>Example: Summarize pixels that intersect buildings of interest</title>
+				<para>This example took 340ms on PostGIS windows 64-bit with all of Boston Buildings
+and aerial Tiles (tiles each 150x150 pixels ~ 134,000 tiles), ~102,000 building records</para>
+				<programlisting>WITH 
+-- our features of interest
+   feat AS (SELECT gid As building_id, geom_26986 As geom FROM buildings AS b 
+    WHERE gid IN(100, 103,150)
+   ),
+-- clip band 2 of raster tiles to boundaries of builds
+-- then get stats for these clipped regions
+   b_stats AS
+	(SELECT  building_id, (stats).*
+FROM (SELECT building_id, ST_SummaryStats(ST_Clip(rast,geom,2)) As stats
+    FROM aerials.boston
+		INNER JOIN feat
+	ON ST_Intersects(feat.geom,rast) 
+ ) As foo
+ )
+-- finally summarize stats
+SELECT building_id, SUM(count) As num_pixels
+  , MIN(min) As min_pval
+  ,  MAX(max) As max_pval
+  , SUM(mean*count)/SUM(count) As avg_pval
+	FROM b_stats
+ WHERE count > 0
+	GROUP BY building_id
+	ORDER BY building_id;
+ building_id | num_pixels | min_pval | max_pval |     avg_pval
+-------------+------------+----------+----------+------------------
+         100 |       1087 |        0 |      252 | 57.5501379944802
+         103 |        655 |        4 |      176 | 65.3145038167939
+         150 |        894 |        0 |      252 | 185.530201342282</programlisting>
+			</refsection>
+			
 			<refsection>
 				<title>Example: Raster coverage</title>
 				<programlisting>
@@ -4633,7 +4668,7 @@ FROM (SELECT band, ST_SummaryStats('o_4_boston','rast', band,true,0.25) As stats
 				<title>See Also</title>
 				<para>
 <xref linkend="RT_ST_Count" />,
-<xref linkend="summarystats" />
+<xref linkend="summarystats" />, <xref linkend="RT_ST_Clip" />
 				</para>
 			</refsection>
 		</refentry>
-- 
2.40.0