| Version 1 (modified by , 12 years ago) ( diff ) |
|---|
Euclidean Distance
ST_EuclideanDistance()
--------------------------------------------------------------------
-- ST_EuclideanDistance - (for one raster) Return a raster generated from
-- a set of raster specifications, which values are
-- the Euclidean distance from pixel to the points
-- (for now only concerning point geometries)
-- from the input source table of 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 text - Pixeltype assigned to the resulting raster. By default
-- to be the pixeltype of the reference raster if not specified.
-- 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_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,
double precision = -1);
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,
double precision = -1)
RETURNS raster AS
$$
DECLARE
maxdist ALIAS FOR $15;
x integer;
y integer;
sourcex float8;
sourcey float8;
sourcexr integer;
sourceyr integer;
sourcesrid integer;
newx float8;
newy float8;
exesql text;
newnodatavalue float8;
newinitialvalue float8;
newrast raster;
newval float8;
newpixeltype text;
-- Assuming reference raster has only one band
band integer := 1;
geom geometry;
BEGIN
-- 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 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' 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_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);
-- Set raster pixel value to zero if pixel is source point
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 raster range, then set pixel value to zero
IF sourcexr <= width AND sourceyr <= height THEN
newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
END IF;
END LOOP;
-- Scan pixels in the raster set pixel value to the computed Eculidean distance
-- to its nearest source point.
-- There is place for optimization here for a more efficient scanning method
FOR x IN 1..width LOOP
FOR y IN 1..height LOOP
newx := ST_Raster2WorldCoordX(newrast, x, y);
newy := ST_Raster2WorldCoordY(newrast, x, y);
-- If pixel is not the source point then compute Eculidean distance
IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN
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));';
EXECUTE exesql INTO newval;
-- If max distance specified
IF maxdist != -1 THEN
-- Any distance exceeds max distance is assigned nodata value
IF newval > maxdist THEN
newval = newnodatavalue;
END IF;
END IF;
END IF;
newrast := ST_SetValue(newrast, band, x, y, newval);
END LOOP;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
-- test
-- 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);
-- 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);
--------------------------------------------------------------------
-- ST_EuclideanDistance - (for one raster) Return a raster generated
-- from a one band reference raster,
-- which values are the Euclidean distance
-- from pixel to the points
-- (for now only concerning point geometries)
-- from the input source table of geometries
-- (points for now, lines or polygons in the future).
-- Arguments
-- refrast raster - Reference raster align with which Euclidean 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.
-- 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_EuclideanDistance(refrast raster,
pixeltype text,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1);
CREATE OR REPLACE FUNCTION ST_EuclideanDistance(refrast raster,
pixeltype text,
sourceschema text,
sourcetable text,
sourcegeomcolumn text,
double precision = -1)
RETURNS raster AS
$$
DECLARE
maxdist ALIAS FOR $6;
width integer;
height integer;
x integer;
y integer;
sourcex float8;
sourcey float8;
sourcexr integer;
sourceyr integer;
sourcesrid integer;
newsrid integer;
newx float8;
newy float8;
exesql text;
newnodatavalue float8;
newinitialvalue float8;
newrast raster;
newval float8;
newpixeltype text;
-- Assuming reference raster has only one band
band integer := 1;
geom geometry;
BEGIN
-- Check if reference raster is NULL
IF refrast IS NULL THEN
RAISE NOTICE 'ST_EuclideanDistance: Reference raster is NULL. Returning NULL';
RETURN NULL;
END IF;
width := ST_Width(refrast);
height := 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(width, height, 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_EuclideanDistance: 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_EuclideanDistance: 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_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);
-- Set raster pixel value to zero if pixel is source point
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 raster range, then set pixel value to zero
IF sourcexr <= width AND sourceyr <= height THEN
newrast := ST_SetValue(newrast, band, sourcexr, sourceyr, 0);
END IF;
END LOOP;
-- Scan pixels in the raster set pixel value to the computed Eculidean distance
-- to its nearest source point.
-- There is place for optimization here for a more efficient scanning method
FOR x IN 1..width LOOP
FOR y IN 1..height LOOP
newx := ST_Raster2WorldCoordX(newrast, x, y);
newy := ST_Raster2WorldCoordY(newrast, x, y);
-- If pixel is not the source point then compute Eculidean distance
IF ST_Value(newrast, band, x, y) != 0 OR ST_Value(newrast, band, x, y) IS NULL THEN
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));';
EXECUTE exesql INTO newval;
-- If max distance specified
IF maxdist != -1 THEN
-- Any distance exceeds max distance is assigned nodata value
IF newval > maxdist THEN
newval := newnodatavalue;
END IF;
END IF;
END IF;
newrast := ST_SetValue(newrast, band, x, y, newval);
END LOOP;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
-- test
-- CREATE TABLE test_eudist AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom') AS rast FROM test_raster);
-- CREATE TABLE test_eudist2 AS (SELECT 1 AS rid,ST_EuclideanDistance(rast,NULL,'public','test_geom','the_geom',4954880) AS rast FROM test_raster);
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';
Note:
See TracWiki
for help on using the wiki.
