File: st_reclass.sql

package info (click to toggle)
postgis 2.1.4%2Bdfsg-3
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 35,424 kB
  • ctags: 8,886
  • sloc: sql: 113,491; ansic: 97,254; xml: 41,127; sh: 11,925; java: 5,662; perl: 3,113; makefile: 2,265; python: 1,198; yacc: 438; lex: 114
file content (70 lines) | stat: -rw-r--r-- 2,885 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
----------------------------------------------------------------------
--
-- $Id: st_reclass.sql 9324 2012-02-27 22:08:12Z pramsey $
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_Reclass() function is already implemented in C. This plpgsql script is provided only as an example. 
-- Defining the plpgsql function below might overwrite the current C implementation and brake other functions dependent on it.
-- Use with caution.
----------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_Reclass(rast raster,
                      band int,
                                      reclassexpr text) 
    RETURNS raster AS 
    $$
    DECLARE
    -- Create a new raster without the band we will reclassify
        newrast raster := ST_DeleteBand(rast, band);

    -- Determine the nodata value
        nodataval float8 := ST_BandNodataValue(rast, band);

    -- First parse of the reclass expression. Split the reclass classes into an array.
        reclarray text[] := string_to_array(reclassexpr, '|');

        -- Determine the number of classes.
        nbreclassstr int := array_length(reclarray, 1);

        -- Initialise the MapAlgebra expression
        maexpr text := 'CASE ';

        -- Temporary container for the two part of the class being parsed.
        reclassstr text[];

        -- Temporary array containing the splited class.
        fromstr text[];

        i int;
    BEGIN
    -- For each classes
        FOR i IN 1..nbreclassstr LOOP
        -- Split the class into an array of classes.
            reclassstr := string_to_array(reclarray[i], ':');
            IF array_length(reclassstr, 1) < 2 THEN
            RAISE EXCEPTION 'ST_Reclass: Invalid reclassification class: "%". Aborting', reclarray[i];
            END IF;
            -- Split the range to reclassify into two
            fromstr := string_to_array(reclassstr[1], '-');
            -- Replace nodata with the nodata value
            IF upper(reclassstr[2]) = 'NODATA' THEN
                reclassstr[2] = nodataval::text;
            END IF;
            -- Build the maexpr expression
            IF fromstr[2] IS NULL OR fromstr[2] = ''  THEN
            maexpr := maexpr || ' WHEN ' || fromstr[1] || ' = rast THEN ' || reclassstr[2] || ' ';
        ELSE
            maexpr := maexpr || ' WHEN ' || fromstr[1] || ' <= rast AND rast < ' || fromstr[2] || ' THEN ' || reclassstr[2] || ' ';
        END IF;
        END LOOP;
        maexpr := maexpr || 'ELSE rast END';
        newrast := ST_AddBand(rast, ST_MapAlgebra(rast, band, maexpr), 1, band);
    RETURN newrast;
    END;
    $$
    LANGUAGE 'plpgsql';

SELECT ST_Value(ST_TestRaster(1, 1, 4),1,1)
SELECT ST_Value(ST_Reclass(ST_TestRaster(1, 1, 4), 1, '1:2|2:2|3-5:10'), 1, 1);