Added neighborhood hillshade on top of ST_MapAlgebraFctNgb. Part of #1318

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

diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 903e376787a2d9385ad3a135840f463f7e112be0..e1f7b2ffcdd2c994374d3b19edadf70c9f122b84 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -1984,6 +1984,10 @@ CREATE OR REPLACE FUNCTION _st_aspect4ma(matrix float[][], nodatamode text, vari
         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);
+        IF dz_dx = 0 AND dz_dy = 0 THEN
+            RETURN -1;
+        END IF;
+
         aspect := atan2(dz_dy, -dz_dx);
         IF aspect > (pi() / 2.0) THEN
             RETURN (5.0 * pi() / 2.0) - aspect;
@@ -2000,6 +2004,48 @@ CREATE OR REPLACE FUNCTION st_aspect(rast raster, band integer, pixeltype text)
     LANGUAGE 'SQL' STABLE;
 
 
+CREATE OR REPLACE FUNCTION _st_hillshade4ma(matrix float[][], nodatamode text, variadic args text[])
+    RETURNS float
+    AS
+    $$
+    DECLARE
+        pwidth float;
+        pheight float;
+        dz_dx float;
+        dz_dy float;
+        zenith float;
+        azimuth float;
+        slope float;
+        aspect float;
+        max_bright float;
+        elevation_scale float;
+    BEGIN
+        pwidth := args[1]::float;
+        pheight := args[2]::float;
+        azimuth := (5.0 * pi() / 2.0) - args[3]::float;
+        zenith := (pi() / 2.0) - args[4]::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);
+        elevation_scale := args[6]::float;
+        slope := atan(sqrt(elevation_scale * pow(dz_dx, 2.0) + pow(dz_dy, 2.0)));
+        aspect := atan2(dz_dy, -dz_dx);
+        max_bright := args[5]::float;
+
+        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 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;
+
+
 -----------------------------------------------------------------------
 -- Get information about the raster
 -----------------------------------------------------------------------

diff --git a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
index 788f0bbb4821fd166879c95d316a0e1c235feea2..abb0e1a63a900f697dea5a2f0af2454f0e4ee956 100644 (file)
--- a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
+++ b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql
@@ -123,10 +123,10 @@ SELECT
 -- test st_aspect, flat plane
 SELECT
   ST_Value(rast, 2, 2) = 1,
-  round(degrees(ST_Value(
+  ST_Value(
     ST_Aspect(rast, 1, NULL), 2, 2
-  ))*10000) = 2700000
-  FROM ST_TestRasterNgb(3, 3, 1) AS rast;
+  ) = -1
+  FROM ST_SetBandNoDataValue(ST_TestRasterNgb(3, 3, 1), NULL) AS rast;
 
 -- test st_aspect, North
 SELECT
@@ -159,3 +159,62 @@ SELECT
     ST_Aspect(rast, 1, NULL), 2, 2
   ))*10000) = 2700000
   FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 1, 2, 0) AS rast;
+
+-- test st_hillshade, flat plane
+SELECT
+  ST_Value(rast, 2, 2) = 1,
+  round(ST_Value(
+    ST_Hillshade(rast, 1, NULL, 0.0, pi()/4.0, 255), 2, 2
+  )*10000) = 1803122
+  FROM ST_TestRasterNgb(3, 3, 1) AS rast;
+
+-- test st_hillshade, known set from: http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=How%20Hillshade%20works
+SELECT
+  round(ST_Value(
+    ST_Hillshade(rast, 1, NULL, radians(315.0), pi()/4.0, 255, 1.0), 2, 2
+  )*10000) = 1540287
+  FROM ST_SetValue(
+    ST_SetValue(
+      ST_SetValue(
+        ST_SetValue(
+          ST_SetValue(
+            ST_SetValue(
+              ST_SetValue(
+                ST_SetValue(
+                  ST_SetValue(
+                    ST_SetScale(ST_TestRasterNgb(3, 3, 1), 5), 1, 1, 2450
+                  ), 2, 1, 2461
+                ), 3, 1, 2483
+              ), 1, 2, 2452
+            ), 2, 2, 2461
+          ), 3, 2, 2483
+        ), 1, 3, 2447
+      ), 2, 3, 2455
+    ), 3, 3, 2477
+  ) AS rast;
+
+-- test st_hillshade, defaults on known set from: http://webhelp.esri.com/arcgiSDEsktop/9.3/index.cfm?TopicName=How%20Hillshade%20works  
+SELECT
+  round(ST_Value(
+    ST_Hillshade(rast, 1, NULL, radians(315.0), pi()/4.0), 2, 2
+  )*10000) = 1540287
+  FROM ST_SetValue(
+    ST_SetValue(
+      ST_SetValue(
+        ST_SetValue(
+          ST_SetValue(
+            ST_SetValue(
+              ST_SetValue(
+                ST_SetValue(
+                  ST_SetValue(
+                    ST_SetScale(ST_TestRasterNgb(3, 3, 1), 5), 1, 1, 2450
+                  ), 2, 1, 2461
+                ), 3, 1, 2483
+              ), 1, 2, 2452
+            ), 2, 2, 2461
+          ), 3, 2, 2483
+        ), 1, 3, 2447
+      ), 2, 3, 2455
+    ), 3, 3, 2477
+  ) AS rast;
+

diff --git a/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected b/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
index c2c9bd8ab55f2d2ede4fe36ced78982be81e6f10..e73d999c6e784c124d1efbe7b4e767fe169d2f1a 100644 (file)
--- a/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
+++ b/raster/test/regress/rt_mapalgebrafctngb_userfunc_expected
@@ -26,3 +26,6 @@ t|t
 t|t
 t|t
 t|t
+t|t
+t
+t