$$
LANGUAGE 'plpgsql' IMMUTABLE;
+CREATE OR REPLACE FUNCTION _st_slope4ma(matrix float[][], nodatamode text, variadic args text[])
+ RETURNS float
+ AS
+ $$
+ DECLARE
+ pwidth float;
+ pheight float;
+ dz_dx float;
+ dz_dy float;
+ BEGIN
+ pwidth := args[1]::float;
+ pheight := args[2]::float;
+ dz_dx := ((matrix[3][1] + 2.0 * matrix[3][2] + matrix[3][3]) - (matrix[1][1] + 2.0 * matrix[1][2] + matrix[1][3])) / (8.0 * pwidth);
+ dz_dy := ((matrix[1][3] + 2.0 * matrix[2][3] + matrix[3][3]) - (matrix[1][1] + 2.0 * matrix[2][1] + matrix[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_mapalgebrafctngb($1, $2, $3, 1, 1, '_st_slope4ma(float[][], text, text[])'::regprocedure, 'value', st_pixelwidth($1)::text, st_pixelheight($1)::text) $$
+ LANGUAGE 'SQL' STABLE;
+
-----------------------------------------------------------------------
-- Get information about the raster
) = 9
FROM ST_SetValue(ST_TestRasterNgb(3, 3, 10), 3, 3, NULL) AS rast;
+-- test st_slope, flat plane
+SELECT
+ ST_Value(rast, 2, 2) = 1,
+ ST_Value(
+ ST_Slope(rast, 1, NULL), 2, 2
+ ) = 0
+ FROM ST_TestRasterNgb(3, 3, 1) AS rast;
+
+-- test st_slope, corner spike
+SELECT
+ ST_Value(rast, 2, 2) = 1,
+ round(degrees(ST_Value(
+ ST_Slope(rast, 1, NULL), 2, 2
+ ))*100000) = 5784902 -- this measurement is in radians
+ FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 1, 1, 10) AS rast;
+
+-- test st_slope, corner droop
+SELECT
+ ST_Value(rast, 2, 2) = 10,
+ round(degrees(ST_Value(
+ ST_Slope(rast, 1, NULL), 2, 2
+ ))*100000) = 5784902 -- this measurement is in radians
+ FROM ST_SetValue(ST_TestRasterNgb(3, 3, 10), 1, 1, 1) AS rast;