From: Bborie Park Date: Tue, 7 May 2013 15:39:58 +0000 (+0000) Subject: Added ST_TRI(raster). Ticket #2164. Thanks to Nathaniel Clay for writing the function... X-Git-Tag: 2.1.0beta2~47 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=a7e0ca4fe3026ee1b9c24f5999e9f914ecf57136;p=postgis Added ST_TRI(raster). Ticket #2164. Thanks to Nathaniel Clay for writing the function and required docs and regression tests. git-svn-id: http://svn.osgeo.org/postgis/trunk@11380 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index c395ba461..0e55bf37e 100644 --- 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 51cdd52b9..8988bf3b1 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -8709,6 +8709,59 @@ GROUP BY t1.rast; + + + ST_TRI + Returns a raster with the calculated Terrain Ruggedness Index. + + + + + + raster ST_TRI + raster rast + integer nband + raster customextent + text pixeltype="32BF" + boolean interpolate_nodata=FALSE + + + + + + Description + + + Terrain Ruggedness Index is calculated by comparing a central pixel with its neighbors, taking the absolute values of the differences, and averaging the result. + + + + 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 8c13476df..ba55bfd3b 100644 --- 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 d1cee1c9a..7780fea8f 100644 --- 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 e5f80b328..aae515840 100644 --- 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