Added ST_TRI(raster). Ticket #2164. Thanks to Nathaniel Clay for writing the function...

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

diff --git a/NEWS b/NEWS
index c395ba4612ff20602616fd8ed69a698ef5a66a0c..0e55bf37ee4bfaaa32184f274b80c1787131d6f1 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)
+  - #2164, ST_TRI 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 51cdd52b9debcaa7da240c9bfb061046cde29d49..8988bf3b1a318e17f8cfc5b8c6d908297071788f 100644 (file)
--- a/doc/reference_raster.xml
+++ b/doc/reference_raster.xml
@@ -8709,6 +8709,59 @@ GROUP BY t1.rast;
                        </refsection>
                </refentry>
 
+               <refentry id="RT_ST_TRI">
+                       <refnamediv>
+                               <refname>ST_TRI</refname>
+                               <refpurpose>Returns a raster with the calculated Terrain Ruggedness Index.</refpurpose>
+                       </refnamediv>
+
+                       <refsynopsisdiv>
+                               <funcsynopsis>
+                                       <funcprototype>
+                                               <funcdef>raster <function>ST_TRI</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>
+                                       Terrain Ruggedness Index is calculated by comparing a central pixel with its neighbors, taking the absolute values of the differences, and averaging the result.
+                               </para>
+
+                               <note>
+                                       <para>This function only supports a focalmean radius of one.</para>
+                               </note>
+
+                               <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_Roughness" />, 
+                                               <xref linkend="RT_ST_TPI" />, 
+                                               <xref linkend="RT_ST_Slope" />, 
+                                               <xref linkend="RT_ST_HillShade" />, 
+                                               <xref linkend="RT_ST_Aspect" />
+                                       </para>
+                       </refsection>
+               </refentry>
+
                <refentry id="RT_ST_Intersection">
                        <refnamediv>
                                <refname>ST_Intersection</refname>

diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in
index 8c13476df08b644a21da9f9058787eb4247733a0..ba55bfd3b9c9ead20104d8981cc68d3cac87152b 100644 (file)
--- a/raster/rt_pg/rtpostgis.sql.in
+++ b/raster/rt_pg/rtpostgis.sql.in
@@ -4241,6 +4241,149 @@ CREATE OR REPLACE FUNCTION st_roughness(
        END;
        $$ LANGUAGE 'plpgsql' IMMUTABLE;
 
+-----------------------------------------------------------------------
+-- ST_TRI
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION _st_tri4ma(value double precision[][][], pos integer[][], VARIADIC userargs text[] DEFAULT NULL)
+       RETURNS double precision
+       AS $$
+       DECLARE
+               x integer;
+               y integer;
+               z integer;
+
+               Z1 double precision;
+               Z2 double precision;
+               Z3 double precision;
+               Z4 double precision;
+               Z5 double precision;
+               Z6 double precision;
+               Z7 double precision;
+               Z8 double precision;
+               Z9 double precision;
+
+               tri 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;
+
+               -------------------------------------------------
+               --|   Z1= Z(-1,1) |  Z2= Z(0,1) | Z3= Z(1,1)  |--
+               -------------------------------------------------
+               --|   Z4= Z(-1,0) |  Z5= Z(0,0) | Z6= Z(1,0)  |--
+               -------------------------------------------------
+               --|   Z7= Z(-1,-1)|  Z8= Z(0,-1)|  Z9= Z(1,-1)|--
+               -------------------------------------------------
+
+               -- _scale width and height units / z units to make z units equal to height width units
+               Z1 := _value[z][1][1];
+               Z2 := _value[z][2][1];
+               Z3 := _value[z][3][1];
+               Z4 := _value[z][1][2];
+               Z5 := _value[z][2][2];
+               Z6 := _value[z][3][2];
+               Z7 := _value[z][1][3];
+               Z8 := _value[z][2][3];
+               Z9 := _value[z][3][3];
+
+               tri := ( abs(Z1 - Z5 ) + abs( Z2 - Z5 ) + abs( Z3 - Z5 ) + abs( Z4 - Z5 ) + abs( Z6 - Z5 ) + abs( Z7 - Z5 ) + abs( Z8 - Z5 ) + abs ( Z9 - Z5 )) / 8;
+               
+               return tri;  
+       END;
+       $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_tri(
+       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;
+
+               -- 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_tri4ma(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 d1cee1c9ab785b53cb3185b7f3ccbbba4ff89bb3..7780fea8f9c1fb0e2b4deb354887d4e90d0f1f22 100644 (file)
--- a/raster/test/regress/rt_elevation_functions.sql
+++ b/raster/test/regress/rt_elevation_functions.sql
@@ -369,5 +369,12 @@ SELECT
 FROM raster_value_arrays
 ORDER BY id;
 
+SELECT
+       id,
+       val,
+       round(_st_tri4ma(val,NULL,NULL)::numeric, 6) as tri
+FROM raster_value_arrays
+ORDER BY id;
+
 DROP TABLE IF EXISTS raster_value_arrays;
 DROP FUNCTION IF EXISTS make_value_array(integer,integer,double precision, double precision, text);

diff --git a/raster/test/regress/rt_elevation_functions_expected b/raster/test/regress/rt_elevation_functions_expected
index e5f80b3282f18e6954ae4d92cc605b78b58971d7..aae5158406d047e27e352493e132bc3b37f1081c 100644 (file)
--- a/raster/test/regress/rt_elevation_functions_expected
+++ b/raster/test/regress/rt_elevation_functions_expected
@@ -507,3 +507,14 @@ DO
 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
+1|{{{1,2,3},{4,5,6},{7,8,9}}}|2.500000
+2|{{{15,14,13},{12,11,10},{9,8,7}}}|2.500000
+3|{{{1,3,5},{7,9,11},{13,15,17}}}|5.000000
+4|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}|
+5|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|1.500000
+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}}}|2.747500
+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}}}|30.000000
+11|{{{2.11,13.11,24.11},{35.11,46.11,57.11},{68.11,79.11,90.11}}}|27.500000