Added ST_Aspect map algebra neighborhood shortcut function. #1318

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

diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 68c344eb18559faea48f272548b9df59593363a9..903e376787a2d9385ad3a135840f463f7e112be0 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -1969,6 +1969,36 @@ CREATE OR REPLACE FUNCTION st_slope(rast raster, band integer, pixeltype text)
     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
+    AS
+    $$
+    DECLARE
+        pwidth float;
+        pheight float;
+        dz_dx float;
+        dz_dy float;
+        aspect 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);
+        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' 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;
+
 
 -----------------------------------------------------------------------
 -- Get information about the raster

diff --git a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
index b0bde6b6de121d0d06eb19afcba1599c6ca9a7e3..788f0bbb4821fd166879c95d316a0e1c235feea2 100644 (file)
--- a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
+++ b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
@@ -109,7 +109,7 @@ 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
+  ))*100000) = 5784902
   FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 1, 1, 10) AS rast;
 
 -- test st_slope, corner droop
@@ -117,5 +117,45 @@ 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
+  ))*100000) = 5784902
   FROM ST_SetValue(ST_TestRasterNgb(3, 3, 10), 1, 1, 1) AS rast;
+
+-- test st_aspect, flat plane
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(degrees(ST_Value(
+    ST_Aspect(rast, 1, NULL), 2, 2
+  ))*10000) = 2700000
+  FROM ST_TestRasterNgb(3, 3, 1) AS rast;
+
+-- test st_aspect, North
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(degrees(ST_Value(
+    ST_Aspect(rast, 1, NULL), 2, 2
+  ))*10000) = 0
+  FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 2, 1, 0) AS rast;
+
+-- test st_aspect, South
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(degrees(ST_Value(
+    ST_Aspect(rast, 1, NULL), 2, 2
+  ))*10000) = 1800000
+  FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 2, 3, 0) AS rast;
+
+-- test st_aspect, East
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(degrees(ST_Value(
+    ST_Aspect(rast, 1, NULL), 2, 2
+  ))*10000) = 900000
+  FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 3, 2, 0) AS rast;
+
+-- test st_aspect, West
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(degrees(ST_Value(
+    ST_Aspect(rast, 1, NULL), 2, 2
+  ))*10000) = 2700000
+  FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 1, 2, 0) AS rast;

diff --git a/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected b/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
index 997bc9fd68147e97b587155476c40fdcde5cebb8..c2c9bd8ab55f2d2ede4fe36ced78982be81e6f10 100644 (file)
--- a/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
+++ b/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
@@ -21,3 +21,8 @@ t|t
 t|t
 t|t
 t|t
+t|t
+t|t
+t|t
+t|t
+t|t