Added ST_Slope map algebra neighborhood shortcut function.

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

raster/rt_pg/rtpostgis.sql.in.c

@@ -1945,6 +1945,30 @@ CREATE OR REPLACE FUNCTION st_range4ma(matrix float[][], nodatamode text, variad
     $$
     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

raster/test/regress/rt_mapalgebrafctngb_userfunc.sql

@@ -96,3 +96,26 @@ SELECT
   ) = 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;

raster/test/regress/rt_mapalgebrafctngb_userfunc_expected

@@ -18,3 +18,6 @@ t|t
 NOTICE:  Neighborhood NODATA behavior defaulting to 'ignore'
 t|t
 t|t
+t|t
+t|t
+t|t