From: Bborie Park Date: Tue, 7 May 2013 15:39:51 +0000 (+0000) Subject: Added ST_Roughness(raster). Ticket #2302. Thanks to Nathaniel Clay for writing the... X-Git-Tag: 2.1.0beta2~48 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=61085d59671c07c6194b9c7a1e56dd32da84f2c3;p=postgis Added ST_Roughness(raster). Ticket #2302. Thanks to Nathaniel Clay for writing the function and required docs and regression tests. git-svn-id: http://svn.osgeo.org/postgis/trunk@11379 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index 47bddd366..c395ba461 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) + - #2302, ST_Roughness for raster (Nathaniel Clay) * Enhancements * diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index 659a307c3..51cdd52b9 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -8329,6 +8329,9 @@ GROUP BY t1.rast; See Also , + , + , + , , @@ -8464,6 +8467,9 @@ GROUP BY t1.rast; See Also , + , + , + , , @@ -8599,6 +8605,9 @@ GROUP BY t1.rast; See Also , + , + , + , , @@ -8617,7 +8626,7 @@ GROUP BY t1.rast; raster ST_TPI raster rast integer nband - raster customextent + raster customextent text pixeltype="32BF" boolean interpolate_nodata=FALSE @@ -8647,6 +8656,53 @@ GROUP BY t1.rast; , , + , + , + , + + + + + + + + ST_Roughness + Returns a raster with the calculated "roughness" of a DEM. + + + + + + raster ST_Roughness + raster rast + integer nband + raster customextent + text pixeltype="32BF" + boolean interpolate_nodata=FALSE + + + + + + Description + Calculates the "roughness" of a DEM, by subtracting the maximum from the minimum for a given area. + Availability: 2.1 + + + + Examples + +-- needs examples + + + + + See Also + + , + , + , + , , diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in index f6a3d008c..8c13476df 100644 --- 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 73f0740d8..d1cee1c9a 100644 --- 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 17221caf0..e5f80b328 100644 --- 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