Added one-raster expression variant for ST_MapAlgebra and appropriate

regression tests.

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

diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 06775e0c90509177aac74257ee13ce5226b42cef..fd57d56e8f973b84c246c3082df0fa690559825d 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -2625,6 +2625,85 @@ CREATE OR REPLACE FUNCTION st_mapalgebra(
        AS $$ SELECT _ST_MapAlgebra(ARRAY[ROW($1, $2), ROW($3, $4)]::rastbandarg[], $5, $6, $9, $10, $7, $8, VARIADIC $11) $$
        LANGUAGE 'sql' STABLE;
 
+CREATE OR REPLACE FUNCTION st_mapalgebra(
+       rast raster, nband integer,
+       pixeltype text,
+       expression text, nodataval double precision DEFAULT NULL
+)
+       RETURNS raster
+       AS $$
+       DECLARE
+               rtn raster;
+               funcname text;
+               funcbody text;
+               funccmd text;
+               expr text;
+       BEGIN
+               -- set name of callback function
+               funcname := 'pg_temp.callbackfunc' || (round(random() * 1000000))::text;
+
+               -- substitute keywords
+               expr := replace(expression, '[rast]', ''' || value[1][1][1] || ''');
+               expr := replace(expr, '[rast.val]', ''' || value[1][1][1] || ''');
+               expr := replace(expr, '[rast.x]', ''' || pos[1][1] || ''');
+               expr := replace(expr, '[rast.y]', ''' || pos[1][2] || ''');
+
+               -- build callback function body
+               funcbody := ' DECLARE val double precision; BEGIN';
+               --funcbody := funcbody || ' RAISE NOTICE ''value = %'', value;';
+               --funcbody := funcbody || ' RAISE NOTICE ''pos = %'', pos;';
+               funcbody := funcbody || ' IF value[1][1][1] IS NULL THEN';
+
+               -- handle NODATA in callback function
+               IF nodataval IS NOT NULL THEN
+                       funcbody := funcbody || ' RETURN ' || nodataval::text || ';';
+               ELSE
+                       funcbody := funcbody || ' RETURN NULL;';
+               END IF;
+               funcbody := funcbody || ' END IF;';
+
+               -- handle value processing in callback function
+               funcbody := funcbody || ' EXECUTE ''SELECT (' || expr || ')::double precision'' INTO val;';
+
+               -- handle val is NULL
+               IF nodataval IS NOT NULL THEN
+                       funcbody := funcbody || ' IF val IS NULL THEN RETURN ' || nodataval::text || '; END IF;';
+               END IF;
+
+               -- finish callback function
+               funcbody := funcbody || ' RETURN val; END;';
+
+               funccmd := 'CREATE OR REPLACE FUNCTION ' || funcname
+                       || '(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL)'
+                       || ' RETURNS double precision AS '
+                       || quote_literal(funcbody)
+                       || ' LANGUAGE ''plpgsql'' STABLE;';
+
+               --RAISE NOTICE 'funccmd = %', funccmd;
+
+               -- create callback function
+               EXECUTE funccmd;
+
+               -- call ST_MapAlgebra
+               rtn := ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], (funcname || '(double precision[][][], integer[][], text[])')::regprocedure, $3, 'FIRST');
+
+               -- drop callback function
+               EXECUTE 'DROP FUNCTION IF EXISTS ' || funcname
+                       || '(double precision[][][], integer[][], text[]);';
+
+               RETURN rtn;
+       END;
+       $$ LANGUAGE 'plpgsql' VOLATILE;
+
+CREATE OR REPLACE FUNCTION st_mapalgebra(
+       rast raster,
+       pixeltype text,
+       expression text, nodataval double precision DEFAULT NULL
+)
+       RETURNS raster
+       AS $$ SELECT st_mapalgebra($1, 1, $2, $3, $4) $$
+       LANGUAGE 'sql' VOLATILE;
+
 -----------------------------------------------------------------------
 -- ST_MapAlgebra callback functions
 -- Should be called with values for distancex and distancey

diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 5c368470f0712ab0f711cba8a9f04b999b0a3f61..c41deb05a256e7b0104e97842e60876c281f0282 100644 (file)
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -99,6 +99,7 @@ TEST_MAPALGEBRA = \
        rt_intersection \
        rt_clip \
        rt_mapalgebra \
+       rt_mapalgebra_expr \
        rt_union \
        rt_invdistweight4ma \
        rt_4ma

diff --git a/raster/test/regress/rt_mapalgebra_expr.sql b/raster/test/regress/rt_mapalgebra_expr.sql
new file mode 100644 (file)
index 0000000..41fa461
--- /dev/null
+++ b/raster/test/regress/rt_mapalgebra_expr.sql
@@ -0,0 +1,164 @@
+SET client_min_messages TO warning;
+
+CREATE OR REPLACE FUNCTION ST_TestRaster(ulx float8, uly float8, val float8) 
+    RETURNS raster AS 
+    $$
+    DECLARE
+    BEGIN
+        RETURN ST_AddBand(ST_MakeEmptyRaster(10, 10, ulx, uly, 1, 1, 0, 0, 0), '32BF', val, -1);
+    END;
+    $$
+    LANGUAGE 'plpgsql';
+
+CREATE OR REPLACE FUNCTION raster_plus_twenty(pixel FLOAT, VARIADIC args TEXT[])
+    RETURNS FLOAT AS 
+    $$
+    BEGIN
+        IF pixel IS NULL THEN
+            RAISE NOTICE 'Pixel value is null.';
+        END IF;
+        RETURN pixel + 20;
+    END;
+    $$ 
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION raster_plus_arg1(pixel FLOAT, VARIADIC args TEXT[])
+    RETURNS FLOAT AS 
+    $$
+    DECLARE
+        x float := 0;
+    BEGIN
+        IF NOT args[1] IS NULL THEN
+            x := args[1]::float;
+        END IF;
+        RETURN pixel + x;
+    END;
+    $$ 
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION raster_polynomial(pixel FLOAT, VARIADIC args TEXT[])
+    RETURNS FLOAT AS
+    $$
+    DECLARE
+        m float := 1;
+        b float := 0;
+    BEGIN
+        IF NOT args[1] is NULL THEN
+            m := args[1]::float;
+        END IF;
+        IF NOT args[2] is NULL THEN
+            b := args[2]::float;
+        END IF;
+        RETURN m * pixel + b;
+    END;
+    $$
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+ CREATE OR REPLACE FUNCTION raster_nullage(pixel FLOAT, VARIADIC args TEXT[])
+    RETURNS FLOAT AS
+    $$
+    BEGIN
+        RETURN NULL;
+    END;
+    $$
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+ CREATE OR REPLACE FUNCTION raster_x_plus_arg(pixel FLOAT, pos INT[], VARIADIC args TEXT[])
+    RETURNS FLOAT AS
+    $$
+    DECLARE
+        x float := 0;
+    BEGIN
+        IF NOT args[1] IS NULL THEN
+            x := args[1]::float;
+        END IF;
+        RETURN pixel + pos[1] + x;
+    END;
+    $$
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+ CREATE OR REPLACE FUNCTION raster_y_plus_arg(pixel FLOAT, pos INT[], VARIADIC args TEXT[])
+    RETURNS FLOAT AS
+    $$
+    DECLARE
+        x float := 0;
+    BEGIN
+        IF NOT args[1] IS NULL THEN
+            x := args[1]::float;
+        END IF;
+        RETURN pixel + pos[2] + x;
+    END;
+    $$
+    LANGUAGE 'plpgsql' IMMUTABLE;
+
+-- Test NULL raster
+SELECT ST_MapAlgebra(NULL, 1, NULL, '[rast] + 20', 2) IS NULL FROM ST_TestRaster(0, 0, -1) rast;
+
+-- Test empty raster
+SELECT 'T1', ST_IsEmpty(ST_MapAlgebra(ST_MakeEmptyRaster(0, 10, 0, 0, 1, 1, 1, 1, 0), 1, NULL, '[rast] + 20', 2));
+
+-- Test hasnoband raster
+SELECT 'T2', ST_HasNoBand(ST_MapAlgebra(ST_MakeEmptyRaster(10, 10, 0, 0, 1, 1, 1, 1, 0), 1, NULL, '[rast] + 20', 2));
+
+-- Test hasnodata value
+SELECT 'T3', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(ST_SetBandNoDataValue(rast, NULL), 1, NULL, '[rast] + 20', 2), 1, 1) FROM ST_TestRaster(0, 0, -1) rast;
+
+-- Test nodata value
+SELECT 'T4', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, NULL, '[rast] + 20', 2), 1, 1) FROM ST_TestRaster(0, 0, -1) rast;
+
+-- Test 'rast' expression
+SELECT 'T5', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, NULL, '[rast]', 2), 1, 1) FROM ST_TestRaster(0, 0, -1) rast;
+SELECT 'T6', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(ST_SetBandNoDataValue(rast, NULL), 1, NULL, '[rast]'), 1, 1) FROM ST_TestRaster(0, 0, -1) rast;
+
+-- Test pixeltype
+SELECT 'T7', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, '4BUI', '[rast] + 20', 2), 1, 1) FROM ST_TestRaster(0, 0, 100) rast;
+SELECT 'T8', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, '4BUId', '[rast] + 20', 2), 1, 1) FROM ST_TestRaster(0, 0, 100) rast;
+SELECT 'T9', ST_Value(rast, 1, 1), ST_Value(ST_MapAlgebra(rast, 1, '2BUI', '[rast] + 20', 2), 1, 1) FROM ST_TestRaster(0, 0, 101) rast;
+
+-- Test '[rast.x]', '[rast.y]' and '[rast.val]' substitutions expression
+SELECT 'T10.'||x||'.'||y, ST_Value(rast, x, y),
+  ST_Value(ST_MapAlgebra(rast, 1, NULL, 'ceil([rast]*[rast.x]/[rast.y]+[rast.val])'), x, y)
+FROM ST_TestRaster(0, 0, 10) rast,
+     generate_series(6, 8) as x,
+     generate_series(2, 4) as y
+ORDER BY x, y;
+
+-- Test divide by zero; blows up the whole query, not just the SPI
+SELECT ST_Value(rast, 1, 1),
+    ST_Value(ST_MapAlgebra(rast, 1, NULL, '[rast]/0'), 1, 1)
+FROM ST_TestRaster(0, 0, 10) rast;
+
+-- Test evaluations to null (see #1523)
+WITH op AS ( select rast AS rin,
+  ST_MapAlgebra(rast, 1, NULL, 'SELECT g from (SELECT NULL::double precision as g) as foo', 2)
+  AS rout FROM ST_TestRaster(0, 0, 10) rast )
+SELECT 'T11.1', ST_Value(rin, 1, 1), ST_Value(rout, 1, 1) FROM op;
+WITH op AS ( select rast AS rin,
+  ST_MapAlgebra(rast, 1, NULL, 'SELECT g from (select [rast],NULL::double precision as g) as foo', 2)
+  AS rout FROM ST_TestRaster(0, 0, 10) rast )
+SELECT 'T11.2', ST_Value(rin, 1, 1), ST_Value(rout, 1, 1) FROM op;
+
+-- Test pracine's new bug #1638
+SELECT 'T12',
+  ST_Value(rast, 1, 2) = 1,
+  ST_Value(rast, 2, 1) = 2,
+  ST_Value(rast, 4, 3) = 4,
+  ST_Value(rast, 3, 4) = 3
+  FROM ST_MapAlgebra(
+    ST_AddBand(
+      ST_MakeEmptyRaster(10, 10, 0, 0, 0.001, 0.001, 0, 0, 4269), 
+      '8BUI'::text, 
+      1, 
+      0
+    ), 
+    '32BUI', 
+    '[rast.x]'
+  ) AS rast; 
+
+DROP FUNCTION ST_TestRaster(ulx float8, uly float8, val float8);
+DROP FUNCTION raster_plus_twenty(pixel FLOAT, VARIADIC args TEXT[]);
+DROP FUNCTION raster_plus_arg1(pixel FLOAT, VARIADIC args TEXT[]);
+DROP FUNCTION raster_polynomial(pixel FLOAT, VARIADIC args TEXT[]);
+DROP FUNCTION raster_nullage(pixel FLOAT, VARIADIC args TEXT[]);
+DROP FUNCTION raster_x_plus_arg(pixel FLOAT, pos INT[], VARIADIC args TEXT[]);
+DROP FUNCTION raster_y_plus_arg(pixel FLOAT, pos INT[], VARIADIC args TEXT[]);

diff --git a/raster/test/regress/rt_mapalgebra_expr_expected b/raster/test/regress/rt_mapalgebra_expr_expected
new file mode 100644 (file)
index 0000000..9eb5fd3
--- /dev/null
+++ b/raster/test/regress/rt_mapalgebra_expr_expected
@@ -0,0 +1,23 @@
+t
+T1|t
+T2|t
+T3||19
+T4||2
+T5||2
+T6||-1
+T7|100|15
+ERROR:  RASTER_nMapAlgebra: Invalid pixel type: 4BUId
+T9|101|3
+T10.6.2|10|40
+T10.6.3|10|30
+T10.6.4|10|25
+T10.7.2|10|45
+T10.7.3|10|33
+T10.7.4|10|27
+T10.8.2|10|50
+T10.8.3|10|36
+T10.8.4|10|30
+ERROR:  division by zero
+T11.1|10|2
+T11.2|10|2
+T12|t|t|t|t