Function necessary for ST_Histogram with a geometry parameter

git-svn-id: http://svn.osgeo.org/postgis/trunk@7056 b70326c6-7e19-0410-871a-916f4a2858ee

diff --git a/raster/scripts/plpgsql/st_geomextent2rastercoord.sql b/raster/scripts/plpgsql/st_geomextent2rastercoord.sql
new file mode 100644 (file)
index 0000000..65f6ba1
--- /dev/null
+++ b/raster/scripts/plpgsql/st_geomextent2rastercoord.sql
@@ -0,0 +1,102 @@
+DROP FUNCTION ST_GeomExtent2RasterCoord(rast raster, geomin geometry);\r
+CREATE OR REPLACE FUNCTION ST_GeomExtent2RasterCoord(rast raster, \r
+                                                       geomin geometry)\r
+    RETURNS int[] AS\r
+    $$\r
+    DECLARE\r
+        geomintersect geometry;\r
+        x1w double precision := 0.0;\r
+        x2w double precision := 0.0;\r
+        y1w double precision := 0.0;\r
+        y2w double precision := 0.0;\r
+        x1 integer := 0;\r
+        x2 integer := 0;\r
+        x3 integer := 0;\r
+        x4 integer := 0;\r
+        y1 integer := 0;\r
+        y2 integer := 0;\r
+        y3 integer := 0;\r
+        y4 integer := 0;\r
+        psize float8;\r
+    BEGIN\r
+\r
+        -- Get the intersection between with the geometry.\r
+        -- We will search for withvalue pixel only in this area.\r
+        geomintersect := st_intersection(geomin, st_convexhull(rast));\r
+--RAISE NOTICE 'geomintersect1=%', astext(geomintersect);\r
+\r
+        -- If the intersection is empty, return false\r
+        IF st_isempty(geomintersect) THEN\r
+            RETURN ARRAY[NULL, NULL, NULL, NULL];\r
+        END IF;\r
+\r
+        -- We create a minimalistic buffer around the intersection in order to scan every pixels\r
+        -- that would touch the edge or intersect with the geometry\r
+        psize := st_scalex(rast) + st_skewy(rast);\r
+        geomintersect := st_buffer(geomintersect, psize / 1000000);\r
+\r
+--RAISE NOTICE 'geomintersect2=%', astext(geomintersect);\r
+\r
+        -- Find the world coordinates of the bounding box of the intersecting area\r
+        x1w := st_xmin(geomintersect);\r
+        y1w := st_ymin(geomintersect);\r
+        x2w := st_xmax(geomintersect);\r
+        y2w := st_ymax(geomintersect);\r
+\r
+--RAISE NOTICE 'x1w=%, y1w=%, x2w=%, y2w=%', x1w, y1w, x2w, y2w;\r
+\r
+        -- Convert world coordinates to raster coordinates\r
+        x1 := st_world2rastercoordx(rast, x1w, y1w);\r
+        y1 := st_world2rastercoordy(rast, x1w, y1w);\r
+        x2 := st_world2rastercoordx(rast, x2w, y1w);\r
+        y2 := st_world2rastercoordy(rast, x2w, y1w);\r
+        x3 := st_world2rastercoordx(rast, x1w, y2w);\r
+        y3 := st_world2rastercoordy(rast, x1w, y2w);\r
+        x4 := st_world2rastercoordx(rast, x2w, y2w);\r
+        y4 := st_world2rastercoordy(rast, x2w, y2w);\r
+\r
+--RAISE NOTICE 'x1=%, y1=%, x2=%, y2=%, x3=%, y3=%, x4=%, y4=%', x1, y1, x2, y2, x3, y3, x4, y4;\r
+\r
+        -- Order the raster coordinates for the upcoming FOR loop.\r
+        x1 := int4smaller(int4smaller(int4smaller(x1, x2), x3), x4);\r
+        y1 := int4smaller(int4smaller(int4smaller(y1, y2), y3), y4);\r
+        x2 := int4larger(int4larger(int4larger(x1, x2), x3), x4);\r
+        y2 := int4larger(int4larger(int4larger(y1, y2), y3), y4);\r
+\r
+        -- Make sure the range is not lower than 1.\r
+        -- This can happen when world coordinate are exactly on the left border\r
+        -- of the raster and that they do not span on more than one pixel.\r
+        x1 := int4smaller(int4larger(x1, 1), st_width(rast));\r
+        y1 := int4smaller(int4larger(y1, 1), st_height(rast));\r
+\r
+        -- Also make sure the range does not exceed the width and height of the raster.\r
+        -- This can happen when world coordinate are exactly on the lower right border\r
+        -- of the raster.\r
+        x2 := int4smaller(x2, st_width(rast));\r
+        y2 := int4smaller(y2, st_height(rast));\r
+        \r
+        RETURN ARRAY[x1, y1, x2, y2];\r
+\r
+    END;\r
+    $$\r
+    LANGUAGE 'plpgsql' IMMUTABLE STRICT;\r
+\r
+--Test\r
+SELECT (gvxy).geom, (gvxy).x::text || ',' || (gvxy).y::text\r
+FROM (SELECT DISTINCT rid, ST_PixelAsPolygons(rast) gvxy\r
+      FROM srtm_22_03_tiled, rect3\r
+      WHERE st_intersects(rast, geom)\r
+     ) foo\r
+\r
+SELECT DISTINCT rid\r
+FROM srtm_22_03_tiled, rect3\r
+WHERE st_intersects(rast, geom)\r
+\r
+\r
+SELECT rid, id, ST_GeomExtent2RasterCoord(rast, geom) \r
+FROM srtm_22_03_tiled, rect3\r
+WHERE st_intersects(rast, geom)\r
+\r
+\r
+SELECT (ext).x1, (ext).y1, (ext).x2, (ext).y2 FROM (SELECT ST_GeomExtent2RasterCoord(rast, geom) ext FROM srtm_22_03_tiled, rect3\r
+WHERE st_intersects(rast, geom)) foo
\ No newline at end of file