<refentry id="RT_ST_HillShade">
<refnamediv>
<refname>ST_HillShade</refname>
- <refpurpose>Returns the hypothetical illumination of an elevation raster band using provided azimuth, altitude, brightness and elevation scale inputs. Useful for visualizing terrain.</refpurpose>
+ <refpurpose>Returns the hypothetical illumination of an elevation raster band using provided azimuth, altitude, brightness and scale inputs.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
<paramdef><type>integer </type> <parameter>band</parameter></paramdef>
<paramdef><type>text </type> <parameter>pixeltype</parameter></paramdef>
- <paramdef><type>double precision </type> <parameter>azimuth</parameter></paramdef>
- <paramdef><type>double precision </type> <parameter>altitude</parameter></paramdef>
+ <paramdef choice="opt"><type>double precision </type> <parameter>azimuth=315</parameter></paramdef>
+ <paramdef choice="opt"><type>double precision </type> <parameter>altitude=45</parameter></paramdef>
<paramdef choice="opt"><type>double precision </type> <parameter>max_bright=255</parameter></paramdef>
- <paramdef choice="opt"><type>double precision </type> <parameter>elevation_scale=1</parameter></paramdef>
+ <paramdef choice="opt"><type>double precision </type> <parameter>scale=1.0</parameter></paramdef>
<paramdef choice="opt"><type>boolean </type> <parameter>interpolate_nodata=FALSE</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<refsection>
<title>Description</title>
- <para>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.</para>
+ <para>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.</para>
+
+ <para>
+ <varname>azimuth</varname> is a value between 0 and 360 degrees measured clockwise from North.
+ </para>
+
+ <para>
+ <varname>altitude</varname> is a value between 0 and 90 degrees where 0 degrees is at the horizon and 90 degrees is directly overhead.
+ </para>
+
+ <para>
+ <varname>max_bright</varname> is a value between 0 and 255 with 0 as no brightness and 255 as max brightness.
+ </para>
+
+ <para>
+ <varname>scale</varname> is the ratio of vertical units to horizontal. For Feet:LatLon use scale=370400, for Meters:LatLon use scale=111120.
+ </para>
<para>
If <varname>interpolate_nodata</varname> is TRUE, values for NODATA pixels from the input raster will be interpolated using <xref linkend="RT_ST_InvDistWeight4ma" /> before computing the hillshade illumination.
</para>
- <para>The hill shade equation is: <programlisting>max_bright * ( (cos(zenith)*cos(slope)) + (sin(zenith)*sin(slope)*cos(azimuth - aspect)) )</programlisting>.</para>
+ <note>
+ <para>
+ For more information about Hillshade, please refer to <ulink url="http://edndoc.esri.com/arcobjects/9.2/net/shared/geoprocessing/spatial_analyst_tools/how_hillshade_works.htm">How hillshade works</ulink>.
+ </para>
+ </note>
+
<para>Availability: 2.0.0 </para>
<para>Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional <varname>interpolate_nodata</varname> function parameter</para>
+ <para>Changed: 2.1.0 In prior versions, azimuth and altitude were expressed in radians. Now, azimuth and altitude are expressed in degrees</para>
</refsection>
<refsection>
- <title>Examples - coming soon</title>
-<programlisting>coming soon</programlisting>
+ <title>Examples</title>
+ <programlisting>
+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)
+ </programlisting>
</refsection>
<refsection>
<refentry id="RT_ST_Aspect">
<refnamediv>
<refname>ST_Aspect</refname>
- <refpurpose>Returns the surface aspect of an elevation raster band. Useful for analyzing terrain.</refpurpose>
+ <refpurpose>Returns the aspect (in degrees by default) of an elevation raster band. Useful for analyzing terrain.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
<funcdef>raster <function>ST_Aspect</function></funcdef>
<paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
- <paramdef><type>integer </type> <parameter>band</parameter></paramdef>
- <paramdef><type>text </type> <parameter>pixeltype</parameter></paramdef>
+ <paramdef choice="opt"><type>integer </type> <parameter>band=1</parameter></paramdef>
+ <paramdef choice="opt"><type>text </type> <parameter>pixeltype=32BF</parameter></paramdef>
+ <paramdef choice="opt"><type>text </type> <parameter>measurement=DEGREES</parameter></paramdef>
<paramdef choice="opt"><type>boolean </type> <parameter>interpolate_nodata=FALSE</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<refsection>
<title>Description</title>
- <para>Returns the surface aspect of an elevation raster band. Utilizes map algebra and applies the aspect equation to neighboring pixels.</para>
+ <para>Returns the aspect (in degrees by default) of an elevation raster band. Utilizes map algebra and applies the aspect equation to neighboring pixels.</para>
<para>
- If <varname>interpolate_nodata</varname> is TRUE, values for NODATA pixels from the input raster will be interpolated using <xref linkend="RT_ST_InvDistWeight4ma" /> before computing the surface aspect.
+ <varname>measurement</varname> indicates the units of the aspect. Possible values are: RADIANS, DEGREES (default).
</para>
- <para>Given the following representation of a 3x3 neighborhood of pixels:</para>
-
- <informaltable rowsep="1" frame="all">
- <tgroup cols="3">
- <tbody>
- <row>
- <entry>A</entry>
- <entry>B</entry>
- <entry>C</entry>
- </row>
- <row>
- <entry>D</entry>
- <entry>E</entry>
- <entry>F</entry>
- </row>
- <row>
- <entry>G</entry>
- <entry>H</entry>
- <entry>I</entry>
- </row>
- </tbody>
- </tgroup>
- </informaltable>
-
- <para>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))</para>
+ <para>
+ When <varname>measurement</varname> = RADIANS, values are between 0 and 2π radians measured clockwise from North.
+ </para>
+
+ <para>
+ When <varname>measurement</varname> = DEGREES, values are between 0 and 360 degrees measured clockwise from North.
+ </para>
+
+ <para>
+ If slope of pixel is zero, aspect of pixel is -1.
+ </para>
+
+ <note>
+ <para>
+ For more information about Slope, Aspect and Hillshade, please refer to <ulink url="http://edndoc.esri.com/arcobjects/9.2/net/shared/geoprocessing/spatial_analyst_tools/how_hillshade_works.htm">ESRI - How hillshade works</ulink> and <ulink url="http://geospatial.intergraph.com/fieldguide/wwhelp/wwhimpl/common/html/wwhelp.htm?context=FieldGuide&file=Aspect_Images.html">ERDAS Field Guide - Aspect Images</ulink>.
+ </para>
+ </note>
<para>Availability: 2.0.0 </para>
<para>Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional <varname>interpolate_nodata</varname> function parameter</para>
-
+ <para>Changed: 2.1.0 In prior versions, return values were in radians. Now, return values default to degrees</para>
</refsection>
<refsection>
- <title>Examples - coming soon</title>
- <programlisting>coming soon</programlisting>
+ <title>Examples</title>
+ <programlisting>
+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)
+ </programlisting>
</refsection>
<refentry id="RT_ST_Slope">
<refnamediv>
<refname>ST_Slope</refname>
- <refpurpose>Returns the surface slope of an elevation raster band. Useful for analyzing terrain.</refpurpose>
+ <refpurpose>Returns the slope (in degrees by default) of an elevation raster band. Useful for analyzing terrain.</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
<funcdef>raster <function>ST_Slope</function></funcdef>
<paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
- <paramdef choice="opt"><type>integer </type> <parameter>band=1</parameter></paramdef>
+ <paramdef choice="opt"><type>integer </type> <parameter>nband=1</parameter></paramdef>
<paramdef choice="opt"><type>text </type> <parameter>pixeltype=32BF</parameter></paramdef>
+ <paramdef choice="opt"><type>text </type> <parameter>measurement=DEGREES</parameter></paramdef>
+ <paramdef choice="opt"><type>double precision </type> <parameter>scale=1.0</parameter></paramdef>
<paramdef choice="opt"><type>boolean </type> <parameter>interpolate_nodata=FALSE</parameter></paramdef>
</funcprototype>
</funcsynopsis>
<refsection>
<title>Description</title>
- <para>Returns the surface slope of an elevation raster band. Utilizes map algebra and applies the slope equation to neighboring pixels.</para>
+ <para>Returns the slope (in degrees by default) of an elevation raster band. Utilizes map algebra and applies the slope equation to neighboring pixels.</para>
+
+ <para>
+ <varname>measurement</varname> indicates the units of the slope. Possible values are: RADIANS, DEGREES (default), PERCENT.
+ </para>
+
+ <para>
+ <varname>scale</varname> is the ratio of vertical units to horizontal. For Feet:LatLon use scale=370400, for Meters:LatLon use scale=111120.
+ </para>
<para>
If <varname>interpolate_nodata</varname> is TRUE, values for NODATA pixels from the input raster will be interpolated using <xref linkend="RT_ST_InvDistWeight4ma" /> before computing the surface slope.
</para>
- <para>Given the following representation of a 3x3 neighborhood of pixels:</para>
-
- <informaltable rowsep="1" frame="all">
- <tgroup cols="3">
- <tbody>
- <row>
- <entry>A</entry>
- <entry>B</entry>
- <entry>C</entry>
- </row>
- <row>
- <entry>D</entry>
- <entry>E</entry>
- <entry>F</entry>
- </row>
- <row>
- <entry>G</entry>
- <entry>H</entry>
- <entry>I</entry>
- </row>
- </tbody>
- </tgroup>
- </informaltable>
-
- <para>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))</para>
+ <note>
+ <para>
+ For more information about Slope, Aspect and Hillshade, please refer to <ulink url="http://edndoc.esri.com/arcobjects/9.2/net/shared/geoprocessing/spatial_analyst_tools/how_hillshade_works.htm">ESRI - How hillshade works</ulink> and <ulink url="http://geospatial.intergraph.com/fieldguide/wwhelp/wwhimpl/common/html/wwhelp.htm?context=FieldGuide&file=Slope_Images.html">ERDAS Field Guide - Slope Images</ulink>.
+ </para>
+ </note>
<para>Availability: 2.0.0 </para>
- <para>Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional <varname>interpolate_nodata</varname> function parameter</para>
+ <para>Enhanced: 2.1.0 Uses ST_MapAlgebra() and added optional <varname>measurement</varname>, <varname>scale</varname>, <varname>interpolate_nodata</varname> function parameters</para>
+ <para>Changed: 2.1.0 In prior versions, return values were in radians. Now, return values default to degrees</para>
</refsection>
<refsection>
- <title>Examples - coming soon</title>
- <programlisting>coming soon</programlisting>
+ <title>Examples</title>
+ <programlisting>
+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)
+
+ </programlisting>
</refsection>
<refsection>
</para>
</refsection>
</refentry>
-
+
<refentry id="RT_ST_Intersection">
<refnamediv>
<refname>ST_Intersection</refname>
-----------------------------------------------------------------------
-- 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);
_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
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
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
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
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
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;
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
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
-- 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
-- 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(
-- 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(
projects / postgis / commitdiff