- #613, ST_SetGeoReference with numerical parameters instead of text
- #2276, ST_AddBand(raster) variant for out-db bands
- #2280, ST_Summary(raster)
+ - #2163, ST_TPI for raster (Nathaniel Clay)
* Enhancements *
Maxime van Noppen,
Michael Fuhr,
Nathan Wagner,
+Nathaniel Clay,
Nikita Shulga,
Norman Vine,
Rafal Magda,
</refsection>
</refentry>
+ <refentry id="RT_ST_TPI">
+ <refnamediv>
+ <refname>ST_TPI</refname>
+ <refpurpose>Returns a raster with the calculated Topographic Position Index.</refpurpose>
+ </refnamediv>
+
+ <refsynopsisdiv>
+ <funcsynopsis>
+ <funcprototype>
+ <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 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 Topographic Position Index, which is defined as the folcal mean with radius of one minus the center cell.</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_TRI" />,
+ <xref linkend="RT_ST_HillShade" />,
+ <xref linkend="RT_ST_Aspect" />
+ </para>
+ </refsection>
+ </refentry>
+
<refentry id="RT_ST_Intersection">
<refnamediv>
<refname>ST_Intersection</refname>
AS $$ SELECT st_hillshade($1, $2, NULL::raster, $3, $4, $5, $6, $7, $8) $$
LANGUAGE 'sql' IMMUTABLE;
+-----------------------------------------------------------------------
+-- ST_TPI
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION _st_tpi4ma(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;
+
+ tpi double precision;
+ mean 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)|--
+ -------------------------------------------------
+
+ 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];
+
+ mean := (Z1 + Z2 + Z3 + Z4 + Z6 + Z7 + Z8 + Z9)/8;
+ tpi := Z5-mean;
+
+ return tpi;
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_tpi(
+ 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_tpi4ma(double precision[][][], integer[][], text[])'::regprocedure,
+ _pixtype,
+ _extenttype, _customextent,
+ 1, 1);
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE;
+
-----------------------------------------------------------------------
-- Get information about the raster
-----------------------------------------------------------------------
DROP TABLE IF EXISTS raster_elevation_out;
DROP TABLE IF EXISTS raster_elevation;
+
+-- ST_TPI
+
+DROP TABLE IF EXISTS raster_value_arrays;
+CREATE TABLE raster_value_arrays (
+ id integer,
+ val double precision[][]
+);
+CREATE OR REPLACE FUNCTION make_value_array(
+ rows integer DEFAULT 3,
+ columns integer DEFAULT 3,
+ start_val double precision DEFAULT 1,
+ step double precision DEFAULT 1,
+ skip_expr text DEFAULT NULL
+)
+ RETURNS double precision[][][]
+ AS $$
+ DECLARE
+ x int;
+ y int;
+ value double precision;
+ values double precision[][][];
+ result boolean;
+ expr text;
+ BEGIN
+ value := start_val;
+
+ values := array_fill(NULL::double precision, ARRAY[1, columns, rows]);
+
+ FOR y IN 1..columns LOOP
+ FOR x IN 1..rows LOOP
+ IF skip_expr IS NULL OR length(skip_expr) < 1 THEN
+ result := TRUE;
+ ELSE
+ expr := replace(skip_expr, '[v]'::text, value::text);
+ EXECUTE 'SELECT (' || expr || ')::boolean' INTO result;
+ END IF;
+
+ IF result IS TRUE THEN
+ values[1][y][x] := value;
+ END IF;
+
+ value := value + step;
+ END LOOP;
+ END LOOP;
+
+ RETURN values;
+ END;
+ $$ LANGUAGE 'plpgsql';
+
+INSERT INTO raster_value_arrays VALUES
+ (1, make_value_array()),
+ (2, make_value_array(3, 3, 15, -1)),
+ (3, make_value_array(3, 3, 1, 2)),
+
+ (4, make_value_array(3, 3, 1, 1, '0')),
+ (5, make_value_array(3, 3, 1, 1, '[v] % 2')),
+ (6, make_value_array(3, 3, 1, 1, '([v] % 2) = 0')),
+ (7, make_value_array(3, 3, 1, 2.1, '([v] NOT IN (7.3, 9.4, 15.7, 17.8))')),
+ (8, make_value_array(3, 3, 0, 3.14, '([v] IN (3.14, 12.56, 25.12))')),
+ (9, make_value_array(3, 3, 1, 1, '[v] > 8')),
+ (10,make_value_array(3,3,3.14,12)),
+ (11,make_value_array(3,3,2.11,11));
+
+
+SELECT
+ id,
+ val,
+ _st_tpi4ma(val,NULL,NULL)::numeric(16, 10) as tpi
+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);
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
+4|{{{NULL,NULL,NULL},{NULL,NULL,NULL},{NULL,NULL,NULL}}}|
+5|{{{1,NULL,3},{NULL,5,NULL},{7,NULL,9}}}|0.0000000000
+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
+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
projects / postgis / commitdiff