Changes between Version 2 and Version 3 of PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/code


Ignore:
Timestamp:
Aug 28, 2012, 10:39:32 AM (12 years ago)
Author:
qliu
Comment:

—

Legend:

Unmodified
Added
Removed
Modified

  • PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/code

    v2 v3  
    1919#!sql
    2020--------------------------------------------------------------------
    21 -- ST_EuclideanDistance - (for one raster) Return a raster generated from
    22 --                        a set of raster specifications, which values are
    23 --                                                the Euclidean distance from pixel to the points
    24 --                        (for now only concerning point geometries)
    25 --                        from the input source table of geometries
     21-- ST_NearestGridCentroid()
     22-- Returns a point geometry representing centroid of the pixel
     23-- intersecting with the geometry passed in to the function.
     24-- Accepts point geometry only.
     25--------------------------------------------------------------------
     26DROP FUNCTION IF EXISTS ST_NearestGridCentroid(raster,geometry);
     27CREATE OR REPLACE FUNCTION ST_NearestGridCentroid(rast raster,point geometry)
     28        RETURNS geometry AS
     29        $$
     30        DECLARE
     31                geomtype text;
     32                pixelcentroid geometry;
     33        BEGIN
     34                -- Accept Point geometry only
     35                geomtype := ST_GeometryType(point);
     36                IF geomtype != 'ST_Point' THEN
     37                        RAISE EXCEPTION 'ST_ST_NearestGridCentroid: Invalid geometry type "%". Aborting.',geomtype;
     38                END IF;
     39                -- Find pixel centroid of the point geometry
     40                pixelcentroid := ST_PixelAsCentroid(rast,ST_World2RasterCoordX(rast,point),ST_World2RasterCoordY(rast,point));
     41                RETURN pixelcentroid;
     42        END;
     43        $$
     44        LANGUAGE 'plpgsql' IMMUTABLE;
     45       
     46--------------------------------------------------------------------
     47-- Custom function called by ST_MapAlgebraFct()
     48-- to calculate Euclidean distance for current pixel.
     49--------------------------------------------------------------------
     50DROP FUNCTION IF EXISTS euclidean_distance_fct(float,integer[],text[]);
     51CREATE OR REPLACE FUNCTION euclidean_distance_fct(val float,pos integer[],VARIADIC args text[])
     52        RETURNS float8 AS
     53        $$
     54        DECLARE
     55                newrast raster;
     56                pixelcentroidgeom geometry;
     57                nearestngb geometry;
     58                exesql text;
     59                maxdist float8;
     60                pixelval float8;
     61        BEGIN
     62                -- Make an empty raster with the raster specifications passed in to the custom function
     63                newrast := ST_MakeEmptyRaster(args[1]::integer,args[2]::integer,
     64                                                                          args[3]::float8,args[4]::float8,
     65                                                                          args[5]::float8,args[6]::float8,
     66                                                                          args[7]::float8,args[8]::float8,
     67                                                                          args[9]::integer);
     68                -- Get pixel centroid point geometry of current pixel
     69                pixelcentroidgeom := ST_PixelAsCentroid(newrast,pos[1],pos[2]);
     70                -- Use KNN index to find the nearest source point
     71                exesql := 'SELECT ' ||  quote_ident(args[13]) || ' FROM ' || quote_ident(args[11]) || '.' || quote_ident(args[12]) || ' ORDER BY ' || quote_ident(args[13]) || ' <-> ST_GeomFromText(' || quote_literal(ST_AsText(pixelcentroidgeom)) || ',' || args[9] ||') LIMIT 1;';
     72                EXECUTE exesql INTO nearestngb;
     73                -- If snap is set True, Euclidean Distance is calculated from current pixel centroid to the centroid of the pixel intersecting with its nearest source point geometry
     74                -- If snap is set False, Euclidean Distance is calculated as actural distance from current pixel centroid to its nearest source point geometry
     75                IF args[14]::boolean THEN
     76                        pixelval := ST_Distance(pixelcentroidgeom,ST_NearestGridCentroid(newrast,nearestngb));
     77                ELSE
     78                        pixelval := ST_Distance(pixelcentroidgeom,nearestngb);
     79                END IF;
     80                -- If max distance is specified
     81                maxdist := args[15];
     82                IF maxdist != -1 THEN
     83                        -- Any distance exceeds max distance is assigned nodata value
     84                        IF pixelval > maxdist THEN
     85                                pixelval := args[10]::float8;
     86                        END IF;
     87                END IF;
     88                RETURN pixelval;
     89        END;
     90        $$
     91        LANGUAGE 'plpgsql' IMMUTABLE;
     92       
     93--------------------------------------------------------------------
     94-- ST_EuclideanDistance - Return a raster generated from a set of raster specifications.
     95--                                                Pixel values are the Euclidean distance from pixel to source points
     96--                        which is a table of point geometries
    2697--                        (points for now, lines or polygons in the future).
    2798-- Arguments
    2899-- width,height,rastulx,rastuly,rastscalex,rastscaley,rastskewx,rastskewy
    29 --  - Dimenssions and georeference of resulted raster.
     100--                      - Dimenssions and georeference of resulted raster.
    30101-- rastsrid - SRID of the resulted raster
    31102-- rastndv - Nodata value for resulted raster
    32 -- pixeltype text - Pixeltype assigned to the resulting raster. By default
    33 --                  to be the pixeltype of the reference raster if not specified.
    34 -- sourceschema text - Database schema of the source geometry table.
    35 -- sourcetable text - Source geometry table name.
    36 -- sourcegeomcolumn text - Source geometry table geometry column name.
    37 -- maxdist double precision - Maximum distance from pixel to the source
    38 --                            when the source is too far from the pixel.
    39 --                            Pixel that exceeds it will be assigned nodatavalue.
     103-- pixeltype - Pixeltype assigned to the resulting raster.
     104--                         As distance raster, we only accept pixel types of:
     105--                         '4BUI','8BSI','8BUI','16BSI','16BUI','32BSI','32BUI','32BF','64BF'.
     106-- sourceschema - Database schema of the source geometry table.
     107-- sourcetable - Database table name of the source geometry table.
     108-- sourcegeomcolumn - Database geometry column name of the source geometry table.
     109-- snaptocentroid - When True, the option will "snap" any source point geometry to the centroid of
     110--                                      the pixel intersecting with it, resulting in a distance value equal to zero.
     111--                                      The Euclidean distance is calculated from the centroid of each pixel to the
     112--                                      centroid of the pixel intersecting with its nearest source point geometry.
     113-- maxdist - Maximum distance from pixel to the source when the source is too far away from the pixel.
     114--                       Pixel that exceeds maximum distance will be assigned nodatavalue.
     115--                       By default to be -1 if not specified.
    40116--------------------------------------------------------------------
    41 DROP FUNCTION IF EXISTS ST_EuclideanDistance(width integer,
    42                                                                                          height integer,
    43                                                                                          rastulx float8,
    44                                                                                          rastuly float8,
    45                                                                                          rastscalex float8,
    46                                                                                          rastscaley float8,
    47                                                                                          rastskewx float8,
    48                                                                                          rastskewy float8,
    49                                                                                          rastsrid integer,
    50                                                                                          rastndv float8,
    51                                                                                          pixeltype text,
    52                                                                                          sourceschema text,
    53                                                                                          sourcetable text,
    54                                                                                          sourcegeomcolumn text,
    55                                                                                          double precision = -1);
     117DROP FUNCTION IF EXISTS ST_EuclideanDistance(integer,integer,float8,float8,float8,float8,float8,float8,integer,float8,text,text,text,text,boolean,float8);
    56118CREATE OR REPLACE FUNCTION ST_EuclideanDistance(width integer,
    57119                                                                                                height integer,
    … …  
    68130                                                                                                sourcetable text,
    69131                                                                                                sourcegeomcolumn text,
    70                                                                                                 double precision = -1)
     132                                                                                                snaptocentroid boolean,
     133                                                                                                float8 = -1)
    71134    RETURNS raster AS
    72135    $$
    73136    DECLARE
    74                 maxdist ALIAS FOR $15;
    75         x integer;
    76         y integer;
    77                 sourcex float8;
    78                 sourcey float8;
    79                 sourcexr integer;
    80                 sourceyr integer;
     137                maxdist ALIAS FOR $16;
    81138                sourcesrid integer;
    82                 newx float8;
    83                 newy float8;
    84                 exesql text;
    85139        newnodatavalue float8;
    86140                newinitialvalue float8;
    87141        newrast raster;
    88         newval float8;
    89142        newpixeltype text;
    90                 -- Assuming reference raster has only one band
    91                 band integer := 1;
    92                 geom geometry;
    93143    BEGIN
    94144                -- Assuming source point table has one SRID for all geometries
    95                 EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid;
     145                EXECUTE 'SELECT DISTINCT ST_SRID(' || quote_ident(sourcegeomcolumn) || ') FROM ' || quote_ident(sourceschema) || '.' || quote_ident(sourcetable) || ';' INTO sourcesrid;
    96146
    97147        -- Create a new empty raster using provided georeference
    … …  
    103153            RETURN newrast;
    104154        END IF;
    105        
     155
    106156        -- Set the new pixeltype
    107157        newpixeltype := pixeltype;
    108                 IF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND
    109            newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN
     158                IF newpixeltype = '1BB' OR newpixeltype = '2BUI' OR (newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND
     159           newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF') THEN
    110160            RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
    111161        END IF;
    … …  
    120170        -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue
    121171        newinitialvalue := newnodatavalue;
    122        
     172
    123173        -- Create the raster receiving all the computed distance values. Initialize it to the new initial value.
    124174        newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
    125175
    126                 -- Set raster pixel value to zero if pixel is source point
    127                 exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';';
    128                 FOR geom IN EXECUTE(exesql) LOOP
    129                         sourcex := ST_X(geom);
    130                         sourcey := ST_Y(geom);
    131                         sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey);
    132                         sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey);
    133                         -- If source point within raster range, then set pixel value to zero
    134                         IF sourcexr <= width AND sourceyr <= height THEN
    135                                 newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
    136                         END IF;
    137                 END LOOP;
    138                
    139                 -- Scan pixels in the raster set pixel value to the computed Eculidean distance
    140                 -- to its nearest source point.
    141                 -- There is place for optimization here for a more efficient scanning method
    142         FOR x IN 1..width LOOP
    143             FOR y IN 1..height LOOP
    144                                 newx := ST_Raster2WorldCoordX(newrast, x, y);
    145                                 newy := ST_Raster2WorldCoordY(newrast, x, y);
    146                                 -- If pixel is not the source point then compute Eculidean distance
    147                                 IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN
    148                                         exesql := 'SELECT ST_Distance(ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || '), (SELECT ' ||  sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ' ORDER BY ' || sourcegeomcolumn || ' <-> ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || ') LIMIT 1));';
    149                                         EXECUTE exesql INTO newval;
    150                                         -- If max distance specified
    151                                         IF maxdist != -1 THEN
    152                                                 -- Any distance exceeds max distance is assigned nodata value
    153                                                 IF newval > maxdist THEN
    154                                                         newval = newnodatavalue;
    155                                                 END IF;
    156                                         END IF;
    157                                 END IF;
    158                                 newrast := ST_SetValue(newrast, band, x, y, newval);
    159             END LOOP;
    160         END LOOP;
    161 
     176                -- Call ST_MapAlgebraFct() with custom EuclideanDistance function to set the new raster to the resulted distance raster
     177                newrast := ST_MapAlgebraFct(newrast,NULL,'euclidean_distance_fct(float,integer[],text[])'::regprocedure,
     178                                                                        width::text,height::text,
     179                                                                        rastulx::text, rastuly::text,
     180                                                                        rastscalex::text,rastscaley::text,
     181                                                                        rastskewx::text, rastskewy::text,
     182                                                                        rastsrid::text, newnodatavalue::text,
     183                                                                        sourceschema,sourcetable,sourcegeomcolumn,
     184                                                                        snaptocentroid::text,maxdist::text);
     185               
    162186        RETURN newrast;
    163187    END;
    164188    $$
    165189    LANGUAGE 'plpgsql';
    166        
    167 -- test
    168 -- CREATE TABLE test_eudist_2 AS (SELECT 1 AS rid,ST_EuclideanDistance(10, 10, 0, 0, 9, 9, 0, 0, 4326,-1,'32BF','public','test_geom','the_geom') AS rast FROM test_raster);
    169 -- CREATE TABLE test_eudist_22 AS (SELECT 1 AS rid,ST_EuclideanDistance(10, 10, 0, 0, 9, 9, 0, 0, 4326,-1,'32BF','public','test_geom','the_geom',4954880) AS rast FROM test_raster);
    170 
    171190
    172191--------------------------------------------------------------------
    173 -- ST_EuclideanDistance - (for one raster) Return a raster generated
    174 --                                                from a one band reference raster,
    175 --                        which values are the Euclidean distance
    176 --                        from pixel to the points
    177 --                                                (for now only concerning point geometries)
    178 --                        from the input source table of geometries
     192-- ST_EuclideanDistance - Return a raster generated from a one-band reference raster.
     193--                                                Pixel values are the Euclidean distance from pixel to source points
     194--                        which is a table of point geometries
    179195--                        (points for now, lines or polygons in the future).
    180196-- Arguments
    181197-- refrast raster -  Reference raster align with which Euclidean distance
    182198--                                       raster will be created.
    183 -- pixeltype text - Pixeltype assigned to the resulting raster. By default
    184 --                  to be the pixeltype of the reference raster if not specified.
    185 -- sourceschema text - Database schema of the source geometry table.
    186 -- sourcetable text - Source geometry table name.
    187 -- sourcegeomcolumn text - Source geometry table geometry column name.
    188 -- maxdist double precision - Maximum distance from pixel to the source
    189 --                            when the source is too far from the pixel.
    190 --                            Pixel that exceeds it will be assigned nodatavalue.
     199-- rastsrid - SRID of the resulted raster
     200-- pixeltype - Pixeltype assigned to the resulting raster.
     201--                         As distance raster, we only accept pixel types of:
     202--                        '4BUI','8BSI','8BUI','16BSI','16BUI','32BSI','32BUI','32BF','64BF'.
     203-- sourceschema - Database schema of the source geometry table.
     204-- sourcetable - Database table name of the source geometry table.
     205-- sourcegeomcolumn - Database geometry column name of the source geometry table.
     206-- snaptocentroid - When True, the option will "snap" any source point geometry to the centroid of
     207--                                      the pixel intersecting with it, resulting in a distance value equal to zero.
     208--                                      The Euclidean distance is calculated from the centroid of each pixel to the
     209--                                      centroid of the pixel intersecting with its nearest source point geometry.
     210-- maxdist - Maximum distance from pixel to the source when the source is too far away from the pixel.
     211--                       Pixel that exceeds maximum distance will be assigned nodatavalue.
     212--                       By default to be -1 if not specified.
    191213--------------------------------------------------------------------
    192 DROP FUNCTION IF EXISTS ST_EuclideanDistance(refrast raster,
    193                                                                                          pixeltype text,
    194                                                                                          sourceschema text,
    195                                                                                          sourcetable text,
    196                                                                                          sourcegeomcolumn text,
    197                                                                                          double precision = -1);
     214DROP FUNCTION IF EXISTS ST_EuclideanDistance(raster,text,text,text,text,boolean,float8);
    198215CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster,
    199216                                                                                                pixeltype text,
    … …  
    201218                                                                                                sourcetable text,
    202219                                                                                                sourcegeomcolumn text,
    203                                                                                                 double precision = -1)
    204     RETURNS raster AS
    205     $$
    206     DECLARE
    207                 maxdist ALIAS FOR $6;
    208         width integer;
    209         height integer;
    210         x integer;
    211         y integer;
    212                 sourcex float8;
    213                 sourcey float8;
    214                 sourcexr integer;
    215                 sourceyr integer;
    216                 sourcesrid integer;
    217                 newsrid integer;
    218                 newx float8;
    219                 newy float8;
    220                 exesql text;
    221         newnodatavalue float8;
    222                 newinitialvalue float8;
    223         newrast raster;
    224         newval float8;
    225         newpixeltype text;
    226                 -- Assuming reference raster has only one band
    227                 band integer := 1;
    228                 geom geometry;
    229     BEGIN
    230         -- Check if reference raster is NULL
    231         IF refrast IS NULL THEN
    232             RAISE NOTICE 'ST_EuclideanDistance: Reference raster is NULL. Returning NULL';
    233             RETURN NULL;
    234         END IF;
    235                
    236         width := ST_Width(refrast);
    237         height := ST_Height(refrast);
    238                 newsrid := ST_SRID(refrast);
    239                 -- Assuming source point table has one SRID for all geometries
    240                 EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid;
    241                
    242         -- Create a new empty raster having the same georeference as the reference raster
    243         newrast := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(refrast), ST_UpperLeftY(refrast), ST_ScaleX(refrast), ST_ScaleY(refrast), ST_SkewX(refrast), ST_SkewY(refrast), newsrid);
    244                
    245         -- Check if the new raster is empty (width = 0 OR height = 0)
    246         IF ST_IsEmpty(newrast) THEN
    247             RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster';
    248             RETURN newrast;
    249         END IF;
    250        
    251         -- Set the new pixeltype
    252         newpixeltype := pixeltype;
    253                 -- If pixeltype not specified then use the pixeltype of the reference raster
    254         IF newpixeltype IS NULL THEN
    255             newpixeltype := ST_BandPixelType(refrast, band);
    256         ELSIF newpixeltype != '1BB' AND newpixeltype != '2BUI' AND newpixeltype != '4BUI' AND newpixeltype != '8BSI' AND newpixeltype != '8BUI' AND
    257                newpixeltype != '16BSI' AND newpixeltype != '16BUI' AND newpixeltype != '32BSI' AND newpixeltype != '32BUI' AND newpixeltype != '32BF' AND newpixeltype != '64BF' THEN
    258             RAISE EXCEPTION 'ST_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
    259         END IF;
    260                
    261         -- Check for notada value
    262         newnodatavalue := ST_BandNodataValue(refrast, band);
    263         IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN
    264             RAISE NOTICE 'ST_EuclideanDistance: Reference raster do 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';
    265             newnodatavalue := ST_MinPossibleValue(newpixeltype);
    266         END IF;
    267         -- We set the initial value of the resulting raster to nodata value.
    268         -- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue
    269         newinitialvalue := newnodatavalue;
    270        
    271         -- Create the raster receiving all the computed distance values. Initialize it to the new initial value.
    272         newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
    273                
    274                 -- Set raster pixel value to zero if pixel is source point
    275                 exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';';
    276                 FOR geom IN EXECUTE(exesql) LOOP
    277                         sourcex := ST_X(geom);
    278                         sourcey := ST_Y(geom);
    279                         sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey);
    280                         sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey);
    281                         -- If source point within raster range, then set pixel value to zero
    282                         IF sourcexr <= width AND sourceyr <= height THEN
    283                                 newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
    284                         END IF;
    285                 END LOOP;
    286                
    287                 -- Scan pixels in the raster set pixel value to the computed Eculidean distance
    288                 -- to its nearest source point.
    289                 -- There is place for optimization here for a more efficient scanning method
    290         FOR x IN 1..width LOOP
    291             FOR y IN 1..height LOOP
    292                                 newx := ST_Raster2WorldCoordX(newrast, x, y);
    293                                 newy := ST_Raster2WorldCoordY(newrast, x, y);
    294                                 -- If pixel is not the source point then compute Eculidean distance
    295                                 IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN
    296                                         exesql := 'SELECT ST_Distance(ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || '), (SELECT ' ||  sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ' ORDER BY ' || sourcegeomcolumn || ' <-> ST_SetSRID(ST_MakePoint(' || newx || ',' || newy || '),' || sourcesrid || ') LIMIT 1));';
    297                                         EXECUTE exesql INTO newval;
    298                                         -- If max distance specified
    299                                         IF maxdist != -1 THEN
    300                                                 -- Any distance exceeds max distance is assigned nodata value
    301                                                 IF newval > maxdist THEN
    302                                                         newval := newnodatavalue;
    303                                                 END IF;
    304                                         END IF;
    305                                 END IF;
    306                                 newrast := ST_SetValue(newrast, band, x, y, newval);
    307             END LOOP;
    308         END LOOP;
    309 
    310         RETURN newrast;
    311     END;
    312     $$
    313     LANGUAGE 'plpgsql';
    314 
    315 -- test
    316 -- CREATE TABLE test_eudist AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom') AS rast FROM test_raster);
    317 -- CREATE TABLE test_eudist2 AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom',4954880) AS rast FROM test_raster);
     220                                                                                                snaptocentroid boolean,
     221                                                                                                float8 = -1)
     222    RETURNS raster
     223    AS $$ SELECT ST_EuclideanDistance(ST_Width($1),ST_Height($1),ST_UpperLeftX($1),ST_UpperLeftY($1),ST_ScaleX($1),ST_ScaleY($1),ST_SkewX($1),ST_SkewY($1),ST_SRID($1),ST_BandNodataValue($1,1),$2,$3,$4,$5,$6,$7) $$
     224    LANGUAGE 'sql';
    318225}}}
    319226