From 596a985b12e024c7e84723e86957b1323361040f Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Fri, 2 Nov 2012 17:22:29 +0000 Subject: [PATCH] Refactored ST_Slope, ST_Aspect and ST_Hillshade() (also their _st_XXX4ma() functions). Detailed docs for all three functions. Outputs now in sync with return from ArcGIS. Ticket is #2077 git-svn-id: http://svn.osgeo.org/postgis/trunk@10619 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 7 +- doc/reference_raster.xml | 225 ++++++---- raster/rt_pg/rtpostgis.sql.in.c | 412 ++++++++++++++---- raster/rt_pg/rtpostgis_drop.sql.in.c | 4 + .../regress/rt_mapalgebrafctngb_userfunc.sql | 34 +- 5 files changed, 510 insertions(+), 172 deletions(-) diff --git a/NEWS b/NEWS index b9783bec6..5d50104d8 100644 --- a/NEWS +++ b/NEWS @@ -8,6 +8,9 @@ PostGIS 2.1.0 - #2026, ST_Union(raster) now unions all bands of all rasters - point variant of ST_SetValue(raster) previously did not check SRID of input geometry and raster. + - ST_Hillshade parameters azimuth and altitude are now in degrees + instead of radians. + - ST_Slope and ST_Aspect return pixel values in degrees instead of radians * Deprecated signatures @@ -44,7 +47,6 @@ PostGIS 2.1.0 (Bborie Park / UC Davis) - ST_Tile(raster) to break up a raster into tiles (Bborie Park / UC Davis) - #739, UpdateRasterSRID() - - #1655, Addition default values for parameters of ST_Slope * Enhancements * - #823, tiger geocoder: Make loader_generate_script download portion @@ -73,6 +75,8 @@ PostGIS 2.1.0 - point variant of ST_SetValue(raster) is now a wrapper around geomval variant of ST_SetValues(rast). - proper support for raster band's isnodata flag in core API and loader. + - #1655, Additional default values for parameters of ST_Slope + - Additional default values for parameters of ST_Aspect and ST_HillShade * Fixes * @@ -92,6 +96,7 @@ PostGIS 2.1.0 - #1981, cleanup of unused variables causing warnings with gcc 4.6+ - #2026, improve performance of distance calculations - #2057, Fixed linking issue for raster2psql to libpq + - #2077, Fixed incorrect values returning from ST_Hillshade() PostGIS 2.0.1 2012/06/22 diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index ee3cb0ef0..e89f87bae 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -7691,7 +7691,7 @@ FROM dummy_rast; ST_HillShade - Returns the hypothetical illumination of an elevation raster band using provided azimuth, altitude, brightness and elevation scale inputs. Useful for visualizing terrain. + Returns the hypothetical illumination of an elevation raster band using provided azimuth, altitude, brightness and scale inputs. @@ -7700,10 +7700,10 @@ FROM dummy_rast; raster rast integer band text pixeltype - double precision azimuth - double precision altitude + double precision azimuth=315 + double precision altitude=45 double precision max_bright=255 - double precision elevation_scale=1 + double precision scale=1.0 boolean interpolate_nodata=FALSE @@ -7712,22 +7712,67 @@ FROM dummy_rast; Description - Returns the hypothetical illumination of an elevation raster band using the azimuth, altitude, brightness, and elevation scale inputs. Utilizes map algebra and applies the hill shade equation to neighboring pixels. + Returns the hypothetical illumination of an elevation raster band using the azimuth, altitude, brightness, and scale inputs. Utilizes map algebra and applies the hill shade equation to neighboring pixels. Return pixel values are between 0 and 255. + + + azimuth is a value between 0 and 360 degrees measured clockwise from North. + + + + altitude is a value between 0 and 90 degrees where 0 degrees is at the horizon and 90 degrees is directly overhead. + + + + max_bright is a value between 0 and 255 with 0 as no brightness and 255 as max brightness. + + + + scale is the ratio of vertical units to horizontal. For Feet:LatLon use scale=370400, for Meters:LatLon use scale=111120. + If interpolate_nodata is TRUE, values for NODATA pixels from the input raster will be interpolated using before computing the hillshade illumination. - The hill shade equation is: max_bright * ( (cos(zenith)*cos(slope)) + (sin(zenith)*sin(slope)*cos(azimuth - aspect)) ). + + + For more information about Hillshade, please refer to How hillshade works. + + + Availability: 2.0.0 Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional interpolate_nodata function parameter + Changed: 2.1.0 In prior versions, azimuth and altitude were expressed in radians. Now, azimuth and altitude are expressed in degrees - Examples - coming soon -coming soon + Examples + +WITH foo AS ( + SELECT ST_SetValues( + ST_AddBand(ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0), 1, '32BF', 0, -9999), + 1, 1, 1, ARRAY[ + [1, 1, 1, 1, 1], + [1, 2, 2, 2, 1], + [1, 2, 3, 2, 1], + [1, 2, 2, 2, 1], + [1, 1, 1, 1, 1] + ]::double precision[][] + ) AS rast +) +SELECT + ST_DumpValues(ST_Hillshade(rast, 1, '32BF')) +FROM foo + st_dumpvalues + +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +----------------------------------------------------------------------- + (1,"{{NULL,NULL,NULL,NULL,NULL},{NULL,251.32763671875,220.749786376953,147.224319458008,NULL},{NULL,220.749786376953,180.312225341797,67.7497863769531,NULL},{NULL,147.224319458008 +,67.7497863769531,43.1210060119629,NULL},{NULL,NULL,NULL,NULL,NULL}}") +(1 row) + @@ -7743,15 +7788,16 @@ FROM dummy_rast; ST_Aspect - Returns the surface aspect of an elevation raster band. Useful for analyzing terrain. + Returns the aspect (in degrees by default) of an elevation raster band. Useful for analyzing terrain. raster ST_Aspect raster rast - integer band - text pixeltype + integer band=1 + text pixeltype=32BF + text measurement=DEGREES boolean interpolate_nodata=FALSE @@ -7760,46 +7806,62 @@ FROM dummy_rast; Description - Returns the surface aspect of an elevation raster band. Utilizes map algebra and applies the aspect equation to neighboring pixels. + Returns the aspect (in degrees by default) of an elevation raster band. Utilizes map algebra and applies the aspect equation to neighboring pixels. - If interpolate_nodata is TRUE, values for NODATA pixels from the input raster will be interpolated using before computing the surface aspect. + measurement indicates the units of the aspect. Possible values are: RADIANS, DEGREES (default). - Given the following representation of a 3x3 neighborhood of pixels: - - - - - - A - B - C - - - D - E - F - - - G - H - I - - - - - - The equation for the pixel aspect of cell E is: atan2((((G + 2H + I) - (A + 2B + C)) / 8), -(((C + 2F + I) - (A + 2D + G)) / 8)) + + When measurement = RADIANS, values are between 0 and 2π radians measured clockwise from North. + + + + When measurement = DEGREES, values are between 0 and 360 degrees measured clockwise from North. + + + + If slope of pixel is zero, aspect of pixel is -1. + + + + + For more information about Slope, Aspect and Hillshade, please refer to ESRI - How hillshade works and ERDAS Field Guide - Aspect Images. + + Availability: 2.0.0 Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional interpolate_nodata function parameter - + Changed: 2.1.0 In prior versions, return values were in radians. Now, return values default to degrees - Examples - coming soon - coming soon + Examples + +WITH foo AS ( + SELECT ST_SetValues( + ST_AddBand(ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0), 1, '32BF', 0, -9999), + 1, 1, 1, ARRAY[ + [1, 1, 1, 1, 1], + [1, 2, 2, 2, 1], + [1, 2, 3, 2, 1], + [1, 2, 2, 2, 1], + [1, 1, 1, 1, 1] + ]::double precision[][] + ) AS rast +) +SELECT + ST_DumpValues(ST_Aspect(rast, 1, '32BF')) +FROM foo + + st_dumpvalues + +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +---------------------------------- + (1,"{{315,341.565063476562,0,18.4349479675293,45},{288.434936523438,315,0,45,71.5650482177734},{270,270,-1,90,90},{251.565048217773,225,180,135,108.434951782227},{225,198.43495178 +2227,180,161.565048217773,135}}") +(1 row) + @@ -7816,15 +7878,17 @@ FROM dummy_rast; ST_Slope - Returns the surface slope of an elevation raster band. Useful for analyzing terrain. + Returns the slope (in degrees by default) of an elevation raster band. Useful for analyzing terrain. raster ST_Slope raster rast - integer band=1 + integer nband=1 text pixeltype=32BF + text measurement=DEGREES + double precision scale=1.0 boolean interpolate_nodata=FALSE @@ -7833,47 +7897,62 @@ FROM dummy_rast; Description - Returns the surface slope of an elevation raster band. Utilizes map algebra and applies the slope equation to neighboring pixels. + Returns the slope (in degrees by default) of an elevation raster band. Utilizes map algebra and applies the slope equation to neighboring pixels. + + + measurement indicates the units of the slope. Possible values are: RADIANS, DEGREES (default), PERCENT. + + + + scale is the ratio of vertical units to horizontal. For Feet:LatLon use scale=370400, for Meters:LatLon use scale=111120. + If interpolate_nodata is TRUE, values for NODATA pixels from the input raster will be interpolated using before computing the surface slope. - Given the following representation of a 3x3 neighborhood of pixels: - - - - - - A - B - C - - - D - E - F - - - G - H - I - - - - - - The equation for the pixel slope of cell E is: atan(sqrt(((c + 2f + i) - (a + 2d + g) / 8)^2 + (((g + 2h + i) - (a + 2b + c)) / 8) ^ 2)) + + + For more information about Slope, Aspect and Hillshade, please refer to ESRI - How hillshade works and ERDAS Field Guide - Slope Images. + + Availability: 2.0.0 - Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional interpolate_nodata function parameter + Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional measurement, scale, interpolate_nodata function parameters + Changed: 2.1.0 In prior versions, return values were in radians. Now, return values default to degrees - Examples - coming soon - coming soon + Examples + +WITH foo AS ( + SELECT ST_SetValues( + ST_AddBand(ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0), 1, '32BF', 0, -9999), + 1, 1, 1, ARRAY[ + [1, 1, 1, 1, 1], + [1, 2, 2, 2, 1], + [1, 2, 3, 2, 1], + [1, 2, 2, 2, 1], + [1, 1, 1, 1, 1] + ]::double precision[][] + ) AS rast +) +SELECT + ST_DumpValues(ST_Slope(rast, 1, '32BF')) +FROM foo + st_dumpvalues + +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ +--------------------------------------------------------------------- + (1,"{{10.0249881744385,21.5681285858154,26.5650520324707,21.5681285858154,10.0249881744385},{21.5681285858154,35.2643890380859,36.8698959350586,35.2643890380859,21.5681285858154}, +{26.5650520324707,36.8698959350586,0,36.8698959350586,26.5650520324707},{21.5681285858154,35.2643890380859,36.8698959350586,35.2643890380859,21.5681285858154},{10.0249881744385,21. +5681285858154,26.5650520324707,21.5681285858154,10.0249881744385}}") +(1 row) + + @@ -7885,7 +7964,7 @@ FROM dummy_rast; - + ST_Intersection diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 3257438d5..8ab85acd7 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -3424,19 +3424,32 @@ CREATE OR REPLACE FUNCTION st_mindist4ma(value double precision[][][], pos integ ----------------------------------------------------------------------- -- ST_Slope +-- Hill Shading and the Reflectance Map +-- by Berthold K.P. Horn +-- Proceedings of the IEEE, Vol 69, No. 1 ----------------------------------------------------------------------- CREATE OR REPLACE FUNCTION _st_slope4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL) RETURNS double precision AS $$ DECLARE - pwidth double precision; - pheight double precision; + x integer; + y integer; + z integer; + + _pixwidth double precision; + _pixheight double precision; + _width double precision; + _height double precision; + _measure text; + _scale double precision; + dz_dx double precision; dz_dy double precision; + slope double precision; + _value double precision[][][]; ndims int; - z int; BEGIN ndims := array_ndims(value); @@ -3449,6 +3462,12 @@ CREATE OR REPLACE FUNCTION _st_slope4ma(value double precision[][][], pos intege _value := value; END IF; + -- only use the first raster passed to this function + IF array_length(_value, 1) > 1 THEN + RAISE NOTICE 'Only using the values from the first raster'; + END IF; + z := array_lower(_value, 1); + IF ( array_lower(_value, 2) != 1 OR array_upper(_value, 2) != 3 OR array_lower(_value, 3) != 1 OR array_upper(_value, 3) != 3 @@ -3456,36 +3475,104 @@ CREATE OR REPLACE FUNCTION _st_slope4ma(value double precision[][][], pos intege RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; END IF; - IF array_length(userargs, 1) < 2 THEN - RAISE EXCEPTION 'At least two elements must be provided for the third parameter'; + IF array_length(userargs, 1) < 6 THEN + RAISE EXCEPTION 'At least six elements must be provided for the third parameter'; END IF; - -- only use the first raster passed to this function - IF array_length(_value, 1) > 1 THEN - RAISE NOTICE 'Only using the values from the first raster'; + _pixwidth := userargs[1]::double precision; + _pixheight := userargs[2]::double precision; + _width := userargs[3]::double precision; + _height := userargs[4]::double precision; + _measure := userargs[5]; + _scale := userargs[6]::double precision; + + /* ArcGIS returns values for edge pixels + -- check that pixel is not edge pixel + IF (pos[1][1] = 1 OR pos[1][2] = 1) OR (pos[1][1] = _width OR pos[1][2] = _height) THEN + RETURN NULL; + ELSEIF _value[z][2][2] IS NULL THEN + */ + -- check that center pixel isn't NODATA + IF _value[z][2][2] IS NULL THEN + RETURN NULL; + -- substitute center pixel for any neighbor pixels that are NODATA + ELSE + FOR y IN 1..3 LOOP + FOR x IN 1..3 LOOP + IF _value[z][y][x] IS NULL THEN + _value[z][y][x] = _value[z][2][2]; + END IF; + END LOOP; + END LOOP; END IF; - z := array_lower(_value, 1); - pwidth := userargs[1]::double precision; - pheight := userargs[2]::double precision; - dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); - dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); + dz_dy := ((_value[z][3][1] + _value[z][3][2] + _value[z][3][2] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][1][2] + _value[z][1][2] + _value[z][1][3])) / _pixheight; + dz_dx := ((_value[z][1][3] + _value[z][2][3] + _value[z][2][3] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][2][1] + _value[z][2][1] + _value[z][3][1])) / _pixwidth; - RETURN atan(sqrt(pow(dz_dx, 2.0) + pow(dz_dy, 2.0))); + slope := sqrt(dz_dx * dz_dx + dz_dy * dz_dy) / (8 * _scale); + + -- output depends on user preference + CASE substring(upper(trim(leading from _measure)) for 3) + -- percentages + WHEN 'PER' THEN + slope := 100.0 * slope; + -- radians + WHEN 'rad' THEN + slope := atan(slope); + -- degrees (default) + ELSE + slope := degrees(atan(slope)); + END CASE; + + RETURN slope; END; $$ LANGUAGE 'plpgsql' IMMUTABLE; -CREATE OR REPLACE FUNCTION st_slope(rast raster, band integer DEFAULT 1, pixeltype text DEFAULT '32BF', interpolate_nodata boolean DEFAULT FALSE) +CREATE OR REPLACE FUNCTION st_slope( + rast raster, nband integer DEFAULT 1, + pixeltype text DEFAULT '32BF', measurement text DEFAULT 'DEGREES', + scale double precision DEFAULT 1.0, interpolate_nodata boolean DEFAULT FALSE +) RETURNS raster AS $$ - SELECT - CASE - WHEN $4 IS FALSE THEN - ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], '_st_slope4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) - ELSE - ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_slope4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) - END - $$ LANGUAGE 'sql' IMMUTABLE; + DECLARE + _rast raster; + _nband integer; + _pixtype text; + _pixwidth double precision; + _pixheight double precision; + _width integer; + _height integer; + BEGIN + IF interpolate_nodata IS TRUE THEN + _rast := ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1); + _nband := 1; + _pixtype := NULL; + ELSE + _rast := rast; + _nband := nband; + _pixtype := pixeltype; + END IF; + + -- get properties + _pixwidth := ST_PixelWidth(_rast); + _pixheight := ST_PixelHeight(_rast); + SELECT width, height INTO _width, _height FROM ST_Metadata(_rast); + + RETURN ST_MapAlgebra( + ARRAY[ROW(_rast, _nband)]::rastbandarg[], + '_st_slope4ma(double precision[][][], integer[][], text[])'::regprocedure, + _pixtype, + 'FIRST', NULL, + 1, 1, + _pixwidth::text, _pixheight::text, + _width::text, _height::text, + $4::text, $5::text + ); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; ----------------------------------------------------------------------- -- ST_Aspect @@ -3494,15 +3581,21 @@ CREATE OR REPLACE FUNCTION _st_aspect4ma(value double precision[][][], pos integ RETURNS double precision AS $$ DECLARE - pwidth double precision; - pheight double precision; + x integer; + y integer; + z integer; + + _width double precision; + _height double precision; + _measure text; + dz_dx double precision; dz_dy double precision; aspect double precision; + halfpi double precision; _value double precision[][][]; ndims int; - z int; BEGIN ndims := array_ndims(value); -- add a third dimension if 2-dimension @@ -3521,8 +3614,8 @@ CREATE OR REPLACE FUNCTION _st_aspect4ma(value double precision[][][], pos integ RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; END IF; - IF array_length(userargs, 1) < 2 THEN - RAISE EXCEPTION 'At least two elements must be provided for the third parameter'; + IF array_length(userargs, 1) < 3 THEN + RAISE EXCEPTION 'At least three elements must be provided for the third parameter'; END IF; -- only use the first raster passed to this function @@ -3531,35 +3624,106 @@ CREATE OR REPLACE FUNCTION _st_aspect4ma(value double precision[][][], pos integ END IF; z := array_lower(_value, 1); - pwidth := userargs[1]::double precision; - pheight := userargs[2]::double precision; + _width := userargs[1]::double precision; + _height := userargs[2]::double precision; + _measure := userargs[3]; - dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); - dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); + /* ArcGIS returns values for edge pixels + -- check that pixel is not edge pixel + IF (pos[1][1] = 1 OR pos[1][2] = 1) OR (pos[1][1] = _width OR pos[1][2] = _height) THEN + RETURN NULL; + ELSEIF _value[z][2][2] IS NULL THEN + */ + -- check that center pixel isn't NODATA + IF _value[z][2][2] IS NULL THEN + RETURN NULL; + -- substitute center pixel for any neighbor pixels that are NODATA + ELSE + FOR y IN 1..3 LOOP + FOR x IN 1..3 LOOP + IF _value[z][y][x] IS NULL THEN + _value[z][y][x] = _value[z][2][2]; + END IF; + END LOOP; + END LOOP; + END IF; + + dz_dy := ((_value[z][3][1] + _value[z][3][2] + _value[z][3][2] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][1][2] + _value[z][1][2] + _value[z][1][3])); + dz_dx := ((_value[z][1][3] + _value[z][2][3] + _value[z][2][3] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][2][1] + _value[z][2][1] + _value[z][3][1])); + + -- aspect is flat IF abs(dz_dx) = 0::double precision AND abs(dz_dy) = 0::double precision THEN RETURN -1; END IF; + -- aspect is in radians aspect := atan2(dz_dy, -dz_dx); - IF aspect > (pi() / 2.0) THEN - RETURN (5.0 * pi() / 2.0) - aspect; + + -- north = 0, pi/2 = east, 3pi/2 = west + halfpi := pi() / 2.0; + IF aspect > halfpi THEN + aspect := (5.0 * halfpi) - aspect; ELSE - RETURN (pi() / 2.0) - aspect; + aspect := halfpi - aspect; + END IF; + + IF aspect = 2 * pi() THEN + aspect := 0.; END IF; + + -- output depends on user preference + CASE substring(upper(trim(leading from _measure)) for 3) + -- radians + WHEN 'rad' THEN + RETURN aspect; + -- degrees (default) + ELSE + RETURN degrees(aspect); + END CASE; + END; $$ LANGUAGE 'plpgsql' IMMUTABLE; -CREATE OR REPLACE FUNCTION st_aspect(rast raster, band integer, pixeltype text, interpolate_nodata boolean DEFAULT FALSE) +CREATE OR REPLACE FUNCTION st_aspect( + rast raster, nband integer DEFAULT 1, + pixeltype text DEFAULT '32BF', measurement text DEFAULT 'DEGREES', + interpolate_nodata boolean DEFAULT FALSE +) RETURNS raster AS $$ - SELECT - CASE - WHEN $4 IS FALSE THEN - ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], '_st_aspect4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) - ELSE - ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_aspect4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text) - END - $$ LANGUAGE 'sql' IMMUTABLE; + DECLARE + _rast raster; + _nband integer; + _pixtype text; + _width integer; + _height integer; + BEGIN + IF interpolate_nodata IS TRUE THEN + _rast := ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1); + _nband := 1; + _pixtype := NULL; + ELSE + _rast := rast; + _nband := nband; + _pixtype := pixeltype; + END IF; + + -- get properties + SELECT width, height INTO _width, _height FROM ST_Metadata(_rast); + + RETURN ST_MapAlgebra( + ARRAY[ROW(_rast, _nband)]::rastbandarg[], + '_st_aspect4ma(double precision[][][], integer[][], text[])'::regprocedure, + _pixtype, + 'FIRST', NULL, + 1, 1, + _width::text, _height::text, + $4::text + ); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; ----------------------------------------------------------------------- -- ST_HillShade @@ -3568,16 +3732,22 @@ CREATE OR REPLACE FUNCTION _st_hillshade4ma(value double precision[][][], pos in RETURNS double precision AS $$ DECLARE - pwidth double precision; - pheight double precision; + _pixwidth double precision; + _pixheight double precision; + _width double precision; + _height double precision; + _azimuth double precision; + _altitude double precision; + _bright double precision; + _scale double precision; + dz_dx double precision; dz_dy double precision; - zenith double precision; azimuth double precision; + zenith double precision; slope double precision; aspect double precision; - max_bright double precision; - elevation_scale double precision; + shade double precision; _value double precision[][][]; ndims int; @@ -3600,8 +3770,8 @@ CREATE OR REPLACE FUNCTION _st_hillshade4ma(value double precision[][][], pos in RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1'; END IF; - IF array_length(userargs, 1) < 6 THEN - RAISE EXCEPTION 'At least six elements must be provided for the third parameter'; + IF array_length(userargs, 1) < 8 THEN + RAISE EXCEPTION 'At least eight elements must be provided for the third parameter'; END IF; -- only use the first raster passed to this function @@ -3610,50 +3780,130 @@ CREATE OR REPLACE FUNCTION _st_hillshade4ma(value double precision[][][], pos in END IF; z := array_lower(_value, 1); - pwidth := userargs[1]::double precision; - pheight := userargs[2]::double precision; + _pixwidth := userargs[1]::double precision; + _pixheight := userargs[2]::double precision; + _width := userargs[3]::double precision; + _height := userargs[4]::double precision; + _azimuth := userargs[5]::double precision; + _altitude := userargs[6]::double precision; + _bright := userargs[7]::double precision; + _scale := userargs[8]::double precision; + + -- check that pixel is not edge pixel + IF (pos[1][1] = 1 OR pos[1][2] = 1) OR (pos[1][1] = _width OR pos[1][2] = _height) THEN + RETURN NULL; + END IF; + + -- clamp azimuth + IF _azimuth < 0. THEN + RAISE NOTICE 'Clamping provided azimuth value % to 0', _azimuth; + _azimuth := 0.; + ELSEIF _azimuth >= 360. THEN + RAISE NOTICE 'Converting provided azimuth value % to be between 0 and 360', _azimuth; + _azimuth := _azimuth - (360. * floor(_azimuth / 360.)); + END IF; + azimuth := 360. - _azimuth + 90.; + IF azimuth >= 360. THEN + azimuth := azimuth - 360.; + END IF; + azimuth := radians(azimuth); + --RAISE NOTICE 'azimuth = %', azimuth; - azimuth := (5.0 * pi() / 2.0) - userargs[3]::double precision; - zenith := (pi() / 2.0) - userargs[4]::double precision; - dz_dy := ((_value[z][3][1] + 2.0 * _value[z][3][2] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][1][2] + _value[z][1][3])) / (8.0 * pwidth); - dz_dx := ((_value[z][1][3] + 2.0 * _value[z][2][3] + _value[z][3][3]) - (_value[z][1][1] + 2.0 * _value[z][2][1] + _value[z][3][1])) / (8.0 * pheight); - elevation_scale := userargs[6]::double precision; - slope := atan(sqrt(elevation_scale * pow(dz_dx, 2.0) + pow(dz_dy, 2.0))); + -- clamp altitude + IF _altitude < 0. THEN + RAISE NOTICE 'Clamping provided altitude value % to 0', _altitude; + _altitude := 0.; + ELSEIF _altitude > 90. THEN + RAISE NOTICE 'Clamping provided altitude value % to 90', _altitude; + _altitude := 90.; + END IF; + zenith := radians(90. - _altitude); + --RAISE NOTICE 'zenith = %', zenith; - -- handle special case of 0, 0 - IF abs(dz_dy) = 0::double precision AND abs(dz_dy) = 0::double precision THEN - -- set to pi as that is the expected PostgreSQL answer in Linux - aspect := pi(); - ELSE - aspect := atan2(dz_dy, -dz_dx); + -- clamp bright + IF _bright < 0. THEN + RAISE NOTICE 'Clamping provided bright value % to 0', _bright; + _bright := 0.; + ELSEIF _bright > 255. THEN + RAISE NOTICE 'Clamping provided bright value % to 255', _bright; + _bright := 255.; END IF; - max_bright := userargs[5]::double precision; - IF aspect < 0 THEN - aspect := aspect + (2.0 * pi()); + dz_dy := ((_value[z][3][1] + _value[z][3][2] + _value[z][3][2] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][1][2] + _value[z][1][2] + _value[z][1][3])) / (8 * _pixheight); + dz_dx := ((_value[z][1][3] + _value[z][2][3] + _value[z][2][3] + _value[z][3][3]) - + (_value[z][1][1] + _value[z][2][1] + _value[z][2][1] + _value[z][3][1])) / (8 * _pixwidth); + + slope := atan(sqrt(dz_dx * dz_dx + dz_dy * dz_dy) / _scale); + + IF dz_dx != 0. THEN + aspect := atan2(dz_dy, -dz_dx); + + IF aspect < 0. THEN + aspect := aspect + (2.0 * pi()); + END IF; + ELSE + IF dz_dy > 0. THEN + aspect := pi() / 2.; + ELSEIF dz_dy < 0. THEN + aspect := (2. * pi()) - (pi() / 2.); + -- set to pi as that is the expected PostgreSQL answer in Linux + ELSE + aspect := pi(); + END IF; END IF; - RETURN max_bright * ((cos(zenith) * cos(slope)) + (sin(zenith) * sin(slope) * cos(azimuth - aspect))); + shade := _bright * ((cos(zenith) * cos(slope)) + (sin(zenith) * sin(slope) * cos(azimuth - aspect))); + RETURN shade; END; $$ LANGUAGE 'plpgsql' IMMUTABLE; CREATE OR REPLACE FUNCTION st_hillshade( - rast raster, band integer, - pixeltype text, - azimuth double precision, altitude double precision, - max_bright double precision DEFAULT 255.0, elevation_scale double precision DEFAULT 1.0, + rast raster, nband integer DEFAULT 1, + pixeltype text DEFAULT '32BF', + azimuth double precision DEFAULT 315.0, altitude double precision DEFAULT 45.0, + max_bright double precision DEFAULT 255.0, scale double precision DEFAULT 1.0, interpolate_nodata boolean DEFAULT FALSE ) RETURNS RASTER AS $$ - SELECT - CASE - WHEN $8 IS FALSE THEN - ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], '_st_hillshade4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text, $4::text, $5::text, $6::text, $7::text) - ELSE - ST_MapAlgebra(ARRAY[ROW(ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1), 1)]::rastbandarg[], '_st_hillshade4ma(double precision[][][], integer[][], text[])'::regprocedure, NULL, 'FIRST', NULL, 1, 1, st_pixelwidth($1)::text, st_pixelheight($1)::text, $4::text, $5::text, $6::text, $7::text) - END - $$ LANGUAGE 'sql' IMMUTABLE; + DECLARE + _rast raster; + _nband integer; + _pixtype text; + _pixwidth double precision; + _pixheight double precision; + _width integer; + _height integer; + BEGIN + IF interpolate_nodata IS TRUE THEN + _rast := ST_MapAlgebra(ARRAY[ROW($1, $2)]::rastbandarg[], 'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure, $3, 'FIRST', NULL, 1, 1); + _nband := 1; + _pixtype := NULL; + ELSE + _rast := rast; + _nband := nband; + _pixtype := pixeltype; + END IF; + + -- get properties + _pixwidth := ST_PixelWidth(_rast); + _pixheight := ST_PixelHeight(_rast); + SELECT width, height, scalex INTO _width, _height FROM ST_Metadata(_rast); + + RETURN ST_MapAlgebra( + ARRAY[ROW(_rast, _nband)]::rastbandarg[], + '_st_hillshade4ma(double precision[][][], integer[][], text[])'::regprocedure, + _pixtype, + 'FIRST', NULL, + 1, 1, + _pixwidth::text, _pixheight::text, + _width::text, _height::text, + $4::text, $5::text, + $6::text, $7::text + ); + END; + $$ LANGUAGE 'plpgsql' IMMUTABLE; ----------------------------------------------------------------------- -- Get information about the raster diff --git a/raster/rt_pg/rtpostgis_drop.sql.in.c b/raster/rt_pg/rtpostgis_drop.sql.in.c index eda7aa374..0a0373cd1 100644 --- a/raster/rt_pg/rtpostgis_drop.sql.in.c +++ b/raster/rt_pg/rtpostgis_drop.sql.in.c @@ -433,8 +433,12 @@ DROP FUNCTION IF EXISTS _st_contains(geometry, raster, integer); DROP FUNCTION IF EXISTS st_addband(raster, raster[], integer); -- function signatures changed +DROP FUNCTION IF EXISTS st_slope(raster, integer, text, boolean); DROP FUNCTION IF EXISTS st_slope(raster, integer, text); +DROP FUNCTION IF EXISTS st_aspect(raster, integer, text, boolean); DROP FUNCTION IF EXISTS st_aspect(raster, integer, text); +DROP FUNCTION IF EXISTS st_hillshade(raster, integer, text, double precision, double precision, double precision, double precision, boolean); +DROP FUNCTION IF EXISTS st_hillshade(raster, integer, text, float, float, float, float, boolean); DROP FUNCTION IF EXISTS st_hillshade(raster, integer, text, float, float, float, float); -- function no longer exists diff --git a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql index f78d9a76f..eb51f1899 100644 --- a/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql +++ b/raster/test/regress/rt_mapalgebrafctngb_userfunc.sql @@ -133,17 +133,17 @@ SELECT -- test st_slope, corner spike SELECT ST_Value(rast, 2, 2) = 1, - round(degrees(ST_Value( + round(ST_Value( ST_Slope(rast, 1, NULL), 2, 2 - ))*100000) = 5784902 + )*100000) = 5784902 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( + round(ST_Value( ST_Slope(rast, 1, NULL), 2, 2 - ))*100000) = 5784902 + )*100000) = 5784902 FROM ST_SetValue(ST_TestRasterNgb(3, 3, 10), 1, 1, 1) AS rast; -- test st_aspect, flat plane @@ -157,47 +157,47 @@ SELECT -- test st_aspect, North SELECT ST_Value(rast, 2, 2) = 1, - round(degrees(ST_Value( + round(ST_Value( ST_Aspect(rast, 1, NULL), 2, 2 - ))*10000) = 0 + )*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( + round(ST_Value( ST_Aspect(rast, 1, NULL), 2, 2 - ))*10000) = 1800000 + )*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( + round(ST_Value( ST_Aspect(rast, 1, NULL), 2, 2 - ))*10000) = 900000 - FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 3, 2, 0) AS rast; + )*10000) = 900000 + FROM ST_SetValue(ST_SetScale(ST_TestRasterNgb(3, 3, 1), 1, -1), 3, 2, 0) AS rast; -- test st_aspect, West SELECT ST_Value(rast, 2, 2) = 1, - round(degrees(ST_Value( + round(ST_Value( ST_Aspect(rast, 1, NULL), 2, 2 - ))*10000) = 2700000 - FROM ST_SetValue(ST_TestRasterNgb(3, 3, 1), 1, 2, 0) AS rast; + )*10000) = 2700000 + FROM ST_SetValue(ST_SetScale(ST_TestRasterNgb(3, 3, 1), 1, -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 + ST_Hillshade(rast, 1, NULL, 0.0, 45., 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 + ST_Hillshade(rast, 1, NULL, 315.0, 45., 255, 1.0), 2, 2 )*10000) = 1540287 FROM ST_SetValue( ST_SetValue( @@ -222,7 +222,7 @@ SELECT -- 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 + ST_Hillshade(rast, 1, NULL, 315.0, 45.), 2, 2 )*10000) = 1540287 FROM ST_SetValue( ST_SetValue( -- 2.50.1