Added ST_Roughness(raster). Ticket #2302. Thanks to Nathaniel Clay for writing the...

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

diff --git a/NEWS b/NEWS
index 47bddd3666b591ada6d92f9efb608a56d9bcc2b9..c395ba4612ff20602616fd8ed69a698ef5a66a0c 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -84,6 +84,7 @@ PostGIS 2.1.0
   - #2276, ST_AddBand(raster) variant for out-db bands
   - #2280, ST_Summary(raster)
   - #2163, ST_TPI for raster (Nathaniel Clay)
+  - #2302, ST_Roughness for raster (Nathaniel Clay)
 
 * Enhancements *
 

diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml
index 659a307c3277b44bbfeeb7ccb627d0d225a0366a..51cdd52b9debcaa7da240c9bfb061046cde29d49 100644 (file)
--- a/doc/reference_raster.xml
+++ b/doc/reference_raster.xml
@@ -8329,6 +8329,9 @@ GROUP BY t1.rast;
                                <title>See Also</title>
                                <para>
                                        <xref linkend="RT_ST_MapAlgebra" />, 
+                                       <xref linkend="RT_ST_TRI" />, 
+                                       <xref linkend="RT_ST_TPI" />, 
+                                       <xref linkend="RT_ST_Roughness" />, 
                                        <xref linkend="RT_ST_Aspect" />, 
                                        <xref linkend="RT_ST_Slope" />
                                </para>
@@ -8464,6 +8467,9 @@ GROUP BY t1.rast;
                                <title>See Also</title>
                                <para>
                                        <xref linkend="RT_ST_MapAlgebra" />, 
+                                       <xref linkend="RT_ST_TRI" />, 
+                                       <xref linkend="RT_ST_TPI" />, 
+                                       <xref linkend="RT_ST_Roughness" />, 
                                        <xref linkend="RT_ST_HillShade" />, 
                                        <xref linkend="RT_ST_Slope" />
                                </para>
@@ -8599,6 +8605,9 @@ GROUP BY t1.rast;
                                <title>See Also</title>
                                <para>
                                        <xref linkend="RT_ST_MapAlgebra" />, 
+                                       <xref linkend="RT_ST_TRI" />, 
+                                       <xref linkend="RT_ST_TPI" />, 
+                                       <xref linkend="RT_ST_Roughness" />, 
                                        <xref linkend="RT_ST_HillShade" />, 
                                        <xref linkend="RT_ST_Aspect" />
                                </para>
@@ -8617,7 +8626,7 @@ GROUP BY t1.rast;
                                                <funcdef>raster <function>ST_TPI</function></funcdef>
                                                <paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
                                                <paramdef><type>integer </type> <parameter>nband</parameter></paramdef>
-                                               <paramdef choice="opt"><type>raster </type> <parameter>customextent</parameter></paramdef>
+                                               <paramdef><type>raster </type> <parameter>customextent</parameter></paramdef>
                                                <paramdef choice="opt"><type>text </type> <parameter>pixeltype="32BF"</parameter> </paramdef>
                                                <paramdef choice="opt"><type>boolean </type> <parameter> interpolate_nodata=FALSE </parameter> </paramdef>
                                        </funcprototype>
@@ -8647,6 +8656,53 @@ GROUP BY t1.rast;
                                <para>
                                        <xref linkend="RT_ST_MapAlgebra" />, 
                                        <xref linkend="RT_ST_TRI" />, 
+                                       <xref linkend="RT_ST_Roughness" />, 
+                                       <xref linkend="RT_ST_Slope" />, 
+                                       <xref linkend="RT_ST_HillShade" />, 
+                                       <xref linkend="RT_ST_Aspect" />
+                               </para>
+                       </refsection>
+               </refentry>
+
+               <refentry id="RT_ST_Roughness">
+                       <refnamediv>
+                               <refname>ST_Roughness</refname>
+                               <refpurpose>Returns a raster with the calculated "roughness" of a DEM.</refpurpose>
+                       </refnamediv>
+
+                       <refsynopsisdiv>
+                               <funcsynopsis>
+                                       <funcprototype>
+                                               <funcdef>raster <function>ST_Roughness</function></funcdef>
+                                               <paramdef><type>raster </type> <parameter>rast</parameter></paramdef>
+                                               <paramdef><type>integer </type> <parameter>nband</parameter></paramdef>
+                                               <paramdef><type>raster </type> <parameter>customextent</parameter></paramdef>
+                                               <paramdef choice="opt"><type>text </type> <parameter>pixeltype="32BF"</parameter> </paramdef>
+                                               <paramdef choice="opt"><type>boolean </type> <parameter> interpolate_nodata=FALSE </parameter> </paramdef>
+                                       </funcprototype>
+                               </funcsynopsis>
+                       </refsynopsisdiv>
+
+                       <refsection>
+                               <title>Description</title>
+                   <para>Calculates the "roughness" of a DEM, by subtracting the maximum from the minimum for a given area.</para>
+                               <para>Availability: 2.1</para>
+                       </refsection>
+
+                       <refsection>
+                               <title>Examples</title>
+                               <programlisting>
+-- needs examples
+                               </programlisting>
+                 </refsection>
+
+                 <refsection>
+                   <title>See Also</title>
+                               <para>
+                                       <xref linkend="RT_ST_MapAlgebra" />, 
+                                       <xref linkend="RT_ST_TRI" />, 
+                                       <xref linkend="RT_ST_TPI" />, 
+                                       <xref linkend="RT_ST_Slope" />, 
                                        <xref linkend="RT_ST_HillShade" />, 
                                        <xref linkend="RT_ST_Aspect" />
                                </para>

diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in
index f6a3d008c1995f5bbc1a379d6a69c5fcf23293a9..8c13476df08b644a21da9f9058787eb4247733a0 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in
+++ b/raster/rt_pg/rtpostgis.sql.in
@@ -4118,6 +4118,129 @@ CREATE OR REPLACE FUNCTION st_tpi(
        END;
        $$ LANGUAGE 'plpgsql' IMMUTABLE;
 
+-----------------------------------------------------------------------
+-- ST_Roughness
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION _st_roughness4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL)
+       RETURNS double precision
+       AS $$
+       DECLARE
+               x integer;
+               y integer;
+               z integer;
+
+               minimum double precision;
+               maximum double precision;
+
+               _value double precision[][][];
+               ndims int;
+       BEGIN
+
+               ndims := array_ndims(value);
+               -- add a third dimension if 2-dimension
+               IF ndims = 2 THEN
+                       _value := _st_convertarray4ma(value);
+               ELSEIF ndims != 3 THEN
+                       RAISE EXCEPTION 'First parameter of function must be a 3-dimension array';
+               ELSE
+                       _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
+               ) THEN
+                       RAISE EXCEPTION 'First parameter of function must be a 1x3x3 array with each of the lower bounds starting from 1';
+               END IF;
+
+               -- 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;
+
+               minimum := _value[z][1][1];
+               maximum := _value[z][1][1];
+
+               FOR Y IN 1..3 LOOP
+                   FOR X IN 1..3 LOOP
+                        IF _value[z][y][x] < minimum THEN
+                           minimum := _value[z][y][x];
+                        ELSIF _value[z][y][x] > maximum THEN
+                           maximum := _value[z][y][x];
+                        END IF;
+                   END LOOP;
+               END LOOP;
+
+               RETURN maximum - minimum;
+       END;
+       $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_roughness(
+       rast raster, nband integer,
+       customextent raster,
+       pixeltype text DEFAULT '32BF', interpolate_nodata boolean DEFAULT FALSE
+)
+       RETURNS raster
+       AS $$
+       DECLARE
+               _rast raster;
+               _nband integer;
+               _pixtype text;
+               _pixwidth double precision;
+               _pixheight double precision;
+               _width integer;
+               _height integer;
+               _customextent raster;
+               _extenttype text;
+       BEGIN
+               _customextent := customextent;
+               IF _customextent IS NULL THEN
+                       _extenttype := 'FIRST';
+               ELSE
+                       _extenttype := 'CUSTOM';
+               END IF;
+
+               IF interpolate_nodata IS TRUE THEN
+                       _rast := ST_MapAlgebra(
+                               ARRAY[ROW(rast, nband)]::rastbandarg[],
+                               'st_invdistweight4ma(double precision[][][], integer[][], text[])'::regprocedure,
+                               pixeltype,
+                               'FIRST', NULL,
+                               1, 1
+                       );
+                       _nband := 1;
+                       _pixtype := NULL;
+               ELSE
+                       _rast := rast;
+                       _nband := nband;
+                       _pixtype := pixeltype;
+               END IF;
+
+               RETURN ST_MapAlgebra(
+                       ARRAY[ROW(_rast, _nband)]::rastbandarg[],
+                       '_st_roughness4ma(double precision[][][], integer[][], text[])'::regprocedure,
+                       _pixtype,
+                       _extenttype, _customextent,
+                       1, 1);
+       END;
+       $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
 -----------------------------------------------------------------------
 -- Get information about the raster
 -----------------------------------------------------------------------

diff --git a/raster/test/regress/rt_elevation_functions.sql b/raster/test/regress/rt_elevation_functions.sql
index 73f0740d8c62f4c94a0ea89407546f5b5cf19c10..d1cee1c9ab785b53cb3185b7f3ccbbba4ff89bb3 100644 (file)
--- a/raster/test/regress/rt_elevation_functions.sql
+++ b/raster/test/regress/rt_elevation_functions.sql
@@ -358,7 +358,14 @@ INSERT INTO raster_value_arrays VALUES
 SELECT
        id,
        val,
-       _st_tpi4ma(val,NULL,NULL)::numeric(16, 10) as tpi
+       round(_st_tpi4ma(val,NULL,NULL)::numeric, 6) as tpi
+FROM raster_value_arrays
+ORDER BY id;
+
+SELECT
+       id,
+       val,
+       round(_st_roughness4ma(val,NULL,NULL)::numeric, 6) as roughness
 FROM raster_value_arrays
 ORDER BY id;
 

diff --git a/raster/test/regress/rt_elevation_functions_expected b/raster/test/regress/rt_elevation_functions_expected
index 17221caf066f4a608456f6ec11ea19fdceef2123..e5f80b3282f18e6954ae4d92cc605b78b58971d7 100644 (file)
--- a/raster/test/regress/rt_elevation_functions_expected
+++ b/raster/test/regress/rt_elevation_functions_expected
@@ -485,14 +485,25 @@ DO
 9|slope|1|3|26.565052
 9|slope|2|3|21.568129
 9|slope|3|3|10.024988
-1|{{{1,2,3},{4,5,6},{7,8,9}}}|0.0000000000
-2|{{{15,14,13},{12,11,10},{9,8,7}}}|0.0000000000
-3|{{{1,3,5},{7,9,11},{13,15,17}}}|0.0000000000
+1|{{{1,2,3},{4,5,6},{7,8,9}}}|0.000000
+2|{{{15,14,13},{12,11,10},{9,8,7}}}|0.000000
+3|{{{1,3,5},{7,9,11},{13,15,17}}}|0.000000
 4|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}|
-5|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|0.0000000000
+5|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|0.000000
 6|{{{NULL,2,NULL},{4,NULL,6},{NULL,8,NULL}}}|
 7|{{{1,3.1,5.2},{NULL,NULL,11.5},{13.6,NULL,NULL}}}|
-8|{{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}}|-0.3925000000
+8|{{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}}|-0.392500
 9|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,9}}}|
-10|{{{3.14,15.14,27.14},{39.14,51.14,63.14},{75.14,87.14,99.14}}}|0.0000000000
-11|{{{2.11,13.11,24.11},{35.11,46.11,57.11},{68.11,79.11,90.11}}}|0.0000000000
+10|{{{3.14,15.14,27.14},{39.14,51.14,63.14},{75.14,87.14,99.14}}}|0.000000
+11|{{{2.11,13.11,24.11},{35.11,46.11,57.11},{68.11,79.11,90.11}}}|0.000000
+1|{{{1,2,3},{4,5,6},{7,8,9}}}|8.000000
+2|{{{15,14,13},{12,11,10},{9,8,7}}}|8.000000
+3|{{{1,3,5},{7,9,11},{13,15,17}}}|16.000000
+4|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}|
+5|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|8.000000
+6|{{{NULL,2,NULL},{4,NULL,6},{NULL,8,NULL}}}|
+7|{{{1,3.1,5.2},{NULL,NULL,11.5},{13.6,NULL,NULL}}}|
+8|{{{NULL,3.14,NULL},{NULL,12.56,NULL},{NULL,NULL,25.12}}}|21.980000
+9|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,9}}}|
+10|{{{3.14,15.14,27.14},{39.14,51.14,63.14},{75.14,87.14,99.14}}}|96.000000
+11|{{{2.11,13.11,24.11},{35.11,46.11,57.11},{68.11,79.11,90.11}}}|88.000000