From 6d4a860c64054fd75e88d0eda322002162ebdb94 Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Fri, 12 Oct 2012 16:06:27 +0000 Subject: [PATCH] Duplicated and refactored the ST_XXX4ma() functions for ST_MapAlgebra usage. Exception for ST_InvDistWeight4ma() and ST_MinDist4ma(), both of which are new for 2.1. Added regression tests as well. git-svn-id: http://svn.osgeo.org/postgis/trunk@10419 b70326c6-7e19-0410-871a-916f4a2858ee --- doc/reference_raster.xml | 22 +- raster/rt_core/rt_api.c | 18 +- raster/rt_pg/rtpostgis.sql.in.c | 756 +++++++++++++++--- raster/test/regress/Makefile.in | 3 +- raster/test/regress/rt_4ma.sql | 82 ++ raster/test/regress/rt_4ma_expected | 17 + raster/test/regress/rt_invdistweight4ma.sql | 41 +- .../test/regress/rt_invdistweight4ma_expected | 31 +- 8 files changed, 820 insertions(+), 150 deletions(-) create mode 100644 raster/test/regress/rt_4ma.sql create mode 100644 raster/test/regress/rt_4ma_expected diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index 727d4939b..089aad6cf 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -3925,7 +3925,7 @@ FROM ( - The 2-D array output can be passed along to any of the raster processing builtin functions, e.g. ST_Min4ma, ST_Sum4ma, ST_Mean4ma. + The 2-D array output can be passed to any of the raster processing builtin functions, e.g. ST_Min4ma, ST_Sum4ma, ST_Mean4ma. Availability: 2.1.0 @@ -9566,9 +9566,9 @@ WHERE rid = 2; double precision ST_InvDistWeight4ma - double precision[][] matrix - text nodatamode - text[] VARIADIC args + double precision[][][] value + integer[][] pos + text[] VARIADIC userargs @@ -9579,7 +9579,7 @@ WHERE rid = 2; Calculate an interpolated value for a pixel using the Inverse Distance Weighted method. - There are two optional parameters that can be passed through args. The first parameter is the power factor (variable k in the equation below) between 0 and 1 used in the Inverse Distance Weighted equation. If not specified, default value is 1. The second parameter is the weight percentage applied only when the value of the pixel of interest is included with the interpolated value from the neighborhood. If not specified and the pixel of interest has a value, that value is returned. + There are two optional parameters that can be passed through userargs. The first parameter is the power factor (variable k in the equation below) between 0 and 1 used in the Inverse Distance Weighted equation. If not specified, default value is 1. The second parameter is the weight percentage applied only when the value of the pixel of interest is included with the interpolated value from the neighborhood. If not specified and the pixel of interest has a value, that value is returned. @@ -9600,7 +9600,7 @@ WHERE rid = 2; - This function is a specialized callback function for use as a callback parameter to . + This function is a specialized callback function for use as a callback parameter to . Availability: 2.1.0 @@ -9632,9 +9632,9 @@ WHERE rid = 2; double precision ST_MinDist4ma - double precision[][] matrix - text nodatamode - text[] VARIADIC args + double precision[][][] value + integer[][] pos + text[] VARIADIC userargs @@ -9651,7 +9651,7 @@ WHERE rid = 2; - This function is a specialized callback function for use as a callback parameter to . + This function is a specialized callback function for use as a callback parameter to . Availability: 2.1.0 @@ -9667,7 +9667,7 @@ WHERE rid = 2; See Also - , + , diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index b07f70382..35de2eb35 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -11496,6 +11496,9 @@ rt_raster_intersects( return 1; } } + else { + RASTER_DEBUGF(4, "GEOSIntersects() returned a 2!!!!"); + } } while (0); @@ -11638,8 +11641,10 @@ rt_raster_intersects( } } } + RASTER_DEBUG(4, "Smaller raster not in the other raster's pixel. Continuing"); } + RASTER_DEBUG(4, "Testing smaller raster vs larger raster"); *intersects = rt_raster_intersects_algorithm( rastS, rastL, bandS, bandL, @@ -11649,6 +11654,7 @@ rt_raster_intersects( if (*intersects) return 1; + RASTER_DEBUG(4, "Testing larger raster vs smaller raster"); *intersects = rt_raster_intersects_algorithm( rastL, rastS, bandL, bandS, @@ -13349,16 +13355,20 @@ rt_raster_iterator( } /* check that all rasters are aligned */ + RASTER_DEBUG(3, "checking alignment of all rasters"); rast = NULL; /* find raster to use as reference */ /* use custom if provided */ - if (extenttype == ET_CUSTOM) + if (extenttype == ET_CUSTOM) { + RASTER_DEBUG(4, "using custom extent as reference raster"); rast = customextent; + } /* use first valid one in _param->raster */ else { for (i = 0; i < itrcount; i++) { if (!_param->isempty[i]) { + RASTER_DEBUGF(4, "using raster at index %d as reference raster", i); rast = _param->raster[i]; break; } @@ -13387,6 +13397,7 @@ rt_raster_iterator( return NULL; } + RASTER_DEBUGF(5, "custom extent alignment: %d", aligned); if (!aligned) break; } @@ -13403,6 +13414,7 @@ rt_raster_iterator( return NULL; } + RASTER_DEBUGF(5, "raster at index %d alignment: %d", i, aligned); /* abort checking since a raster isn't aligned */ if (!aligned) @@ -13665,6 +13677,7 @@ rt_raster_iterator( /* input raster's X,Y */ x = _x - (int) _param->offset[i][0]; y = _y - (int) _param->offset[i][1]; + RASTER_DEBUGF(4, "source pixel (x, y) = (%d, %d)", x, y); _param->arg->src_pixel[i][0] = x; _param->arg->src_pixel[i][1] = y; @@ -13673,6 +13686,8 @@ rt_raster_iterator( npixels = NULL; status = 0; if (distancex > 0 && distancey > 0) { + RASTER_DEBUG(4, "getting neighborhood"); + status = rt_band_get_nearest_pixel( _param->band[i], x, y, @@ -13697,6 +13712,7 @@ rt_raster_iterator( (x >= 0 && x < rt_band_get_width(_param->band[i])) && (y >= 0 && y < rt_band_get_height(_param->band[i])) ) { + RASTER_DEBUG(4, "getting value of POI"); if (rt_band_get_pixel( _param->band[i], x, y, diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 9a7542b46..06775e0c9 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -2265,7 +2265,7 @@ CREATE OR REPLACE FUNCTION st_mapalgebrafctngb( LANGUAGE 'c' IMMUTABLE; ----------------------------------------------------------------------- --- Neighborhood MapAlgebra processing functions. +-- ST_MapAlgebraFctNgb() Neighborhood MapAlgebra processing functions. ----------------------------------------------------------------------- CREATE OR REPLACE FUNCTION st_max4ma(matrix float[][], nodatamode text, variadic args text[]) RETURNS float AS @@ -2433,10 +2433,12 @@ CREATE OR REPLACE FUNCTION _st_slope4ma(matrix float[][], nodatamode text, varia $$ LANGUAGE 'plpgsql' IMMUTABLE; +/* CREATE OR REPLACE FUNCTION st_slope(rast raster, band integer, pixeltype text) RETURNS RASTER AS $$ SELECT st_mapalgebrafctngb($1, $2, $3, 1, 1, '_st_slope4ma(float[][], text, text[])'::regprocedure, 'value', st_pixelwidth($1)::text, st_pixelheight($1)::text) $$ LANGUAGE 'sql' STABLE; +*/ CREATE OR REPLACE FUNCTION _st_aspect4ma(matrix float[][], nodatamode text, variadic args text[]) RETURNS float @@ -2467,10 +2469,12 @@ CREATE OR REPLACE FUNCTION _st_aspect4ma(matrix float[][], nodatamode text, vari $$ LANGUAGE 'plpgsql' IMMUTABLE; +/* CREATE OR REPLACE FUNCTION st_aspect(rast raster, band integer, pixeltype text) RETURNS RASTER AS $$ SELECT st_mapalgebrafctngb($1, $2, $3, 1, 1, '_st_aspect4ma(float[][], text, text[])'::regprocedure, 'value', st_pixelwidth($1)::text, st_pixelheight($1)::text) $$ LANGUAGE 'sql' STABLE; +*/ CREATE OR REPLACE FUNCTION _st_hillshade4ma(matrix float[][], nodatamode text, variadic args text[]) @@ -2515,10 +2519,12 @@ CREATE OR REPLACE FUNCTION _st_hillshade4ma(matrix float[][], nodatamode text, v $$ LANGUAGE 'plpgsql' IMMUTABLE; +/* CREATE OR REPLACE FUNCTION st_hillshade(rast raster, band integer, pixeltype text, azimuth float, altitude float, max_bright float DEFAULT 255.0, elevation_scale float DEFAULT 1.0) RETURNS RASTER AS $$ SELECT st_mapalgebrafctngb($1, $2, $3, 1, 1, '_st_hillshade4ma(float[][], text, text[])'::regprocedure, 'value', st_pixelwidth($1)::text, st_pixelheight($1)::text, $4::text, $5::text, $6::text, $7::text) $$ LANGUAGE 'sql' STABLE; +*/ CREATE OR REPLACE FUNCTION st_distinct4ma(matrix float[][], nodatamode TEXT, VARIADIC args TEXT[]) RETURNS float AS @@ -2530,15 +2536,410 @@ CREATE OR REPLACE FUNCTION st_stddev4ma(matrix float[][], nodatamode TEXT, VARIA $$ SELECT stddev(unnest) FROM unnest($1) $$ LANGUAGE 'sql' IMMUTABLE; +----------------------------------------------------------------------- +-- n-Raster ST_MapAlgebra +----------------------------------------------------------------------- + +CREATE OR REPLACE FUNCTION _st_mapalgebra( + rastbandargset rastbandarg[], + callbackfunc regprocedure, + pixeltype text DEFAULT NULL, + distancex integer DEFAULT 0, distancey integer DEFAULT 0, + extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, + VARIADIC userargs text[] DEFAULT NULL +) + RETURNS raster + AS 'MODULE_PATHNAME', 'RASTER_nMapAlgebra' + LANGUAGE 'c' STABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra( + rastbandargset rastbandarg[], + callbackfunc regprocedure, + pixeltype text DEFAULT NULL, + extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, + distancex integer DEFAULT 0, distancey integer DEFAULT 0, + VARIADIC userargs text[] DEFAULT NULL +) + RETURNS raster + AS $$ SELECT _ST_MapAlgebra($1, $2, $3, $6, $7, $4, $5, VARIADIC $8) $$ + LANGUAGE 'sql' STABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra( + rast raster, nband int[], + callbackfunc regprocedure, + pixeltype text DEFAULT NULL, + extenttype text DEFAULT 'FIRST', customextent raster DEFAULT NULL, + distancex integer DEFAULT 0, distancey integer DEFAULT 0, + VARIADIC userargs text[] DEFAULT NULL +) + RETURNS raster + AS $$ + DECLARE + x int; + argset rastbandarg[]; + BEGIN + IF $2 IS NULL OR array_ndims($2) < 1 OR array_length($2, 1) < 1 THEN + RAISE EXCEPTION 'Populated 1D array must be provided for nband'; + RETURN NULL; + END IF; + + FOR x IN array_lower($2, 1)..array_upper($2, 1) LOOP + IF $2[x] IS NULL THEN + CONTINUE; + END IF; + + argset := argset || ROW($1, $2[x])::rastbandarg; + END LOOP; + + IF array_length(argset, 1) < 1 THEN + RAISE EXCEPTION 'Populated 1D array must be provided for nband'; + RETURN NULL; + END IF; + + RETURN _ST_MapAlgebra(argset, $3, $4, $7, $8, $5, $6, VARIADIC $9); + END; + $$ LANGUAGE 'plpgsql' STABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra( + rast raster, nband int, + callbackfunc regprocedure, + pixeltype text DEFAULT NULL, + extenttype text DEFAULT 'FIRST', customextent raster DEFAULT NULL, + distancex integer DEFAULT 0, distancey integer DEFAULT 0, + VARIADIC userargs text[] DEFAULT NULL +) + RETURNS raster + AS $$ SELECT _ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], $3, $4, $7, $8, $5, $6, VARIADIC $9) $$ + LANGUAGE 'sql' STABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra( + rast1 raster, nband1 int, + rast2 raster, nband2 int, + callbackfunc regprocedure, + pixeltype text DEFAULT NULL, + extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, + distancex integer DEFAULT 0, distancey integer DEFAULT 0, + VARIADIC userargs text[] DEFAULT NULL +) + RETURNS raster + AS $$ SELECT _ST_MapAlgebra(ARRAY[ROW($1, $2), ROW($3, $4)]::rastbandarg[], $5, $6, $9, $10, $7, $8, VARIADIC $11) $$ + LANGUAGE 'sql' STABLE; + +----------------------------------------------------------------------- +-- ST_MapAlgebra callback functions +-- Should be called with values for distancex and distancey +-- These functions are meant for one raster +----------------------------------------------------------------------- + +-- helper function to convert 2D array to 3D array +CREATE OR REPLACE FUNCTION _st_convertarray4ma(value double precision[][]) + RETURNS double precision[][][] + AS $$ + DECLARE + _value double precision[][][]; + x int; + y int; + BEGIN + IF array_ndims(value) != 2 THEN + RAISE EXCEPTION 'Function parameter must be a 2-dimension array'; + END IF; + + _value := array_fill(NULL::double precision, ARRAY[1, array_length(value, 1), array_length(value, 2)]::int[], ARRAY[1, array_lower(value, 1), array_lower(value, 2)]::int[]); + + -- row + FOR y IN array_lower(value, 1)..array_upper(value, 1) LOOP + -- column + FOR x IN array_lower(value, 2)..array_upper(value, 2) LOOP + _value[1][y][x] = value[y][x]; + END LOOP; + END LOOP; + + RETURN _value; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT; + +CREATE OR REPLACE FUNCTION st_max4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + _value double precision[][][]; + max double precision; + x int; + y int; + z int; + ndims int; + BEGIN + max := '-Infinity'::double precision; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + -- raster + FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP + -- row + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + -- column + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP + IF _value[z][y][x] IS NULL THEN + IF array_length(userargs, 1) > 0 THEN + _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision; + ELSE + CONTINUE; + END IF; + END IF; + + IF _value[z][y][x] > max THEN + max := _value[z][y][x]; + END IF; + END LOOP; + END LOOP; + END LOOP; + + IF max = '-Infinity'::double precision THEN + RETURN NULL; + END IF; + + RETURN max; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_min4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + _value double precision[][][]; + min double precision; + x int; + y int; + z int; + ndims int; + BEGIN + min := 'Infinity'::double precision; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + -- raster + FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP + -- row + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + -- column + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP + IF _value[z][y][x] IS NULL THEN + IF array_length(userargs, 1) > 0 THEN + _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision; + ELSE + CONTINUE; + END IF; + END IF; + + IF _value[z][y][x] < min THEN + min := _value[z][y][x]; + END IF; + END LOOP; + END LOOP; + END LOOP; + + IF min = 'Infinity'::double precision THEN + RETURN NULL; + END IF; + + RETURN min; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_sum4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + _value double precision[][][]; + sum double precision; + x int; + y int; + z int; + ndims int; + BEGIN + sum := 0; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + -- raster + FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP + -- row + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + -- column + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP + IF _value[z][y][x] IS NULL THEN + IF array_length(userargs, 1) > 0 THEN + _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision; + ELSE + CONTINUE; + END IF; + END IF; + + sum := sum + _value[z][y][x]; + END LOOP; + END LOOP; + END LOOP; + + RETURN sum; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_mean4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + _value double precision[][][]; + sum double precision; + count int; + x int; + y int; + z int; + ndims int; + BEGIN + sum := 0; + count := 0; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + -- raster + FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP + -- row + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + -- column + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP + IF _value[z][y][x] IS NULL THEN + IF array_length(userargs, 1) > 0 THEN + _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision; + ELSE + CONTINUE; + END IF; + END IF; + + sum := sum + _value[z][y][x]; + count := count + 1; + END LOOP; + END LOOP; + END LOOP; + + IF count < 1 THEN + RETURN NULL; + END IF; + + RETURN sum / count::double precision; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_range4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + _value double precision[][][]; + min double precision; + max double precision; + x int; + y int; + z int; + ndims int; + BEGIN + min := 'Infinity'::double precision; + max := '-Infinity'::double precision; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + -- raster + FOR z IN array_lower(_value, 1)..array_upper(_value, 1) LOOP + -- row + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + -- column + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP + IF _value[z][y][x] IS NULL THEN + IF array_length(userargs, 1) > 0 THEN + _value[z][y][x] = userargs[array_lower(userargs, 1)]::double precision; + ELSE + CONTINUE; + END IF; + END IF; + + IF _value[z][y][x] < min THEN + min := _value[z][y][x]; + END IF; + IF _value[z][y][x] > max THEN + max := _value[z][y][x]; + END IF; + END LOOP; + END LOOP; + END LOOP; + + IF max = '-Infinity'::double precision OR min = 'Infinity'::double precision THEN + RETURN NULL; + END IF; + + RETURN max - min; + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_distinct4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ SELECT COUNT(DISTINCT unnest)::double precision FROM unnest($1) $$ + LANGUAGE 'sql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_stddev4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ SELECT stddev(unnest) FROM unnest($1) $$ + LANGUAGE 'sql' IMMUTABLE; + -- Inverse distance weight equation is based upon Equation 3.51 from the book -- Spatial Analysis A Guide for Ecologists -- by Marie-Josee Fortin and Mark Dale -- published by Cambridge University Press -- ISBN 0-521-00973-1 -CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodatamode text, VARIADIC args text[]) +CREATE OR REPLACE FUNCTION st_invdistweight4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) RETURNS double precision AS $$ DECLARE + _value double precision[][][]; + ndims int; + k double precision DEFAULT 1.; _k double precision DEFAULT 1.; z double precision[]; @@ -2546,6 +2947,7 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata _d double precision; z0 double precision; + _z integer; x integer; y integer; @@ -2559,17 +2961,28 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata max_dx double precision; max_dy double precision; BEGIN --- RAISE NOTICE 'matrix = %', matrix; --- RAISE NOTICE 'args = %', args; +-- RAISE NOTICE 'value = %', value; +-- RAISE NOTICE 'userargs = %', userargs; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; - -- make sure matrix is 2 dimensions - IF array_ndims(matrix) != 2 THEN - RAISE EXCEPTION 'Neighborhood array must be of two dimensions'; + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; END IF; + _z := array_lower(_value, 1); -- width and height (0-based) - w := array_upper(matrix, 2) - array_lower(matrix, 2); - h := array_upper(matrix, 1) - array_lower(matrix, 1); + h := array_upper(_value, 2) - array_lower(_value, 2); + w := array_upper(_value, 3) - array_lower(_value, 3); -- max distance from center pixel max_dx := w / 2; @@ -2590,14 +3003,14 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata END IF; -- center pixel's coordinates - cx := max_dx + array_lower(matrix, 2); - cy := max_dy + array_lower(matrix, 1); + cy := max_dy + array_lower(_value, 2); + cx := max_dx + array_lower(_value, 3); -- RAISE NOTICE 'cx, cy = %, %', cx, cy; - -- if args provided, only use the first two args - IF args IS NOT NULL AND array_ndims(args) = 1 THEN + -- if userargs provided, only use the first two args + IF userargs IS NOT NULL AND array_ndims(userargs) = 1 THEN -- first arg is power factor - k := args[array_lower(args, 1)]::double precision; + k := userargs[array_lower(userargs, 1)]::double precision; IF k IS NULL THEN k := _k; ELSEIF k < 0. THEN @@ -2610,8 +3023,8 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata -- second arg is what to do if center pixel has a value -- this will be a weight to apply for the center pixel - IF array_length(args, 1) > 1 THEN - cw := abs(args[array_lower(args, 1) + 1]::double precision); + IF array_length(userargs, 1) > 1 THEN + cw := abs(userargs[array_lower(userargs, 1) + 1]::double precision); IF cw IS NOT NULL THEN IF cw < 0. THEN RAISE NOTICE 'Weight (< 0) of center pixel value must be between 0 and 1. Defaulting to 0'; @@ -2627,24 +3040,24 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata k = abs(k) * -1; -- center pixel value - cv := matrix[cy][cx]; + cv := _value[_z][cy][cx]; -- check to see if center pixel has value -- RAISE NOTICE 'cw = %', cw; IF cw IS NULL AND cv IS NOT NULL THEN - RETURN matrix[cy][cx]; + RETURN cv; END IF; - FOR y IN array_lower(matrix, 1)..array_upper(matrix, 1) LOOP - FOR x IN array_lower(matrix, 2)..array_upper(matrix, 2) LOOP --- RAISE NOTICE 'matrix[%][%] = %', y, x, matrix[y][x]; + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP +-- RAISE NOTICE 'value[%][%][%] = %', _z, y, x, _value[_z][y][x]; -- skip NODATA values and center pixel - IF matrix[y][x] IS NULL OR (x = cx AND y = cy) THEN + IF _value[_z][y][x] IS NULL OR (x = cx AND y = cy) THEN CONTINUE; END IF; - z := z || matrix[y][x]; + z := z || _value[_z][y][x]; -- use pythagorean theorem _d := sqrt(power(cx - x, 2) + power(cy - y, 2)); @@ -2682,13 +3095,17 @@ CREATE OR REPLACE FUNCTION st_invdistweight4ma(matrix double precision[], nodata END; $$ LANGUAGE 'plpgsql' IMMUTABLE; -CREATE OR REPLACE FUNCTION st_mindist4ma(matrix double precision[], nodatamode text, VARIADIC args text[]) +CREATE OR REPLACE FUNCTION st_mindist4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) RETURNS double precision AS $$ DECLARE + _value double precision[][][]; + ndims int; + d double precision DEFAULT NULL; _d double precision; + z integer; x integer; y integer; @@ -2701,14 +3118,26 @@ CREATE OR REPLACE FUNCTION st_mindist4ma(matrix double precision[], nodatamode t max_dx double precision; max_dy double precision; BEGIN - -- make sure matrix is 2 dimensions - IF array_ndims(matrix) != 2 THEN - RAISE EXCEPTION 'Neighborhood array must be of two dimensions'; + + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; END IF; + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; + END IF; + z := array_lower(_value, 1); + -- width and height (0-based) - w := array_upper(matrix, 2) - array_lower(matrix, 2); - h := array_upper(matrix, 1) - array_lower(matrix, 1); + h := array_upper(_value, 2) - array_lower(_value, 2); + w := array_upper(_value, 3) - array_lower(_value, 3); -- max distance from center pixel max_dx := w / 2; @@ -2727,22 +3156,22 @@ CREATE OR REPLACE FUNCTION st_mindist4ma(matrix double precision[], nodatamode t END IF; -- center pixel's coordinates - cx := max_dx + array_lower(matrix, 2); - cy := max_dy + array_lower(matrix, 1); + cy := max_dy + array_lower(_value, 2); + cx := max_dx + array_lower(_value, 3); -- center pixel value - cv := matrix[cy][cx]; + cv := _value[z][cy][cx]; -- check to see if center pixel has value IF cv IS NOT NULL THEN RETURN 0.; END IF; - FOR y IN array_lower(matrix, 1)..array_upper(matrix, 1) LOOP - FOR x IN array_lower(matrix, 2)..array_upper(matrix, 2) LOOP + FOR y IN array_lower(_value, 2)..array_upper(_value, 2) LOOP + FOR x IN array_lower(_value, 3)..array_upper(_value, 3) LOOP -- skip NODATA values and center pixel - IF matrix[y][x] IS NULL OR (x = cx AND y = cy) THEN + IF _value[z][y][x] IS NULL OR (x = cx AND y = cy) THEN CONTINUE; END IF; @@ -2762,93 +3191,214 @@ CREATE OR REPLACE FUNCTION st_mindist4ma(matrix double precision[], nodatamode t $$ LANGUAGE 'plpgsql' IMMUTABLE; ----------------------------------------------------------------------- --- n-Raster ST_MapAlgebra +-- ST_Slope ----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION _st_slope4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + pwidth double precision; + pheight double precision; + dz_dx double precision; + dz_dy double precision; + + _value double precision[][][]; + ndims int; + z int; + BEGIN -CREATE OR REPLACE FUNCTION _st_mapalgebra( - rastbandargset rastbandarg[], - callbackfunc regprocedure, - pixeltype text DEFAULT NULL, - distancex integer DEFAULT 0, distancey integer DEFAULT 0, - extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, - VARIADIC userargs text[] DEFAULT NULL -) - RETURNS raster - AS 'MODULE_PATHNAME', 'RASTER_nMapAlgebra' - LANGUAGE 'c' STABLE; + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; -CREATE OR REPLACE FUNCTION st_mapalgebra( - rastbandargset rastbandarg[], - callbackfunc regprocedure, - pixeltype text DEFAULT NULL, - extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, - distancex integer DEFAULT 0, distancey integer DEFAULT 0, - VARIADIC userargs text[] DEFAULT NULL -) - RETURNS raster - AS $$ SELECT _ST_MapAlgebra($1, $2, $3, $6, $7, $4, $5, VARIADIC $8) $$ - LANGUAGE 'sql' STABLE; + IF ( + array_lower(_value, 2) != 1 OR array_upper(_value, 2) != 3 OR + array_lower(_value, 3) != 1 OR array_upper(_value, 3) != 3 + ) THEN + RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; + END IF; -CREATE OR REPLACE FUNCTION st_mapalgebra( - rast raster, nband int[], - callbackfunc regprocedure, - pixeltype text DEFAULT NULL, - extenttype text DEFAULT 'FIRST', customextent raster DEFAULT NULL, - distancex integer DEFAULT 0, distancey integer DEFAULT 0, - VARIADIC userargs text[] DEFAULT NULL -) + IF array_length(userargs, 1) < 2 THEN + RAISE EXCEPTION 'At least two elements must be provided for the third parameter'; + END IF; + + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; + END IF; + z := array_lower(_value, 1); + + pwidth := userargs[1]::double precision; + pheight := userargs[2]::double precision; + dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); + dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); + + RETURN atan(sqrt(pow(dz_dx, 2.0) + pow(dz_dy, 2.0))); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_slope(rast raster, band integer, pixeltype text) RETURNS raster + AS $$ SELECT ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_slope4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) $$ + LANGUAGE 'sql' IMMUTABLE; + +----------------------------------------------------------------------- +-- ST_Aspect +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION _st_aspect4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision AS $$ DECLARE - x int; - argset rastbandarg[]; + pwidth double precision; + pheight double precision; + dz_dx double precision; + dz_dy double precision; + aspect double precision; + + _value double precision[][][]; + ndims int; + z int; BEGIN - IF $2 IS NULL OR array_ndims($2) < 1 OR array_length($2, 1) < 1 THEN - RAISE EXCEPTION 'Populated 1D array must be provided for nband'; - RETURN NULL; + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; END IF; - FOR x IN array_lower($2, 1)..array_upper($2, 1) LOOP - IF $2[x] IS NULL THEN - CONTINUE; - END IF; + IF ( + array_lower(_value, 2) != 1 OR array_upper(_value, 2) != 3 OR + array_lower(_value, 3) != 1 OR array_upper(_value, 3) != 3 + ) THEN + RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; + END IF; - argset := argset || ROW($1, $2[x])::rastbandarg; - END LOOP; + IF array_length(userargs, 1) < 2 THEN + RAISE EXCEPTION 'At least two elements must be provided for the third parameter'; + END IF; - IF array_length(argset, 1) < 1 THEN - RAISE EXCEPTION 'Populated 1D array must be provided for nband'; - RETURN NULL; + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; END IF; + z := array_lower(_value, 1); - RETURN _ST_MapAlgebra(argset, $3, $4, $7, $8, $5, $6, VARIADIC $9); + pwidth := userargs[1]::double precision; + pheight := userargs[2]::double precision; + + dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); + dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); + IF abs(dz_dx) = 0::double precision AND abs(dz_dy) = 0::double precision THEN + RETURN -1; + END IF; + + aspect := atan2(dz_dy, -dz_dx); + IF aspect > (pi() / 2.0) THEN + RETURN (5.0 * pi() / 2.0) - aspect; + ELSE + RETURN (pi() / 2.0) - aspect; + END IF; END; - $$ LANGUAGE 'plpgsql' STABLE; + $$ LANGUAGE 'plpgsql' IMMUTABLE; -CREATE OR REPLACE FUNCTION st_mapalgebra( - rast raster, nband int, - callbackfunc regprocedure, - pixeltype text DEFAULT NULL, - extenttype text DEFAULT 'FIRST', customextent raster DEFAULT NULL, - distancex integer DEFAULT 0, distancey integer DEFAULT 0, - VARIADIC userargs text[] DEFAULT NULL -) +CREATE OR REPLACE FUNCTION st_aspect(rast raster, band integer, pixeltype text) RETURNS raster - AS $$ SELECT _ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], $3, $4, $7, $8, $5, $6, VARIADIC $9) $$ - LANGUAGE 'sql' STABLE; + AS $$ SELECT ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_aspect4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) $$ + LANGUAGE 'sql' IMMUTABLE; -CREATE OR REPLACE FUNCTION st_mapalgebra( - rast1 raster, nband1 int, - rast2 raster, nband2 int, - callbackfunc regprocedure, - pixeltype text DEFAULT NULL, - extenttype text DEFAULT 'INTERSECTION', customextent raster DEFAULT NULL, - distancex integer DEFAULT 0, distancey integer DEFAULT 0, - VARIADIC userargs text[] DEFAULT NULL +----------------------------------------------------------------------- +-- ST_HillShade +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION _st_hillshade4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) + RETURNS double precision + AS $$ + DECLARE + pwidth double precision; + pheight double precision; + dz_dx double precision; + dz_dy double precision; + zenith double precision; + azimuth double precision; + slope double precision; + aspect double precision; + max_bright double precision; + elevation_scale double precision; + + _value double precision[][][]; + ndims int; + z int; + BEGIN + ndims := array_ndims(value); + -- add a third dimension if 2-dimension + IF ndims = 2 THEN + _value := _st_convertarray4ma(value); + ELSEIF ndims != 3 THEN + RAISE EXCEPTION 'First parameter of function must be a 3-dimension array'; + ELSE + _value := value; + END IF; + + IF ( + array_lower(_value, 2) != 1 OR array_upper(_value, 2) != 3 OR + array_lower(_value, 3) != 1 OR array_upper(_value, 3) != 3 + ) THEN + RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; + END IF; + + IF array_length(userargs, 1) < 6 THEN + RAISE EXCEPTION 'At least six elements must be provided for the third parameter'; + END IF; + + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; + END IF; + z := array_lower(_value, 1); + + pwidth := userargs[1]::double precision; + pheight := userargs[2]::double precision; + + azimuth := (5.0 * pi() / 2.0) - userargs[3]::double precision; + zenith := (pi() / 2.0) - userargs[4]::double precision; + dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); + dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); + elevation_scale := userargs[6]::double precision; + slope := atan(sqrt(elevation_scale * pow(dz_dx, 2.0) + pow(dz_dy, 2.0))); + + -- handle special case of 0, 0 + IF abs(dz_dy) = 0::double precision AND abs(dz_dy) = 0::double precision THEN + -- set to pi as that is the expected PostgreSQL answer in Linux + aspect := pi(); + ELSE + aspect := atan2(dz_dy, -dz_dx); + END IF; + max_bright := userargs[5]::double precision; + + IF aspect < 0 THEN + aspect := aspect + (2.0 * pi()); + END IF; + + RETURN max_bright * ((cos(zenith) * cos(slope)) + (sin(zenith) * sin(slope) * cos(azimuth - aspect))); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_hillshade( + rast raster, band integer, + pixeltype text, + azimuth double precision, altitude double precision, max_bright double precision DEFAULT 255.0, elevation_scale double precision DEFAULT 1.0 ) - RETURNS raster - AS $$ SELECT _ST_MapAlgebra(ARRAY[ROW($1, $2), ROW($3, $4)]::rastbandarg[], $5, $6, $9, $10, $7, $8, VARIADIC $11) $$ - LANGUAGE 'sql' STABLE; + RETURNS RASTER + AS $$ SELECT ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_hillshade4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text, $4::text, $5::text, $6::text, $7::text) $$ + LANGUAGE 'sql' IMMUTABLE; ----------------------------------------------------------------------- -- Get information about the raster diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 74839735a..5c368470f 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -100,7 +100,8 @@ TEST_MAPALGEBRA = \ rt_clip \ rt_mapalgebra \ rt_union \ - rt_invdistweight4ma + rt_invdistweight4ma \ + rt_4ma TEST_GIST = \ rt_above \ diff --git a/raster/test/regress/rt_4ma.sql b/raster/test/regress/rt_4ma.sql new file mode 100644 index 000000000..08332791d --- /dev/null +++ b/raster/test/regress/rt_4ma.sql @@ -0,0 +1,82 @@ +DROP TABLE IF EXISTS raster_value_arrays; +CREATE TABLE raster_value_arrays ( + id integer, + val double precision[][] +); +CREATE OR REPLACE FUNCTION make_value_array( + rows integer DEFAULT 3, + columns integer DEFAULT 3, + start_val double precision DEFAULT 1, + step double precision DEFAULT 1, + skip_expr text DEFAULT NULL +) + RETURNS double precision[][][] + AS $$ + DECLARE + x int; + y int; + value double precision; + values double precision[][][]; + result boolean; + expr text; + BEGIN + value := start_val; + + values := array_fill(NULL::double precision, ARRAY[1, columns, rows]); + + FOR y IN 1..columns LOOP + FOR x IN 1..rows LOOP + IF skip_expr IS NULL OR length(skip_expr) < 1 THEN + result := TRUE; + ELSE + expr := replace(skip_expr, '[v]'::text, value::text); + EXECUTE 'SELECT (' || expr || ')::boolean' INTO result; + END IF; + + IF result IS TRUE THEN + values[1][y][x] := value; + END IF; + + value := value + step; + END LOOP; + END LOOP; + + RETURN values; + END; + $$ LANGUAGE 'plpgsql'; + +INSERT INTO raster_value_arrays VALUES + (1, make_value_array()), + (2, make_value_array(5, 5)), + (3, make_value_array(5, 5, 100)), + (4, make_value_array(3, 3, 15, -1)), + (5, make_value_array(5, 5, 15, -1)), + (6, make_value_array(3, 3, 1, 2)), + (7, make_value_array(5, 5, 1, 3)), + + (10, make_value_array(3, 3, 1, 1, '0')), + (11, make_value_array(5, 5, 1, 1, '0')), + (12, make_value_array(3, 3, 1, 1, '[v] % 2')), + (13, make_value_array(5, 5, 1, 1, '[v] % 2')), + (14, make_value_array(3, 3, 1, 1, '([v] % 2) = 0')), + (15, make_value_array(5, 5, 1, 1, '([v] % 2) = 0')), + (16, make_value_array(3, 3, 1, 2.1, '([v] NOT IN (7.3, 9.4, 15.7, 17.8))')), + (17, make_value_array(3, 3, 0, 3.14, '([v] IN (3.14, 12.56, 25.12))')), + (18, make_value_array(3, 3, 1, 1, '[v] > 8')) +; + +SELECT + id, + val, + round(st_distinct4ma(val, NULL)::numeric, 6) AS distinct4ma, + round(st_max4ma(val, NULL)::numeric, 6) AS max4ma, + round(st_mean4ma(val, NULL)::numeric, 6) AS mean4ma, + round(st_min4ma(val, NULL)::numeric, 6) AS min4ma, + round(st_range4ma(val, NULL)::numeric, 6) AS range4ma, + round(st_stddev4ma(val, NULL)::numeric, 6) AS stddev4ma, + round(st_sum4ma(val, NULL)::numeric, 6) AS sum4ma +FROM raster_value_arrays +ORDER BY id; + +DROP TABLE IF EXISTS raster_value_arrays; +DROP FUNCTION IF EXISTS make_value_array(integer, integer, double precision, double precision, text); diff --git a/raster/test/regress/rt_4ma_expected b/raster/test/regress/rt_4ma_expected new file mode 100644 index 000000000..d6ec1c2b2 --- /dev/null +++ b/raster/test/regress/rt_4ma_expected @@ -0,0 +1,17 @@ +NOTICE: table "raster_value_arrays" does not exist, skipping +1|{{{1,2,3},{4,5,6},{7,8,9}}}|9.000000|9.000000|5.000000|1.000000|8.000000|2.738613|45.000000 +2|{{{1,2,3,4,5},{6,7,8,9,10},{11,12,13,14,15},{16,17,18,19,20},{21,22,23,24,25}}}|25.000000|25.000000|13.000000|1.000000|24.000000|7.359801|325.000000 +3|{{{100,101,102,103,104},{105,106,107,108,109},{110,111,112,113,114},{115,116,117,118,119},{120,121,122,123,124}}}|25.000000|124.000000|112.000000|100.000000|24.000000|7.359801|2800.000000 +4|{{{15,14,13},{12,11,10},{9,8,7}}}|9.000000|15.000000|11.000000|7.000000|8.000000|2.738613|99.000000 +5|{{{15,14,13,12,11},{10,9,8,7,6},{5,4,3,2,1},{0,-1,-2,-3,-4},{-5,-6,-7,-8,-9}}}|25.000000|15.000000|3.000000|-9.000000|24.000000|7.359801|75.000000 +6|{{{1,3,5},{7,9,11},{13,15,17}}}|9.000000|17.000000|9.000000|1.000000|16.000000|5.477226|81.000000 +7|{{{1,4,7,10,13},{16,19,22,25,28},{31,34,37,40,43},{46,49,52,55,58},{61,64,67,70,73}}}|25.000000|73.000000|37.000000|1.000000|72.000000|22.079402|925.000000 +10|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}|0.000000||||||0.000000 +11|{{{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}}|0.000000||||||0.000000 +12|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|5.000000|9.000000|5.000000|1.000000|8.000000|3.162278|25.000000 +13|{{{1,NULL,3,NULL,5},{NULL,7,NULL,9,NULL},{11,NULL,13,NULL,15},{NULL,17,NULL,19,NULL},{21,NULL,23,NULL,25}}}|13.000000|25.000000|13.000000|1.000000|24.000000|7.788881|169.000000 +14|{{{NULL,2,NULL},{4,NULL,6},{NULL,8,NULL}}}|4.000000|8.000000|5.000000|2.000000|6.000000|2.581989|20.000000 +15|{{{NULL,2,NULL,4,NULL},{6,NULL,8,NULL,10},{NULL,12,NULL,14,NULL},{16,NULL,18,NULL,20},{NULL,22,NULL,24,NULL}}}|12.000000|24.000000|13.000000|2.000000|22.000000|7.211103|156.000000 +16|{{{1,3.1,5.2},{NULL,NULL,11.5},{13.6,NULL,NULL}}}|5.000000|13.600000|6.880000|1.000000|12.600000|5.435715|34.400000 +17|{{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}}|3.000000|25.120000|13.606667|3.140000|21.980000|11.027318|40.820000 +18|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,9}}}|1.000000|9.000000|9.000000|9.000000|0.000000||9.000000 diff --git a/raster/test/regress/rt_invdistweight4ma.sql b/raster/test/regress/rt_invdistweight4ma.sql index 926fe4138..8573d7031 100644 --- a/raster/test/regress/rt_invdistweight4ma.sql +++ b/raster/test/regress/rt_invdistweight4ma.sql @@ -10,19 +10,19 @@ CREATE OR REPLACE FUNCTION make_value_array( step double precision DEFAULT 1, skip_expr text DEFAULT NULL ) - RETURNS double precision[][] + RETURNS double precision[][][] AS $$ DECLARE x int; y int; value double precision; - values double precision[][]; + values double precision[][][]; result boolean; expr text; BEGIN value := start_val; - values := array_fill(NULL::double precision, ARRAY[columns, rows]); + values := array_fill(NULL::double precision, ARRAY[1, columns, rows]); FOR y IN 1..columns LOOP FOR x IN 1..rows LOOP @@ -34,7 +34,7 @@ CREATE OR REPLACE FUNCTION make_value_array( END IF; IF result IS TRUE THEN - values[y][x] := value; + values[1][y][x] := value; END IF; value := value + step; @@ -61,26 +61,29 @@ INSERT INTO raster_value_arrays VALUES (14, make_value_array(3, 3, 1, 1, '([v] % 2) = 0')), (15, make_value_array(5, 5, 1, 1, '([v] % 2) = 0')), (16, make_value_array(3, 3, 1, 2.1, '([v] NOT IN (7.3, 9.4, 15.7, 17.8))')), - (17, make_value_array(3, 3, 0, 3.14, '([v] IN (3.14, 12.56, 25.12))')) + (17, make_value_array(3, 3, 0, 3.14, '([v] IN (3.14, 12.56, 25.12))')), + (18, make_value_array(3, 3, 1, 1, '[v] > 8')) ; + SELECT id, val, - round(st_invdistweight4ma(val, NULL::text, NULL)::numeric, 6) AS test1, - round(st_invdistweight4ma(val, NULL::text, '0')::numeric, 6) AS test2, - round(st_invdistweight4ma(val, NULL::text, '0.5')::numeric, 6) AS test3, - round(st_invdistweight4ma(val, NULL::text, '0.9')::numeric, 6) AS test4, - round(st_invdistweight4ma(val, NULL::text, '1')::numeric, 6) AS test5, - round(st_invdistweight4ma(val, NULL::text, '0.9', '1')::numeric, 6) AS test6, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.9')::numeric, 6) AS test7, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.75')::numeric, 6) AS test8, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.5')::numeric, 6) AS test9, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.25')::numeric, 6) AS test10, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.1')::numeric, 6) AS test11, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.01')::numeric, 6) AS test12, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0.001')::numeric, 6) AS test13, - round(st_invdistweight4ma(val, NULL::text, '0.9', '0')::numeric, 6) AS test14 + round(st_invdistweight4ma(val, NULL, NULL)::numeric, 6) AS idw1, + round(st_invdistweight4ma(val, NULL, '0')::numeric, 6) AS idw2, + round(st_invdistweight4ma(val, NULL, '0.5')::numeric, 6) AS idw3, + round(st_invdistweight4ma(val, NULL, '0.9')::numeric, 6) AS idw4, + round(st_invdistweight4ma(val, NULL, '1')::numeric, 6) AS idw5, + round(st_invdistweight4ma(val, NULL, '0.9', '1')::numeric, 6) AS idw6, + round(st_invdistweight4ma(val, NULL, '0.9', '0.9')::numeric, 6) AS idw7, + round(st_invdistweight4ma(val, NULL, '0.9', '0.75')::numeric, 6) AS idw8, + round(st_invdistweight4ma(val, NULL, '0.9', '0.5')::numeric, 6) AS idw9, + round(st_invdistweight4ma(val, NULL, '0.9', '0.25')::numeric, 6) AS idw10, + round(st_invdistweight4ma(val, NULL, '0.9', '0.1')::numeric, 6) AS idw11, + round(st_invdistweight4ma(val, NULL, '0.9', '0.01')::numeric, 6) AS idw12, + round(st_invdistweight4ma(val, NULL, '0.9', '0.001')::numeric, 6) AS idw13, + round(st_invdistweight4ma(val, NULL, '0.9', '0')::numeric, 6) AS idw14, + round(st_mindist4ma(val, NULL)::numeric, 6) AS mindist4ma FROM raster_value_arrays ORDER BY id; diff --git a/raster/test/regress/rt_invdistweight4ma_expected b/raster/test/regress/rt_invdistweight4ma_expected index 48dffcebe..c21aebe79 100644 --- a/raster/test/regress/rt_invdistweight4ma_expected +++ b/raster/test/regress/rt_invdistweight4ma_expected @@ -1,16 +1,17 @@ NOTICE: table "raster_value_arrays" does not exist, skipping -1|{{1,2,3},{4,5,6},{7,8,9}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000 -2|{{1,2,3,4,5},{6,7,8,9,10},{11,12,13,14,15},{16,17,18,19,20},{21,22,23,24,25}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000 -3|{{100,101,102,103,104},{105,106,107,108,109},{110,111,112,113,114},{115,116,117,118,119},{120,121,122,123,124}}|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000 -4|{{15,14,13},{12,11,10},{9,8,7}}|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000 -5|{{15,14,13,12,11},{10,9,8,7,6},{5,4,3,2,1},{0,-1,-2,-3,-4},{-5,-6,-7,-8,-9}}|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000 -6|{{1,3,5},{7,9,11},{13,15,17}}|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000 -7|{{1,4,7,10,13},{16,19,22,25,28},{31,34,37,40,43},{46,49,52,55,58},{61,64,67,70,73}}|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000 -10|{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}|||||||||||||| -11|{{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}|||||||||||||| -12|{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000 -13|{{1,NULL,3,NULL,5},{NULL,7,NULL,9,NULL},{11,NULL,13,NULL,15},{NULL,17,NULL,19,NULL},{21,NULL,23,NULL,25}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000 -14|{{NULL,2,NULL},{4,NULL,6},{NULL,8,NULL}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000 -15|{{NULL,2,NULL,4,NULL},{6,NULL,8,NULL,10},{NULL,12,NULL,14,NULL},{16,NULL,18,NULL,20},{NULL,22,NULL,24,NULL}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000 -16|{{1,3.1,5.2},{NULL,NULL,11.5},{13.6,NULL,NULL}}|6.939697|6.880000|6.909550|6.933641|6.939697|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641 -17|{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}|12.560000|12.560000|12.560000|12.560000|12.560000|12.560000|12.546978|12.527446|12.494891|12.462337|12.442804|12.431085|12.429913|12.429783 +1|{{{1,2,3},{4,5,6},{7,8,9}}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|0.000000 +2|{{{1,2,3,4,5},{6,7,8,9,10},{11,12,13,14,15},{16,17,18,19,20},{21,22,23,24,25}}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|0.000000 +3|{{{100,101,102,103,104},{105,106,107,108,109},{110,111,112,113,114},{115,116,117,118,119},{120,121,122,123,124}}}|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|112.000000|0.000000 +4|{{{15,14,13},{12,11,10},{9,8,7}}}|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|11.000000|0.000000 +5|{{{15,14,13,12,11},{10,9,8,7,6},{5,4,3,2,1},{0,-1,-2,-3,-4},{-5,-6,-7,-8,-9}}}|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|3.000000|0.000000 +6|{{{1,3,5},{7,9,11},{13,15,17}}}|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|0.000000 +7|{{{1,4,7,10,13},{16,19,22,25,28},{31,34,37,40,43},{46,49,52,55,58},{61,64,67,70,73}}}|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|37.000000|0.000000 +10|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}||||||||||||||| +11|{{{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}}||||||||||||||| +12|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|0.000000 +13|{{{1,NULL,3,NULL,5},{NULL,7,NULL,9,NULL},{11,NULL,13,NULL,15},{NULL,17,NULL,19,NULL},{21,NULL,23,NULL,25}}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|0.000000 +14|{{{NULL,2,NULL},{4,NULL,6},{NULL,8,NULL}}}|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|5.000000|1.000000 +15|{{{NULL,2,NULL,4,NULL},{6,NULL,8,NULL,10},{NULL,12,NULL,14,NULL},{16,NULL,18,NULL,20},{NULL,22,NULL,24,NULL}}}|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|13.000000|1.000000 +16|{{{1,3.1,5.2},{NULL,NULL,11.5},{13.6,NULL,NULL}}}|6.939697|6.880000|6.909550|6.933641|6.939697|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|6.933641|1.000000 +17|{{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}}|12.560000|12.560000|12.560000|12.560000|12.560000|12.560000|12.546978|12.527446|12.494891|12.462337|12.442804|12.431085|12.429913|12.429783|0.000000 +18|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,9}}}|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|9.000000|1.414214 -- 2.40.0