From: David Zwarg Date: Thu, 1 Dec 2011 04:29:24 +0000 (+0000) Subject: Added neighborhood hillshade on top of ST_MapAlgebraFctNgb. Part of #1318 X-Git-Tag: 2.0.0alpha1~591 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=bc11c03245a420692a5a7328e2741527013ea98b;p=postgis 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 903e37678..e1f7b2ffc 100644 --- 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 788f0bbb4..abb0e1a63 100644 --- 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 c2c9bd8ab..e73d999c6 100644 --- 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