From 1ad6b96d8d2f9c864a8838509d39381939e33bb4 Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Tue, 7 May 2013 15:39:40 +0000 Subject: [PATCH] Added ST_TPI(raster). Ticket #2163. Thanks to Nathaniel Clay for writing the function and required docs and regression tests. git-svn-id: http://svn.osgeo.org/postgis/trunk@11378 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + doc/introduction.xml | 1 + doc/reference_raster.xml | 48 ++++++ raster/rt_pg/rtpostgis.sql.in | 144 ++++++++++++++++++ .../test/regress/rt_elevation_functions.sql | 74 +++++++++ .../regress/rt_elevation_functions_expected | 11 ++ 6 files changed, 279 insertions(+) diff --git a/NEWS b/NEWS index f80f9d187..47bddd366 100644 --- a/NEWS +++ b/NEWS @@ -83,6 +83,7 @@ PostGIS 2.1.0 - #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 * diff --git a/doc/introduction.xml b/doc/introduction.xml index bf8c1411b..813f7c501 100644 --- a/doc/introduction.xml +++ b/doc/introduction.xml @@ -205,6 +205,7 @@ Maxime Guillaud, Maxime van Noppen, Michael Fuhr, Nathan Wagner, +Nathaniel Clay, Nikita Shulga, Norman Vine, Rafal Magda, diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index 38bbb2240..659a307c3 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -8605,6 +8605,54 @@ GROUP BY t1.rast; + + + ST_TPI + Returns a raster with the calculated Topographic Position Index. + + + + + + raster ST_TPI + raster rast + integer nband + raster customextent + text pixeltype="32BF" + boolean interpolate_nodata=FALSE + + + + + + Description + + Calculates the Topographic Position Index, which is defined as the folcal mean with radius of one minus the center cell. + + + This function only supports a focalmean radius of one. + + Availability: 2.1 + + + + Examples + +-- needs examples + + + + + See Also + + , + , + , + + + + + ST_Intersection diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in index 66e605c0b..f6a3d008c 100644 --- a/raster/rt_pg/rtpostgis.sql.in +++ b/raster/rt_pg/rtpostgis.sql.in @@ -3974,6 +3974,150 @@ CREATE OR REPLACE FUNCTION st_hillshade( 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 ----------------------------------------------------------------------- diff --git a/raster/test/regress/rt_elevation_functions.sql b/raster/test/regress/rt_elevation_functions.sql index f73fd08b6..73f0740d8 100644 --- a/raster/test/regress/rt_elevation_functions.sql +++ b/raster/test/regress/rt_elevation_functions.sql @@ -290,3 +290,77 @@ ORDER BY 2, 1, 4, 3; 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); diff --git a/raster/test/regress/rt_elevation_functions_expected b/raster/test/regress/rt_elevation_functions_expected index d2b9ed26c..17221caf0 100644 --- a/raster/test/regress/rt_elevation_functions_expected +++ b/raster/test/regress/rt_elevation_functions_expected @@ -485,3 +485,14 @@ 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 +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 -- 2.50.1