From: Bborie Park Date: Tue, 16 Oct 2012 19:55:26 +0000 (+0000) Subject: Added one-raster expression variant for ST_MapAlgebra and appropriate X-Git-Tag: 2.1.0beta2~530 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=2558772ba542c614d7d57dc5889fde3ae0681a71;p=postgis 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 06775e0c9..fd57d56e8 100644 --- 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 5c368470f..c41deb05a 100644 --- 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 index 000000000..41fa46190 --- /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 index 000000000..9eb5fd363 --- /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