plpgsql implementation for st_clip.sql(raster, geom)

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

diff --git a/raster/scripts/plpgsql/st_clip.sql b/raster/scripts/plpgsql/st_clip.sql
index 53e6f7908c4da9662f32b640749c0edb54c6e259..19b6a2e43da473563201c1530656b3065ca5ef31 100644 (file)
--- a/raster/scripts/plpgsql/st_clip.sql
+++ b/raster/scripts/plpgsql/st_clip.sql
@@ -1,4 +1,4 @@
-----------------------------------------------------------------------
+----------------------------------------------------------------------
 --
 -- $Id$
 --
@@ -34,4 +34,87 @@ CREATE OR REPLACE FUNCTION ST_Clip(rast raster, x int, y int, width int, height
     $$
     LANGUAGE 'plpgsql';
 
+CREATE OR REPLACE FUNCTION ST_Clip(rast raster, band int, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false) 
+    RETURNS raster AS 
+    $$
+    DECLARE
+        sourceraster raster := rast;
+        newrast raster;
+        geomrast raster;
+        numband int := ST_Numbands(rast);
+        bandstart int;
+        bandend int;
+        newextent text;
+        newnodata float8;
+        bandi int;
+    BEGIN
+        IF rast IS NULL THEN
+            RETURN null;
+        END IF;
+        IF geom IS NULL THEN
+            RETURN rast;
+        END IF;
+        IF band IS NULL THEN
+            bandstart := 1;
+            bandend := numband;
+        ELSEIF ST_HasNoBand(rast, band) THEN
+            bandstart := numband;
+            bandend := numband;
+        ELSE
+            bandstart := band;
+            bandend := band;
+        END IF;
+        IF nodata IS NULL THEN
+            newnodata := coalesce(ST_BandNodataValue(rast, bandstart), ST_MinPossibleVal(ST_BandPixelType(rast, bandstart)));
+        ELSE
+            newnodata := nodata;
+        END IF;
+        newextent := CASE WHEN trimraster THEN 'INTERSECTION' ELSE 'FIRST' END;
+ RAISE NOTICE 'bandstart=%, bandend=%', bandstart, bandend;
+        -- Convert the geometry to a raster
+        geomrast := ST_AsRaster(geom, rast, ST_BandPixelType(rast, band), 1, newnodata);  
+
+        -- Compute the first raster band
+        sourceraster := ST_SetBandNodataValue(sourceraster, bandstart, newnodata);
+        newrast := ST_MapAlgebraExpr(sourceraster, bandstart, geomrast, 1, 'rast1');
+        
+        FOR bandi IN bandstart+1..bandend LOOP
+ RAISE NOTICE 'bandi=%', bandi;
+            -- for each band we must determine the nodata value
+            IF nodata IS NULL THEN
+                newnodata := coalesce(ST_BandNodataValue(sourceraster, bandi), ST_MinPossibleVal(ST_BandPixelType(rast, bandi)));
+            ELSE
+                newnodata := nodata;
+            END IF;
+            sourceraster := ST_SetBandNodataValue(sourceraster, bandi, newnodata);
+            newrast := ST_AddBand(newrast, ST_MapAlgebraExpr(sourceraster, bandi, geomrast, 1, 'rast1'));
+        END LOOP;
+        RETURN newrast;
+    END;
+    $$
+    LANGUAGE 'plpgsql';
+
+
+-- Test
+
+CREATE OR REPLACE FUNCTION ST_TestRaster(h integer, w integer, val float8) 
+    RETURNS raster AS 
+    $$
+    DECLARE
+    BEGIN
+        RETURN ST_AddBand(ST_MakeEmptyRaster(h, w, 0, 0, 1, 1, 0, 0, 0), '32BF', val, -1);
+    END;
+    $$
+    LANGUAGE 'plpgsql';
+
+SELECT ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4)) rast
+
+-- Test one band raster
+SELECT ST_AsBinary((gv).geom), (gv).val 
+FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4))) gv
+
+-- Test two bands raster
+SELECT ST_AsBinary((gv).geom), (gv).val 
+FROM ST_PixelAsPolygons(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), null, ST_Buffer(ST_MakePoint(8, 5), 4)), 2) gv
+