PL/pgSQL script
Euclidean Distance
ST_EuclideanDistance()
--------------------------------------------------------------------
-- ST_NearestGridCentroid()
-- Returns a point geometry representing centroid of the pixel
-- intersecting with the geometry passed in to the function.
-- Accepts point geometry only.
--------------------------------------------------------------------
DROP FUNCTION IF EXISTS ST_NearestGridCentroid(raster,geometry);
CREATE OR REPLACE FUNCTION ST_NearestGridCentroid(rast raster,point geometry)
RETURNS geometry AS
$$
DECLARE
geomtype text;
pixelcentroid geometry;
BEGIN
-- Accept Point geometry only
geomtype := ST_GeometryType(point);
IF geomtype != 'ST_Point' THEN
RAISE EXCEPTION 'ST_ST_NearestGridCentroid: Invalid geometry type "%". Aborting.',geomtype;
END IF;
-- Find pixel centroid of the point geometry
pixelcentroid := ST_PixelAsCentroid(rast,ST_World2RasterCoordX(rast,point),ST_World2RasterCoordY(rast,point));
RETURN pixelcentroid;
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;
--------------------------------------------------------------------
-- Custom function called by ST_MapAlgebraFct()
-- to calculate Euclidean distance for current pixel.
--------------------------------------------------------------------
DROP FUNCTION IF EXISTS euclidean_distance_fct(float,integer[],text[]);
CREATE OR REPLACE FUNCTION euclidean_distance_fct(val float,pos integer[],VARIADIC args text[])
RETURNS float8 AS
$$
DECLARE
newrast raster;
pixelcentroidgeom geometry;
nearestngb geometry;
exesql text;
maxdist float8;
pixelval float8;
BEGIN
-- Make an empty raster with the raster specifications passed in to the custom function
newrast := ST_MakeEmptyRaster(args[1]::integer,args[2]::integer,
args[3]::float8,args[4]::float8,
args[5]::float8,args[6]::float8,
args[7]::float8,args[8]::float8,
args[9]::integer);
-- Get pixel centroid point geometry of current pixel
pixelcentroidgeom := ST_PixelAsCentroid(newrast,pos[1],pos[2]);
-- Use KNN index to find the nearest source point
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;';
EXECUTE exesql INTO nearestngb;
-- 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
-- If snap is set False, Euclidean Distance is calculated as actural distance from current pixel centroid to its nearest source point geometry
IF args[14]::boolean THEN
pixelval := ST_Distance(pixelcentroidgeom,ST_NearestGridCentroid(newrast,nearestngb));
ELSE
pixelval := ST_Distance(pixelcentroidgeom,nearestngb);
END IF;
-- If max distance is specified
maxdist := args[15];
IF maxdist != -1 THEN
-- Any distance exceeds max distance is assigned nodata value
IF pixelval > maxdist THEN
pixelval := args[10]::float8;
END IF;
END IF;
RETURN pixelval;
END;
$$
LANGUAGE 'plpgsql' IMMUTABLE;
--------------------------------------------------------------------
-- ST_EuclideanDistance - Return a raster generated from a set of raster specifications.
-- Pixel values are the Euclidean distance from pixel to source points
-- which is a table of point geometries
-- (points for now, lines or polygons in the future).
-- Arguments
-- width,height,rastulx,rastuly,rastscalex,rastscaley,rastskewx,rastskewy
-- - Dimenssions and georeference of resulted raster.
-- rastsrid - SRID of the resulted raster
-- rastndv - Nodata value for resulted raster
-- pixeltype - Pixeltype assigned to the resulting raster.
-- As distance raster, we only accept pixel types of:
-- '4BUI','8BSI','8BUI','16BSI','16BUI','32BSI','32BUI','32BF','64BF'.
-- sourceschema - Database schema of the source geometry table.
-- sourcetable - Database table name of the source geometry table.
-- sourcegeomcolumn - Database geometry column name of the source geometry table.
-- snaptocentroid - When True, the option will "snap" any source point geometry to the centroid of
-- the pixel intersecting with it, resulting in a distance value equal to zero.
-- The Euclidean distance is calculated from the centroid of each pixel to the
-- centroid of the pixel intersecting with its nearest source point geometry.
-- maxdist - Maximum distance from pixel to the source when the source is too far away from the pixel.
-- Pixel that exceeds maximum distance will be assigned nodatavalue.
-- By default to be -1 if not specified.
--------------------------------------------------------------------
DROP FUNCTION IF EXISTS ST_EuclideanDistance(integer,integer,float8,float8,float8,float8,float8,float8,integer,float8,text,text,text,text,boolean,float8);
CREATE OR REPLACE FUNCTION ST_EuclideanDistance(width integer,
height integer,
rastulx float8,
rastuly float8,
rastscalex float8,
rastscaley float8,
rastskewx float8,
rastskewy float8,
rastsrid integer,
rastndv float8,
pixeltype text,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
snaptocentroid boolean,
float8 = -1)
RETURNS raster AS
$$
DECLARE
maxdist ALIAS FOR $16;
sourcesrid integer;
newnodatavalue float8;
newinitialvalue float8;
newrast raster;
newpixeltype text;
BEGIN
-- Assuming source point table has one SRID for all geometries
EXECUTE 'SELECT DISTINCT ST_SRID(' || quote_ident(sourcegeomcolumn) || ') FROM ' || quote_ident(sourceschema) || '.' || quote_ident(sourcetable) || ';' INTO sourcesrid;
-- Create a new empty raster using provided georeference
newrast := ST_MakeEmptyRaster(width, height, rastulx, rastuly, rastscalex,rastscaley, rastskewx, rastskewy, rastsrid);
-- Check if the new raster is empty (width = 0 OR height = 0)
IF ST_IsEmpty(newrast) THEN
RAISE NOTICE 'ST_EuclideanDistance: Raster is empty. Returning an empty raster';
RETURN newrast;
END IF;
-- Set the new pixeltype
newpixeltype := pixeltype;
IF newpixeltype = '1BB' OR newpixeltype = '2BUI' OR (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_EuclideanDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
END IF;
-- Check for notada value
newnodatavalue := rastndv;
IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN
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';
newnodatavalue := ST_MinPossibleValue(newpixeltype);
END IF;
-- We set the initial value of the resulting raster to nodata value.
-- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue
newinitialvalue := newnodatavalue;
-- Create the raster receiving all the computed distance values. Initialize it to the new initial value.
newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
-- Call ST_MapAlgebraFct() with custom EuclideanDistance function to set the new raster to the resulted distance raster
newrast := ST_MapAlgebraFct(newrast,NULL,'euclidean_distance_fct(float,integer[],text[])'::regprocedure,
width::text,height::text,
rastulx::text, rastuly::text,
rastscalex::text,rastscaley::text,
rastskewx::text, rastskewy::text,
rastsrid::text, newnodatavalue::text,
sourceschema,sourcetable,sourcegeomcolumn,
snaptocentroid::text,maxdist::text);
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
--------------------------------------------------------------------
-- ST_EuclideanDistance - Return a raster generated from a one-band reference raster.
-- Pixel values are the Euclidean distance from pixel to source points
-- which is a table of point geometries
-- (points for now, lines or polygons in the future).
-- Arguments
-- refrast raster - Reference raster align with which Euclidean distance
-- raster will be created.
-- rastsrid - SRID of the resulted raster
-- pixeltype - Pixeltype assigned to the resulting raster.
-- As distance raster, we only accept pixel types of:
-- '4BUI','8BSI','8BUI','16BSI','16BUI','32BSI','32BUI','32BF','64BF'.
-- sourceschema - Database schema of the source geometry table.
-- sourcetable - Database table name of the source geometry table.
-- sourcegeomcolumn - Database geometry column name of the source geometry table.
-- snaptocentroid - When True, the option will "snap" any source point geometry to the centroid of
-- the pixel intersecting with it, resulting in a distance value equal to zero.
-- The Euclidean distance is calculated from the centroid of each pixel to the
-- centroid of the pixel intersecting with its nearest source point geometry.
-- maxdist - Maximum distance from pixel to the source when the source is too far away from the pixel.
-- Pixel that exceeds maximum distance will be assigned nodatavalue.
-- By default to be -1 if not specified.
--------------------------------------------------------------------
DROP FUNCTION IF EXISTS ST_EuclideanDistance(raster,text,text,text,text,boolean,float8);
CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster,
pixeltype text,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
snaptocentroid boolean,
float8 = -1)
RETURNS raster
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) $$
LANGUAGE 'sql';
Cost Distance
ST_CostDistance()
--------------------------------------------------------------------
-- ST_CostDistance - Return a raster generated from a one band reference
-- raster a table and a one band cost raster, which values
-- are the the least accumulative cost of traversing
-- the cost surface from source points to each cell.
-- (points for now, lines or polygons in the future).
-- Arguments
-- refrast raster - Reference raster align with which cost distance
-- raster will be created.
-- pixeltype text - Pixeltype assigned to the resulting raster. By default
-- to be the pixeltype of the reference raster if not specified.
-- costrast raster - Cost raster in which each cell value represents the
-- cost of traversing this cell.
-- sourceschema text - Database schema of the source geometry table.
-- sourcetable text - Source geometry table name.
-- sourcegeomcolumn text - Source geometry table geometry column name.
-- maxdist double precision - Maximum distance from pixel to the source
-- when the source is too far from the pixel.
-- Pixel that exceeds it will be assigned nodatavalue.
--------------------------------------------------------------------
-- Create array type which stores the column, row, and cost value of each pixel in the iteration.
CREATE TYPE p_node AS (id integer,x integer,y integer,val float8);
-- Function to get a set of elements from an array
DROP FUNCTION IF EXISTS explode_array(p_node[]);
CREATE OR REPLACE FUNCTION explode_array(p_node[])
RETURNS SETOF p_node AS
$$
SELECT $1[i] FROM GENERATE_SERIES(1,ARRAY_UPPER($1,1)) AS i;
$$
LANGUAGE 'sql' IMMUTABLE;
DROP FUNCTION IF EXISTS ST_CostDistance(refrast raster,
pixeltype text,
costrast raster,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1);
CREATE OR REPLACE FUNCTION ST_CostDistance(refrast raster,
pixeltype text,
costrast raster,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1)
RETURNS raster AS
$$
DECLARE
maxdist ALIAS FOR $7;
newrast raster;
newwidth integer;
newheight integer;
newsrid integer;
newpixeltype text;
newnodatavalue float8;
newinitialvalue float8;
x integer;
y integer;
sourcex float8;
sourcey float8;
sourcexr integer;
sourceyr integer;
sourcesrid integer;
sourceboxgeom geometry;
costarray p_node[] = '{}';
node p_node;
costval float8;
-- Assuming reference raster has only one band
band integer := 1;
nx integer;
ny integer;
ncoordx float8;
ncoordy float8;
ncostval float8;
nnewcostval float8;
ncostrastx integer;
ncostrasty integer;
ncostrastval float8;
arrayid integer := 0;
tempnode p_node[];
minval float8;
minnodeid int;
exesql text;
geom geometry;
flag boolean := False;
BEGIN
-- Check if reference raster is NULL
IF refrast IS NULL THEN
RAISE NOTICE 'ST_CostDistance: Reference raster is NULL. Returning NULL';
RETURN NULL;
END IF;
-- Check if cost raster is NULL
IF costrast IS NULL THEN
RAISE NOTICE 'ST_CostDistance: Cost raster is NULL. Returning NULL';
RETURN NULL;
END IF;
newwidth := ST_Width(refrast);
newheight := ST_Height(refrast);
newsrid := ST_SRID(refrast);
-- Assuming source point table has one SRID for all geometries
EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid;
-- Create a new empty raster having the same georeference as the reference raster
newrast := ST_MakeEmptyRaster(newwidth, newheight, ST_UpperLeftX(refrast), ST_UpperLeftY(refrast), ST_ScaleX(refrast), ST_ScaleY(refrast), ST_SkewX(refrast), ST_SkewY(refrast), newsrid);
-- Check if the new raster is empty (width = 0 OR height = 0)
IF ST_IsEmpty(newrast) THEN
RAISE NOTICE 'ST_CostDistance: Raster is empty. Returning an empty raster';
RETURN newrast;
END IF;
-- Set the new pixeltype
newpixeltype := pixeltype;
-- If pixeltype not specified then use the pixeltype of the reference raster
IF newpixeltype IS NULL THEN
newpixeltype := ST_BandPixelType(refrast, 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_CostDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
END IF;
-- Check for notada value
newnodatavalue := ST_BandNodataValue(refrast, band);
IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN
RAISE NOTICE 'ST_CostDistance: 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';
newnodatavalue := ST_MinPossibleValue(newpixeltype);
END IF;
-- We set the initial value of the resulting raster to nodata value.
-- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue
newinitialvalue := newnodatavalue;
-- Create the raster receiving all the computed cost distance values. Initialize it to the new initial value.
newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
-- Check if source points boundary intersects with extent of cost raster, if not, return new raster initiated with no data value.
exesql := 'SELECT ST_CONVEXHULL(ST_COLLECT(' || sourcegeomcolumn || ')) FROM ' || sourceschema || '.' || sourcetable || ';';
EXECUTE exesql INTO sourceboxgeom;
SELECT ST_Intersects(costrast,sourceboxgeom) INTO flag;
IF NOT flag THEN
RAISE NOTICE 'ST_CostDistance: Source points are out side of the cost raster. Returning raster with no data values.';
RETURN newrast;
END IF;
-- Set raster pixel value to nodatavalue or zero if pixel is source point;
-- Initialize cost_array with source points that are within the extent of cost raster.
exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';';
FOR geom IN EXECUTE(exesql) LOOP
sourcex := ST_X(geom);
sourcey := ST_Y(geom);
sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey);
sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey);
-- If source point within new raster extent
IF sourcexr <= width AND sourceyr <= height THEN
-- If source point within cost raster extent, set new raster pixel value to zero
-- since there is no accumulative cost to return to itself;
-- Add source point to cost array.
IF ST_Within(ST_SetSRID(ST_MakePoint(sourcex,sourcey),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
arrayid := arrayid + 1;
node.id := arrayid;
node.x := sourcexr;
node.y := sourceyr;
node.val := costval;
costarray := ARRAY_APPEND(costarray,node);
END IF;
END IF;
END LOOP;
-- Loop through the cost array, pop up node with least accumulative cost value;
-- Calculate cost values for its neighbors (8 directions);
flag := True;
WHILE flag LOOP
-- End loop when cost values in the cost array are all NULL
SELECT ARRAY_AGG(cv) FROM (SELECT DISTINCT explode_array(costarray) AS cv) AS foo INTO tempnode;
IF ARRAY_LENGTH(tempnode,1) = 1 AND tempnode[1] = NULL THEN
flag := False;
ELSE
SELECT MIN(costarray[i].val) FROM GENERATE_SERIES(1,ARRAY_UPPER(costarray,1)) AS i INTO minval;
SELECT costarray[i].id FROM GENERATE_SERIES(1,ARRAY_UPPER(costarray,1)) AS i WHERE costarray[i].val = minval INTO minnodeid;
-- Pop up node with min cost value, set its val in the array to be NULL
x := costarray[minnodeid].x;
y := costarray[minnodeid].y;
costval := ST_Value(newrast,band,x,y);
costarray[minnodeid] := NULL;
-- Calculate cost distance for neighbors in perpendicular directions
-- neighbor(x,y+1)
nx := x;
ny := y+1;
IF ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x,y-1)
nx := x;
ny := y-1;
IF ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
-- If max distance specified
IF maxdist != -1 THEN
-- Any distance exceeds max distance is assigned nodata value
IF nnewcostval > maxdist THEN
nnewcostval := newnodatavalue;
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
if nnewcostval <> newnodatavalue THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y)
nx := x+1;
ny := y;
IF nx <= newwidth THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x-1,y)
nx := x-1;
ny := y;
IF nx >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- Calculate cost distance for neighbors in diagonal directions
-- neighbor(x-1,y+1)
nx := x-1;
ny := y+1;
IF nx >= 1 AND ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x-1,y-1)
nx := x-1;
ny := y-1;
IF nx >= 1 AND ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y+1)
nx := x+1;
ny := y+1;
IF nx <= newwidth AND ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y-1)
nx := x+1;
ny := y-1;
IF nx <= newwidth AND ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
END IF;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
--------------------------------------------------------------------
-- ST_CostDistance - Return a raster generated from a set of raster specifications
-- raster a table and a one band cost raster, which values
-- are the the least accumulative cost of traversing
-- the cost surface from source points to each cell.
-- (points for now, lines or polygons in the future).
-- Arguments
-- width,height,rastulx,rastuly,rastscalex,rastscaley,rastskewx,rastskewy
-- - Dimenssions and georeference of resulted raster.
-- rastsrid - SRID of the resulted raster
-- rastndv - Nodata value for resulted raster
-- pixeltype text - Pixeltype assigned to the resulting raster. By default
-- to be the pixeltype of the reference raster if not specified.
-- costrast raster - Cost raster in which each cell value represents the
-- cost of traversing this cell.
-- sourceschema text - Database schema of the source geometry table.
-- sourcetable text - Source geometry table name.
-- sourcegeomcolumn text - Source geometry table geometry column name.
-- maxdist double precision - Maximum distance from pixel to the source
-- when the source is too far from the pixel.
-- Pixel that exceeds it will be assigned nodatavalue.
--------------------------------------------------------------------
DROP FUNCTION IF EXISTS ST_CostDistance(width integer,
height integer,
rastulx float8,
rastuly float8,
rastscalex float8,
rastscaley float8,
rastskewx float8,
rastskewy float8,
rastsrid integer,
rastndv float8,
pixeltype text,
costrast raster,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1);
CREATE OR REPLACE FUNCTION ST_CostDistance(width integer,
height integer,
rastulx float8,
rastuly float8,
rastscalex float8,
rastscaley float8,
rastskewx float8,
rastskewy float8,
rastsrid integer,
rastndv float8,
pixeltype text,
costrast raster,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1)
RETURNS raster AS
$$
DECLARE
maxdist ALIAS FOR $7;
newrast raster;
newwidth integer;
newheight integer;
newsrid integer;
newpixeltype text;
newnodatavalue float8;
newinitialvalue float8;
x integer;
y integer;
sourcex float8;
sourcey float8;
sourcexr integer;
sourceyr integer;
sourcesrid integer;
sourceboxgeom geometry;
costarray p_node[] = '{}';
node p_node;
costval float8;
-- Assuming reference raster has only one band
band integer := 1;
nx integer;
ny integer;
ncoordx float8;
ncoordy float8;
ncostval float8;
nnewcostval float8;
ncostrastx integer;
ncostrasty integer;
ncostrastval float8;
arrayid integer := 0;
tempnode p_node[];
minval float8;
minnodeid int;
exesql text;
geom geometry;
flag boolean := False;
BEGIN
-- Check if cost raster is NULL
IF costrast IS NULL THEN
RAISE NOTICE 'ST_CostDistance: Cost raster is NULL. Returning NULL';
RETURN NULL;
END IF;
newwidth := width;
newheight := height;
newsrid := rastsrid;
-- Assuming source point table has one SRID for all geometries
EXECUTE 'SELECT DISTINCT ST_SRID(' || sourcegeomcolumn || ') FROM ' || sourceschema || '.' || sourcetable || ';' INTO sourcesrid;
-- Create a new empty raster having the same georeference as the reference raster
newrast := ST_MakeEmptyRaster(width, height, rastulx, rastuly, rastscalex,rastscaley, rastskewx, rastskewy, rastsrid);
-- Check if the new raster is empty (width = 0 OR height = 0)
IF ST_IsEmpty(newrast) THEN
RAISE NOTICE 'ST_CostDistance: Raster is empty. Returning an empty raster';
RETURN newrast;
END IF;
-- Set the new pixeltype
newpixeltype := pixeltype;
-- If pixeltype not specified then use the pixeltype of the reference raster
IF newpixeltype IS NULL THEN
newpixeltype := ST_BandPixelType(refrast, 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_CostDistance: Invalid pixeltype "%". Aborting.', newpixeltype;
END IF;
-- Check for notada value
newnodatavalue := ST_BandNodataValue(refrast, band);
IF newnodatavalue IS NULL OR newnodatavalue < ST_MinPossibleValue(newpixeltype) OR newnodatavalue > (-ST_MinPossibleValue(newpixeltype) - 1) THEN
RAISE NOTICE 'ST_CostDistance: 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';
newnodatavalue := ST_MinPossibleValue(newpixeltype);
END IF;
-- We set the initial value of the resulting raster to nodata value.
-- If nodatavalue is null then the raster will be initialise to ST_MinPossibleValue
newinitialvalue := newnodatavalue;
-- Create the raster receiving all the computed cost distance values. Initialize it to the new initial value.
newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
-- Check if source points boundary intersects with extent of cost raster, if not, return new raster initiated with no data value.
exesql := 'SELECT ST_CONVEXHULL(ST_COLLECT(' || sourcegeomcolumn || ')) FROM ' || sourceschema || '.' || sourcetable || ';';
EXECUTE exesql INTO sourceboxgeom;
SELECT ST_Intersects(costrast,sourceboxgeom) INTO flag;
IF NOT flag THEN
RAISE NOTICE 'ST_CostDistance: Source points are out side of the cost raster. Returning raster with no data values.';
RETURN newrast;
END IF;
-- Set raster pixel value to nodatavalue or zero if pixel is source point;
-- Initialize cost_array with source points that are within the extent of cost raster.
exesql := 'SELECT ' || sourcegeomcolumn || ' FROM ' || sourceschema || '.' || sourcetable || ';';
FOR geom IN EXECUTE(exesql) LOOP
sourcex := ST_X(geom);
sourcey := ST_Y(geom);
sourcexr := ST_World2RasterCoordX(newrast, sourcex, sourcey);
sourceyr := ST_World2RasterCoordY(newrast, sourcex, sourcey);
-- If source point within new raster extent
IF sourcexr <= width AND sourceyr <= height THEN
-- If source point within cost raster extent, set new raster pixel value to zero
-- since there is no accumulative cost to return to itself;
-- Add source point to cost array.
IF ST_Within(ST_SetSRID(ST_MakePoint(sourcex,sourcey),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
arrayid := arrayid + 1;
node.id := arrayid;
node.x := sourcexr;
node.y := sourceyr;
node.val := costval;
costarray := ARRAY_APPEND(costarray,node);
END IF;
END IF;
END LOOP;
-- Loop through the cost array, pop up node with least accumulative cost value;
-- Calculate cost values for its neighbors (8 directions);
flag := True;
WHILE flag LOOP
-- End loop when cost values in the cost array are all NULL
SELECT ARRAY_AGG(cv) FROM (SELECT DISTINCT explode_array(costarray) AS cv) AS foo INTO tempnode;
IF ARRAY_LENGTH(tempnode,1) = 1 AND tempnode[1] = NULL THEN
flag := False;
ELSE
SELECT MIN(costarray[i].val) FROM GENERATE_SERIES(1,ARRAY_UPPER(costarray,1)) AS i INTO minval;
SELECT costarray[i].id FROM GENERATE_SERIES(1,ARRAY_UPPER(costarray,1)) AS i WHERE costarray[i].val = minval INTO minnodeid;
-- Pop up node with min cost value, set its val in the array to be NULL
x := costarray[minnodeid].x;
y := costarray[minnodeid].y;
costval := ST_Value(newrast,band,x,y);
costarray[minnodeid] := NULL;
-- Calculate cost distance for neighbors in perpendicular directions
-- neighbor(x,y+1)
nx := x;
ny := y+1;
IF ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x,y-1)
nx := x;
ny := y-1;
IF ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
-- If max distance specified
IF maxdist != -1 THEN
-- Any distance exceeds max distance is assigned nodata value
IF nnewcostval > maxdist THEN
nnewcostval := newnodatavalue;
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
if nnewcostval <> newnodatavalue THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y)
nx := x+1;
ny := y;
IF nx <= newwidth THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x-1,y)
nx := x-1;
ny := y;
IF nx >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- Calculate cost distance for neighbors in diagonal directions
-- neighbor(x-1,y+1)
nx := x-1;
ny := y+1;
IF nx >= 1 AND ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x-1,y-1)
nx := x-1;
ny := y-1;
IF nx >= 1 AND ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y+1)
nx := x+1;
ny := y+1;
IF nx <= newwidth AND ny <= newheight THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
-- neighbor(x+1,y-1)
nx := x+1;
ny := y-1;
IF nx <= newwidth AND ny >= 1 THEN
ncoordx := ST_Raster2WorldCoordX(newrast, x, y);
ncoordy := ST_Raster2WorldCoordY(newrast, x, y);
IF ST_Within(ST_SetSRID(ST_MakePoint(ncoordx,ncoordy),ST_SRID(costrast)),ST_Envelope(costrast)) THEN
ncostrastx := ST_World2RasterCoordX(costrast,ncoordx,ncoordy);
ncostrasty := ST_World2RasterCoordY(costrast,ncoordx,ncoordy);
ncostrastval := ST_Value(costrast,band,ncostrastx,ncostrasty);
nnewcostval := sqrt(2) * (costval + ncostrastval)/2;
ncostval := ST_Value(newrast,band,nx,ny);
IF ncostval <> newnodatavalue THEN
IF nnewcostval < ncostval THEN
ncostval := nnewcostval;
arrayid := arrayid + 1;
node.id := arrayid;
node.x := nx;
node.y := ny;
node.val := ncostval;
costarray := ARRAY_APPEND(costarray,node);
newrast := ST_SetValue(newrast, band, nx, ny, ncostval);
END IF;
END IF;
ELSE
newrast := ST_SetValue(newrast, band, nx, ny, newnodatavalue);
END IF;
END IF;
END IF;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
Last modified
12 years ago
Last modified on 08/28/12 10:39:32
Note:
See TracWiki
for help on using the wiki.
