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
$$
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
$$
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[])
$$
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
$$ 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[];
_d double precision;
z0 double precision;
+ _z integer;
x integer;
y integer;
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;
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
-- 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';
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));
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;
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;
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;
$$ 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
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
projects / postgis / commitdiff