From: David Zwarg Date: Fri, 11 Nov 2011 18:50:18 +0000 (+0000) Subject: Added pl/pgsql prototype of the neighborhood map algebra function. X-Git-Tag: 2.0.0alpha1~721 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=1f559f918e02642b58b3a7eea9d082cae60b19e6;p=postgis Added pl/pgsql prototype of the neighborhood map algebra function. git-svn-id: http://svn.osgeo.org/postgis/trunk@8139 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/scripts/plpgsql/st_mapalgebrafctngb.sql b/raster/scripts/plpgsql/st_mapalgebrafctngb.sql new file mode 100644 index 000000000..46cca2aa9 --- /dev/null +++ b/raster/scripts/plpgsql/st_mapalgebrafctngb.sql @@ -0,0 +1,308 @@ +---------------------------------------------------------------------- +-- $Id:$ +-- +-- Copyright (c) 2011 David Zwarg +-- +---------------------------------------------------------------------- + +-- +-- Helper method to get the smallest value in a raster, based on the pixeltype. +-- +CREATE OR REPLACE FUNCTION ST_MinPossibleVal(pixeltype text) + RETURNS float8 AS + $$ + DECLARE + newval int := 0; + BEGIN + newval := CASE + WHEN pixeltype = '1BB' THEN 0 + WHEN pixeltype = '2BUI' THEN 0 + WHEN pixeltype = '4BUI' THEN 0 + WHEN pixeltype = '8BUI' THEN 0 + WHEN pixeltype = '8BSI' THEN -128 + WHEN pixeltype = '16BUI' THEN 0 + WHEN pixeltype = '16BSI' THEN -32768 + WHEN pixeltype = '32BUI' THEN 0 + WHEN pixeltype = '32BSI' THEN -2147483648 + WHEN pixeltype = '32BF' THEN -2147483648 -- Could not find a function returning the smallest real yet + WHEN pixeltype = '64BF' THEN -2147483648 -- Could not find a function returning the smallest float8 yet + END; + RETURN newval; + END; + $$ + LANGUAGE 'plpgsql'; + +-- +--Test rasters +-- +CREATE OR REPLACE FUNCTION ST_TestRaster(h integer, w integer, val float8) + RETURNS raster AS + $$ + DECLARE + BEGIN + RETURN ST_AddBand(ST_MakeEmptyRaster(h, w, 0, 0, 1, 1, 0, 0, -1), '32BF', val, -1); + END; + $$ + LANGUAGE 'plpgsql'; + +-------------------------------------------------------------------- +-- ST_MapAlgebraFctNgb - (one raster version) Return a raster which values +-- are the result of a PLPGSQL user function involving a +-- neighborhood of values from the input raster band. +-- Arguments +-- rast raster - Raster on which the user function is evaluated. +-- band integer - Band number of the raster to be evaluated. Default to 1. +-- pixeltype text - Pixeltype assigned to the resulting raster. User function +-- results are truncated to this type. Default to the +-- pixeltype of the first raster. +-- ngbwidth integer - The width of the neighborhood, in cells. +-- ngbheight integer - The heigh of the neighborhood, in cells. +-- userfunction regprocedure - PLPGSQL user function to apply to with neighborhod pixels. +-- args variadic text[] - Arguments to pass into the user function. +-------------------------------------------------------------------- +DROP FUNCTION IF EXISTS ST_MapAlgebraFctNgb(rast raster, band integer, pixeltype text, ngbwidth integer, ngbheight integer, userfunction text, nodatamode text, variadic args text[]); +CREATE OR REPLACE FUNCTION ST_MapAlgebraFctNgb(rast raster, band integer, pixeltype text, ngbwidth integer, ngbheight integer, userfunction text, nodatamode text, variadic args text[]) + RETURNS raster AS + $$ + DECLARE + width integer; + height integer; + x integer; + y integer; + r float; + newrast raster; + newval float8; + newpixeltype text; + newnodatavalue float; + newinitialvalue float; + neighborhood float[][]; + nrow float[]; + valuesub boolean; + BEGIN + -- Check if raster is NULL + IF rast IS NULL THEN + RAISE NOTICE 'ST_MapAlgebraFctNgb: Raster is NULL. Returning NULL'; + RETURN NULL; + END IF; + + width := ST_Width(rast); + height := ST_Height(rast); + + -- Create a new empty raster having the same georeference as the provided raster + newrast := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(rast), ST_UpperLeftY(rast), ST_ScaleX(rast), ST_ScaleY(rast), ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast)); + + -- If this new raster is empty (width = 0 OR height = 0) then there is nothing to compute and we return it right now + IF ST_IsEmpty(newrast) THEN + RAISE NOTICE 'ST_MapAlgebraFctNgb: Raster is empty. Returning an empty raster'; + RETURN newrast; + END IF; + + -- Check if rast has the required band. Otherwise return a raster without band + IF ST_HasNoBand(rast, band) THEN + RAISE NOTICE 'ST_MapAlgebraFctNgb: Raster does not have the required band. Returning a raster without a band'; + RETURN newrast; + END IF; + + -- Set the new pixeltype + newpixeltype := pixeltype; + IF newpixeltype IS NULL THEN + newpixeltype := ST_BandPixelType(rast, band); + ELSIF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND + newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN + RAISE EXCEPTION 'ST_MapAlgebraFctNgb: Invalid pixeltype "%". Aborting.', newpixeltype; + END IF; + + -- Check for notada value + newnodatavalue := ST_BandNodataValue(rast, band); + IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleVal(newpixeltype) OR newnodatavalue > (-ST_MinPossibleVal(newpixeltype) - 1) THEN + RAISE NOTICE 'ST_MapAlgebraFctNgb: Source raster does not have a nodata value or is out of range for the new raster pixeltype, nodata value for new raster set to the min value possible'; + newnodatavalue := ST_MinPossibleVal(newpixeltype); + END IF; + + -- We set the initial value of the future band to nodata value. + -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleVal + -- but all the values should be recomputed anyway. + newinitialvalue := newnodatavalue; + + -- Optimization: If the raster is only filled with nodata value return right now a raster filled with the nodatavalueexpr + IF ST_BandIsNoData(rast, band) THEN + RETURN ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue); + END IF; + + --Create the raster receiving all the computed values. Initialize it to the new initial value. + newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue); + + IF nodatamode = 'value' THEN + valuesub := TRUE; + ELSIF nodatamode != 'ignore' AND nodatamode != 'NULL' AND NOT nodatamode ~ E'^[\-+]?[0-9]*(\.[0-9]+(e[\-+]?[0-9]+)?)?$' THEN + RAISE NOTICE 'ST_MapAlgebraFctNgb: ''nodatamode'' parameter must be one of: value, ignore, NULL, or a numerical value.'; + RETURN NULL; + END IF; + + -- Computation loop + FOR x IN (1+ngbwidth)..(width-ngbwidth) LOOP + FOR y IN (1+ngbheight)..(height-ngbheight) LOOP + + -- create a matrix of the pixel values in the neighborhood + neighborhood := ARRAY[[]]::float[]; + FOR u IN (x-ngbwidth)..(x+ngbwidth) LOOP + nrow := ARRAY[]::float[]; + FOR v in (y-ngbheight)..(y+ngbheight) LOOP + nrow := nrow || ARRAY[ ST_Value(rast, band, u, v)::float ]; + END LOOP; + neighborhood := neighborhood || ARRAY[ nrow ]; + END LOOP; + + -- RAISE NOTICE 'Neighborhood: %', neighborhood; + + IF valuesub IS TRUE THEN + nodatamode := ST_Value(rast, band, x, y); + -- special case where NULL needs to be quoted + IF nodatamode IS NULL THEN + nodatamode := 'NULL'; + END IF; + END IF; + + EXECUTE 'SELECT ' ||userfunction || '(' || quote_literal(neighborhood) || ', ' || quote_literal(nodatamode) || ', NULL)' INTO r; + + newrast = ST_SetValue(newrast, band, x, y, r); + END LOOP; + END LOOP; + RETURN newrast; + END; + $$ + LANGUAGE 'plpgsql'; + +-- +-- A simple 'callback' user function that sums up all the values in a neighborhood. +-- +CREATE OR REPLACE FUNCTION ST_Sum(matrix float[][], nodatamode text, variadic args text[]) + RETURNS float AS + $$ + DECLARE + x1 integer; + x2 integer; + y1 integer; + y2 integer; + sum float; + BEGIN + sum := 0; + FOR x in array_lower(matrix, 1)..array_upper(matrix, 1) LOOP + FOR y in array_lower(matrix, 2)..array_upper(matrix, 2) LOOP + IF matrix[x][y] IS NULL THEN + CASE nodatamode + WHEN 'NULL' THEN + RETURN NULL; + WHEN 'ignore' THEN + matrix[x][y] := 0; + ELSE -- if the nodatamode was 'value', it was already substituted in the main ST_MapAlgebraFctNgb call. + -- otherwise, this is a literal value that should replace this nodata value + matrix[x][y] := nodatamode::float; + END CASE; + END IF; + sum := sum + matrix[x][y]; + END LOOP; + END LOOP; + RETURN sum; + END; + $$ + LANGUAGE 'plpgsql'; + +-- Tests +-- Test NULL Raster. Should be true. +SELECT ST_MapAlgebraFctNgb(NULL, 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL) IS NULL FROM ST_TestRaster(0, 0, -1) rast; + +-- Test empty Raster. Should be true. +SELECT ST_IsEmpty(ST_MapAlgebraFctNgb(ST_MakeEmptyRaster(0, 10, 0, 0, 1, 1, 1, 1, -1), 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL)); + +-- Test has no band raster. Should be true +SELECT ST_HasNoBand(ST_MapAlgebraFctNgb(ST_MakeEmptyRaster(10, 10, 0, 0, 1, 1, 1, 1, -1), 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL)); + +-- Test has no nodata value. Should return null and 7. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb( + ST_SetBandNoDataValue(rast, NULL), 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL + ), 2, 2) = 7 + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 2, 2, NULL) AS rast; + +-- Test NULL nodatamode. Should return null and null. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) IS NULL + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 2, 2, NULL) AS rast; + +-- Test ignore nodatamode. Should return null and 8. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', 'ignore', NULL), 2, 2 + ) = 8 + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 2, 2, NULL) AS rast; + +-- Test value nodatamode. Should return null and null. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', 'value', NULL), 2, 2 + ) IS NULL + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 2, 2, NULL) AS rast; + +-- Test value nodatamode. Should return null and 9. +SELECT + ST_Value(rast, 1, 1) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', 'value', NULL), 2, 2 + ) = 9 + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 1, 1, NULL) AS rast; + +-- Test value nodatamode. Should return null and 0. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', '-8', NULL), 2, 2 + ) = 0 + FROM ST_SetValue(ST_TestRaster(3, 3, 1), 2, 2, NULL) AS rast; + +-- Test ST_Sum user function. Should be 1 and 9. +SELECT + ST_Value(rast, 2, 2) = 1, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) = 9 + FROM ST_TestRaster(3, 3, 1) AS rast; + +-- Test ST_Sum user function on a no nodata value raster. Should be null and -1. +SELECT + ST_Value(rast, 2, 2) IS NULL, + ST_Value( + ST_MapAlgebraFctNgb(ST_SetBandNoDataValue(rast, NULL), 1, NULL, 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) = -1 + FROM ST_SetValue(ST_TestRaster(3, 3, 0), 2, 2, NULL) AS rast; + +-- Test pixeltype 1. Should return 2 and 15. +SELECT + ST_Value(rast, 2, 2) = 2, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, '4BUI', 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) = 15 + FROM ST_TestRaster(3, 3, 2) AS rast; + +-- Test pixeltype 1. Should return an error. +SELECT + ST_Value(rast, 2, 2), + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, '4BUId', 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) + FROM ST_TestRaster(3, 3, 2) AS rast; + +-- Test pixeltype 1. Should return 1 and 3. +SELECT + ST_Value(rast, 2, 2) = 1, + ST_Value( + ST_MapAlgebraFctNgb(rast, 1, '2BUI', 1, 1, 'ST_Sum', 'NULL', NULL), 2, 2 + ) = 3 + FROM ST_TestRaster(3, 3, 1) AS rast;