self-hosted/postgis
self-hosted/postgis
1
0
Fork 0
mirror of https://git.osgeo.org/gitea/postgis/postgis synced 2024-10-24 09:02:37 +00:00
Code Issues Packages Projects Releases Wiki Activity

--General enhancement to the script functions. Still a lot of work to do...

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@6147 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Pierre Racine 2010-10-28 23:08:17 +00:00
parent a0415e4f38
commit 6c1adcafe6
12 changed files with 1715 additions and 186 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

401
raster/scripts/plpgsql/_MapAlgebraParts.sql Normal file
View file

@ -0,0 +1,401 @@
----------------------------------------------------------------------
--
-- $Id: _MapAlgebraParts.sql 6127 2010-10-25 16:06:00Z jorgearevalo $
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- Note: The functions found in this file are for exclusive usage of ST_MapAlgebra2
CREATE OR REPLACE FUNCTION max(a int, b int)
RETURNS int
AS 'SELECT CASE WHEN $1 < $2 THEN $2 ELSE $1 END'
LANGUAGE 'SQL' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION min(a int, b int)
RETURNS int
AS 'SELECT CASE WHEN $1 < $2 THEN $1 ELSE $2 END'
LANGUAGE 'SQL' IMMUTABLE STRICT;
--DROP FUNCTION _MapAlgebraParts(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int);
CREATE OR REPLACE FUNCTION _MapAlgebraParts(r1x int,
r1y int,
r1w int,
r1h int,
r2x int,
r2y int,
r2w int,
r2h int)
RETURNS int[]
AS $$
DECLARE
z11x int;
z11y int;
z11w int;
z11h int;
z12x int;
z12y int;
z12w int;
z12h int;
z13x int;
z13y int;
z13w int;
z13h int;
z14x int;
z14y int;
z14w int;
z14h int;
z21x int;
z21y int;
z21w int;
z21h int;
z22x int;
z22y int;
z22w int;
z22h int;
z23x int;
z23y int;
z23w int;
z23h int;
z24x int;
z24y int;
z24w int;
z24h int;
zcx int;
zcy int;
zcw int;
zch int;
BEGIN
z11x := r1x;
z11y := r1y;
z11w := r1w;
z11h := min(max(0, r2y - r1y), r1h);
z12x := r1x;
z12y := z11y + z11h;
z12w := min(max(0, r2x - r1x), r1w);
z12h := max(0, min(max(0, r2y + r2h - (r1y + z11h)), z11y + r1h - z12y));
z13x := max(min(r1x + r1w, r2x + r2w), r1x);
z13y := z12y;
z13w := min(max(0, r1x + r1w - (r2x + r2w)), r1w);
z13h := z12h;
z14x := r1x;
z14y := z12y + z12h;
z14w := r1w;
z14h := min(max(0, r1y + r1h - (r2y + r2h)), r1h);
z21x := r2x;
z21y := r2y;
z21w := r2w;
z21h := min(max(0, r1y - r2y), r2h);
z22x := r2x;
z22y := z21y + z21h;
z22w := min(max(0, r1x - r2x), r2w);
z22h := max(0, min(max(0, r1y + r1h - (r2y + z21h)), z21y + r2h - z22y));
z23x := max(min(r2x + r2w, r1x + r1w), r2x);
z23y := z22y;
z23w := min(max(0, r2x + r2w - (r1x + r1w)), r2w);
z23h := z22h;
z24x := r2x;
z24y := z22y + z22h;
z24w := r2w;
z24h := min(max(0, r2y + r2h - (r1y + r1h)), r2h);
zcx := z12x + z12w;
zcy := z12y;
zcw := z13x - (z12x + z12w);
zch := z14y - z12y;
-- Merge z11 with z12 if they are continuous parts of the same vertical bar
IF z12h > 0 AND z12x = z11x AND z12w = z11w THEN
z12y := z11y;
z12h := z11h + z12h;
z11h := 0;
END IF;
-- Merge z11 with z13 if they are continuous parts of the same vertical bar
IF z13h > 0 AND z13x = z11x AND z13w = z11w THEN
z13y := z11y;
z13h := z11h + z13h;
z11h := 0;
END IF;
-- Merge z12 with z14 if they are continuous parts of the same vertical bar
IF z14h > 0 AND z14x = z12x AND z14w = z12w THEN
z14y := z12y;
z14h := z12h + z14h;
z12h := 0;
END IF;
-- Merge z13 with z14 if they are continuous parts of the same vertical bar
IF z14h > 0 AND z14x = z13x AND z14w = z13w THEN
z14y := z13y;
z14h := z13h + z14h;
z13h := 0;
END IF;
-- Merge z21 with z22 if they are continuous parts of the same vertical bar
IF z22h > 0 AND z22x = z21x AND z22w = z21w THEN
z22y := z21y;
z22h := z21h + z22h;
z21h := 0;
END IF;
-- Merge z21 with z23 if they are continuous parts of the same vertical bar
IF z23h > 0 AND z23x = z21x AND z23w = z21w THEN
z23y := z21y;
z23h := z21h + z23h;
z21h := 0;
END IF;
-- Merge z22 with z24 if they are continuous parts of the same vertical bar
IF z24h > 0 AND z24x = z22x AND z24w = z22w THEN
z24y := z22y;
z24h := z22h + z24h;
z22h := 0;
END IF;
-- Merge z23 with z24 if they are continuous parts of the same vertical bar
IF z24h > 0 AND z24x = z23x AND z24w = z23w THEN
z24y := z23y;
z24h := z23h + z24h;
z23h := 0;
END IF;
RETURN ARRAY[z11x, z11y, z11w, z11h, z12x, z12y, z12w, z12h, z13x, z13y, z13w, z13h, z14x, z14y, z14w, z14h, z21x, z21y, z21w, z21h, z22x, z22y, z22w, z22h, z23x, z23y, z23w, z23h, z24x, z24y, z24w, z24h, zcx, zcy, zcw, zch];
END;
$$
LANGUAGE 'plpgsql';
--DROP FUNCTION _MapAlgebraPartsGeom(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int);
CREATE OR REPLACE FUNCTION _MapAlgebraPartsGeom(nx int,
ny int,
r1x int,
r1y int,
r1w int,
r1h int,
r2x int,
r2y int,
r2w int,
r2h int)
RETURNS SETOF geometry AS
$$
DECLARE
BEGIN
RETURN NEXT ST_MakeBox2D(ST_Point(10 * ny + r1x, -10 * nx + 5 - r1y), ST_Point(10 * ny + r1x + r1w, -10 * nx + 5 - (r1y + r1h)))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(10 * ny + r2x, -10 * nx + 5 - r2y), ST_Point(10 * ny + r2x + r2w, -10 * nx + 5 - (r2y + r2h)))::geometry;
RETURN;
END;
$$
LANGUAGE 'plpgsql';
CREATE OR REPLACE FUNCTION _MapAlgebraAllPartsGeom(r1x int,
r1y int,
r1w int,
r1h int,
r2x int,
r2y int,
r2w int,
r2h int)
RETURNS SETOF geometry AS
$$
DECLARE
z int[];
BEGIN
z := _MapAlgebraParts(r1x, r1y, r1w, r1h, r2x, r2y, r2w, r2h);
RETURN NEXT ST_MakeBox2D(ST_Point(z[1], z[2]), ST_Point(z[1] + z[3], z[2] + z[4]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[5], z[6]), ST_Point(z[5] + z[7], z[6] + z[8]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[9], z[10]), ST_Point(z[9] + z[11], z[10] + z[12]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[13], z[14]), ST_Point(z[13] + z[15], z[14] + z[16]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[17], z[18]), ST_Point(z[17] + z[19], z[18] + z[20]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[21], z[22]), ST_Point(z[21] + z[23], z[22] + z[24]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[25], z[26]), ST_Point(z[25] + z[27], z[26] + z[28]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[29], z[30]), ST_Point(z[29] + z[31], z[30] + z[32]))::geometry;
RETURN NEXT ST_MakeBox2D(ST_Point(z[33], z[34]), ST_Point(z[33] + z[35], z[34] + z[36]))::geometry;
RETURN;
END;
$$
LANGUAGE 'plpgsql';
SELECT asbinary(_MapAlgebraAllPartsGeom(0, 0, 1, 1, 1, 0, 1, 1))
CREATE OR REPLACE FUNCTION X1W1X2W2()
RETURNS SETOF record AS
$$
DECLARE
x1w1x2w2 record;
BEGIN
x1w1x2w2 := (0, 3-4, 2, 0-4, 2);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (1, 2-4, 3, 0-4, 2);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (2, 2-4, 3, 0-4, 3);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (3, 2-4, 3, 0-4, 5);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (4, 1-4, 3, 0-4, 5);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (5, 1-4, 4, 1-4, 2);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (6, 1-4, 3, 1-4, 3);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (7, 1-4, 2, 1-4, 4);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (8, 0-4, 5, 1-4, 3);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (9, 0-4, 5, 1-4, 4);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (10, 0-4, 3, 2-4, 3);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (11, 0-4, 3, 3-4, 2);
RETURN NEXT x1w1x2w2;
x1w1x2w2 := (12, 0-4, 2, 3-4, 2);
RETURN NEXT x1w1x2w2;
RETURN;
END;
$$
LANGUAGE 'plpgsql';
CREATE OR REPLACE FUNCTION Y1H1Y2H2()
RETURNS SETOF record AS
$$
DECLARE
y1h1y2h2 record;
BEGIN
y1h1y2h2 := (0, 3-4, 2, 0-4, 2);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (1, 2-4, 3, 0-4, 2);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (2, 2-4, 3, 0-4, 3);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (3, 2-4, 3, 0-4, 5);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (4, 1-4, 3, 0-4, 5);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (5, 1-4, 4, 1-4, 2);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (6, 1-4, 3, 1-4, 3);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (7, 1-4, 2, 1-4, 4);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (8, 0-4, 5, 1-4, 3);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (9, 0-4, 5, 1-4, 4);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (10, 0-4, 3, 2-4, 3);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (11, 0-4, 3, 3-4, 2);
RETURN NEXT y1h1y2h2;
y1h1y2h2 := (12, 0-4, 2, 3-4, 2);
RETURN NEXT y1h1y2h2;
RETURN;
END;
$$
LANGUAGE 'plpgsql';
_MapAlgebraParts(r1x, r1y, r1w, r1h, r2x, r2y, r2w, r2h)
SELECT nx, x1, w1, x2, w2 FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int);
SELECT nx, ny, x1, w1, x2, w2, y1, h1, y2, h2
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int);
SELECT nx, ny, x1, w1, x2, w2, y1, h1, y2, h2, asbinary(_MapAlgebraPartsGeom(nx, ny, x1, y1, w1, h1, x2, y2, w2, h2))
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int);
SELECT asbinary(_MapAlgebraPartsGeom(nx, ny, x1, y1, w1, h1, x2, y2, w2, h2))
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int);
-- First series of zones covering raster 1
SELECT nx, ny, map[1], map[2], map[3], map[4],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[1], -10 * nx + 5 - map[2]), ST_Point(10 * ny + map[1] + map[3], -10 * nx + 5 - (map[2] + map[4])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Second series of zones covering raster 1
SELECT nx, ny, map[5], map[6], map[7], map[8],
asbinary(ST_Point(10 * ny + map[5], -10 * nx + 5 - map[6])::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
SELECT nx, ny, map[5], map[6], map[7], map[8],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[5], -10 * nx + 5 - map[6]), ST_Point(10 * ny + map[5] + map[7], -10 * nx + 5 - (map[6] + map[8])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Third series of zones covering raster 1
SELECT nx, ny, map[9], map[10], map[11], map[12],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[9], -10 * nx + 5 - map[10]), ST_Point(10 * ny + map[9] + map[11], -10 * nx + 5 - (map[10] + map[12])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Fourth series of zones covering raster 1
SELECT nx, ny, map[13], map[14], map[15], map[16],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[13], -10 * nx + 5 - map[14]), ST_Point(10 * ny + map[13] + map[15], -10 * nx + 5 - (map[14] + map[16])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- First series of zones covering raster 2
SELECT nx, ny, map[17], map[18], map[19], map[20],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[17], -10 * nx + 5 - map[18]), ST_Point(10 * ny + map[17] + map[19], -10 * nx + 5 - (map[18] + map[20])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Second series of zones covering raster 2
SELECT nx, ny, map[21], map[22], map[23], map[24],
asbinary(ST_Point(10 * ny + map[21], -10 * nx + 5 - map[22])::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
SELECT nx, ny, map[21], map[22], map[23], map[24],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[21], -10 * nx + 5 - map[22]), ST_Point(10 * ny + map[21] + map[23], -10 * nx + 5 - (map[22] + map[24])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Third series of zones covering raster 2
SELECT nx, ny, map[25], map[26], map[27], map[28],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[25], -10 * nx + 5 - map[26]), ST_Point(10 * ny + map[25] + map[27], -10 * nx + 5 - (map[26] + map[28])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Fourth series of zones covering raster 2
SELECT nx, ny, map[29], map[30], map[31], map[32],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[28], -10 * nx + 5 - map[29]), ST_Point(10 * ny + map[28] + map[30], -10 * nx + 5 - (map[29] + map[31])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;
-- Common zone
SELECT nx, ny, map[33], map[34], map[35], map[36],
asbinary(ST_MakeBox2D(ST_Point(10 * ny + map[33], -10 * nx + 5 - map[34]), ST_Point(10 * ny + map[33] + map[35], -10 * nx + 5 - (map[34] + map[36])))::geometry)
FROM (
SELECT nx, ny, x1, y1, w1, h1, x2, y2, w2, h2, _MapAlgebraParts(x1, y1, w1, h1, x2, y2, w2, h2) as map
FROM X1W1X2W2() as (nx int, x1 int, w1 int, x2 int, w2 int), Y1H1Y2H2() as (ny int, y1 int, h1 int, y2 int, h2 int)
) as foo;

54
raster/scripts/plpgsql/st_addband.sql Normal file
View file

@ -0,0 +1,54 @@
----------------------------------------------------------------------
--
-- $Id: st_addband.sql 6127 2010-10-25 16:06:00Z jorgearevalo $
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_AddBand function found in this file still many improvements before being implemented in C.
CREATE OR REPLACE FUNCTION ST_AddBand(rast1 raster, rast2 raster, band int, index int)
RETURNS raster AS
$$
DECLARE
newraster raster;
newnodatavalue int;
newindex int := index;
newband int := band;
x int;
y int;
BEGIN
IF ST_Width(rast1) != ST_Width(rast2) OR ST_Height(rast1) != ST_Height(rast2) THEN
RAISE EXCEPTION 'ST_AddBand: Attempting to add a band with different width or height';
END IF;
IF newindex < 1 THEN
newindex := 1;
END IF;
IF newindex IS NULL OR newindex > ST_NumBands(rast1) THEN
newindex := ST_NumBands(rast1) + 1;
END IF;
IF newband < 1 THEN
newband := 1;
END IF;
IF newband IS NULL OR newband > ST_NumBands(rast2) THEN
newband := ST_NumBands(rast2);
END IF;
IF newband = 0 THEN
RETURN rast1;
END IF;
newnodatavalue := ST_BandNodataValue(rast2, newband);
newraster := ST_AddBand(rast1, newindex, ST_BandPixelType(rast2, newband), newnodatavalue, newnodatavalue);
FOR x IN 1..ST_Width(rast2) LOOP
FOR y IN 1..ST_Height(rast2) LOOP
newraster := ST_SetValue(newraster, newindex, x, y, ST_Value(rast2, newband, x, y));
END LOOP;
END LOOP;
RETURN newraster;
END;
$$
LANGUAGE 'plpgsql';
--Test
SELECT ST_NumBands(ST_AddBand(ST_TestRaster(0, 0, 1), ST_TestRaster(0, 0, 1), 1, 2))

7
raster/scripts/plpgsql/st_asraster.sql
View file

@ -1,8 +1,13 @@
----------------------------------------------------------------------
--
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_AsRaster function found in this file still many improvements before being implemented in C.
CREATE OR REPLACE FUNCTION ST_AsRaster(geom geometry, rast raster, pixeltype text, nodatavalue float8, val float8)
RETURNS raster AS
$$

40
raster/scripts/plpgsql/st_clip.sql
View file

@ -1,31 +1,35 @@
----------------------------------------------------------------------
--
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_Clip function found in this file still many improvements before being implemented in C.
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, x int, y int, width int, height int)
RETURNS raster AS
$$
DECLARE
newrast raster := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(rast), ST_UpperLeftY(rast),
ST_PixelSizeX(rast), ST_PixelSizeY(rast), ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast));
newrast raster := ST_MakeEmptyRaster(width, height, ST_UpperLeftX(rast), ST_UpperLeftY(rast),
ST_PixelSizeX(rast), ST_PixelSizeY(rast), ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast));
numband int := ST_Numbands(rast);
band int;
cx int;
cy int;
newwidth int := CASE WHEN x + width > ST_Width(rast) THEN ST_Width(rast) - x ELSE width END;
newheight int := CASE WHEN y + height > ST_Height(rast) THEN ST_Height(rast) - y ELSE height END;
band int;
cx int;
cy int;
newwidth int := CASE WHEN x + width > ST_Width(rast) THEN ST_Width(rast) - x ELSE width END;
newheight int := CASE WHEN y + height > ST_Height(rast) THEN ST_Height(rast) - y ELSE height END;
BEGIN
FOR b IN 1..numband LOOP
newrast := ST_AddBand(newrast, ST_PixelType(rast, band), ST_NodataValue(rast, band), ST_NodataValue(rast, band));
FOR cx IN 1..newwidth LOOP
FOR cy IN 1..newheight LOOP
newrast := ST_SetValue(newrast, band, cx, cy, ST_Value(rast, band, cx + x - 1, cy + y - 1));
END LOOP;
END LOOP;
END LOOP;
RETURN newrast;
FOR b IN 1..numband LOOP
newrast := ST_AddBand(newrast, ST_PixelType(rast, band), ST_NodataValue(rast, band), ST_NodataValue(rast, band));
FOR cx IN 1..newwidth LOOP
FOR cy IN 1..newheight LOOP
newrast := ST_SetValue(newrast, band, cx, cy, ST_Value(rast, band, cx + x - 1, cy + y - 1));
END LOOP;
END LOOP;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';

12
raster/scripts/plpgsql/st_deleteband.sql
View file

@ -1,8 +1,18 @@
----------------------------------------------------------------------
--
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The ST_DeleteBand function found in this file still need enhancement before being implemented in C.
-- NOTE: ST_DeleteBand(rast raster, band int) is dependent on
-- ST_AddBand(rast1 raster, rast2 raster, band int, index int)
-- to be found in the script/plpgsql folder
CREATE OR REPLACE FUNCTION ST_DeleteBand(rast raster,
band int)
RETURNS raster AS

211
raster/scripts/plpgsql/st_mapalgebra.sql
View file

@ -1,8 +1,14 @@
----------------------------------------------------------------------
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- NOTE: The one raster version of ST_MapAlgebra found in this file is ready to be being implemented in C
-- NOTE: The ST_SameAlignment function found in this files is is ready to be being implemented in C
-- NOTE: The two raster version of ST_MapAlgebra found in this file is to be replaced by the optimized version found in st_mapalgebra_optimized.sql
--------------------------------------------------------------------
-- ST_MapAlgebra - (one raster version) Return a raster which values
-- are the result of an SQL expression involving pixel
@ -34,6 +40,8 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, band integer, expression t
initexpr text;
initndvexpr text;
newpixeltype text;
skipcomputation int := 0;
newhasnodatavalue boolean := FALSE;
BEGIN
-- Check if raster is NULL
IF rast IS NULL THEN
@ -63,10 +71,12 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, band integer, expression t
IF ST_BandHasNodataValue(rast, band) THEN
newnodatavalue := ST_BandNodataValue(rast, band);
ELSE
newnodatavalue := NULL;
RAISE NOTICE 'ST_MapAlgebra: Source raster do not have a nodata value, nodata value for new raster set to the min value possible';
newnodatavalue := ST_MinPossibleVal(newrast);
END IF;
-- We set the initial value of the future band to nodata value.
-- If nodatavalue is null then the raster will be initialise to 0 and all the values will be recomputed.
-- If nodatavalue is null then the raster will be initialise to ST_MinPossibleVal
-- but all the values should be recomputed anyway.
newinitialvalue := newnodatavalue;
-- Set the new pixeltype
@ -78,44 +88,79 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, band integer, expression t
RAISE EXCEPTION 'ST_MapAlgebra: Invalid pixeltype "%". Aborting.', newpixeltype;
END IF;
initexpr := 'SELECT ' || expression;
initndvexpr := 'SELECT ' || nodatavalueexpr;
initexpr := 'SELECT ' || trim(upper(expression));
initndvexpr := 'SELECT ' || trim(upper(nodatavalueexpr));
-- First optimization: Recompute the initial value if a nodatavalueexpr is provided so we can then initialise the raster with this value.
-- This avoid having to compute nodata values one by one in the main computing loop
IF NOT nodatavalueexpr IS NULL AND NOT newnodatavalue IS NULL THEN
newexpr := replace(initndvexpr, 'rast', newnodatavalue::text);
--RAISE NOTICE '111 initexpr=%, newnodatavalue=%', initexpr,newnodatavalue;
-- Optimization: If a nodatavalueexpr is provided, recompute the initial value
-- so we can then initialise the raster with this value and skip the computation
-- of nodata values one by one in the main computing loop
IF NOT nodatavalueexpr IS NULL THEN
newexpr := replace(initndvexpr, 'RAST', newnodatavalue::text);
--RAISE NOTICE '222 newexpr=%', newexpr;
EXECUTE newexpr INTO newinitialvalue;
END IF;
-- Second optimization: If the raster is only filled with nodata value return right now a raster filled with the nodatavalueexpr
--RAISE NOTICE '333';
-- Optimization: If the raster is only filled with nodata value return right now a raster filled with the nodatavalueexpr
IF ST_BandIsNoData(rast, band) THEN
RETURN ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
END IF;
-- Third optimization: If expression resume to 'rast' and nodatavalueexpr is NULL or also equal to 'RAST', we can just return the band from the original raster
IF trim(upper(expression)) = 'RAST' AND (nodatavalueexpr IS NULL OR trim(upper(nodatavalueexpr)) = 'RAST') THEN
RETURN ST_AddBand(newrast, rast, band, 1); -- To be implemented
-- Optimization: If expression resume to 'RAST' and nodatavalueexpr is NULL or also equal to 'RAST',
-- we can just return the band from the original raster
IF initexpr = 'SELECT RAST' AND (nodatavalueexpr IS NULL OR initndvexpr = 'SELECT RAST') THEN
RETURN ST_AddBand(newrast, rast, band, 1); -- To be implemented in C
END IF;
-- Optimization: If expression resume to a constant (it does not contain RAST)
IF position('RAST' in initexpr) = 0 THEN
--RAISE NOTICE '444';
EXECUTE initexpr INTO newval;
--RAISE NOTICE '555';
skipcomputation := 1;
IF nodatavalueexpr IS NULL THEN
-- Compute the new value, set it and we will return after creating the new raster
newinitialvalue := newval;
skipcomputation := 2;
ELSEIF newval = newinitialvalue THEN
-- Return the new raster as it will be before computing pixel by pixel
skipcomputation := 2;
END IF;
END IF;
--Create the raster receiving all the computed values. Initialize it to the new initial value.
newrast := ST_AddBand(newrast, newpixeltype, newinitialvalue, newnodatavalue);
-- Optimization: If expression is NULL, or all the pixels could be set in a one step, return the initialised raster now
IF expression IS NULL OR skipcomputation = 2 THEN
RETURN newrast;
END IF;
FOR x IN 1..width LOOP
FOR y IN 1..height LOOP
r := ST_Value(rast, band, x, y);
-- We compute a value only for the withdata value pixel since the nodata value have already been set by the first optimization
-- We compute a value only for the withdata value pixel since the nodata value have already been set by the first optimization
IF NOT r IS NULL AND (newnodatavalue IS NULL OR r != newnodatavalue) THEN -- The second part of the test can be removed once ST_Value return NULL for a nodata value
newexpr := replace(initexpr, 'rast', r::text);
EXECUTE newexpr INTO newval;
IF newval IS NULL THEN
newval := newnodatavalue;
IF skipcomputation = 0 THEN
newexpr := replace(initexpr, 'RAST', r::text);
--RAISE NOTICE '666 newexpr=%', newexpr;
EXECUTE newexpr INTO newval;
IF newval IS NULL THEN
newval := newnodatavalue;
END IF;
END IF;
newrast = ST_SetValue(newrast, band, x, y, newval);
ELSE
newhasnodatavalue := TRUE;
END IF;
END LOOP;
END LOOP;
RETURN newrast;
RETURN ST_SetBandHasNodataValue(newrast, 1, newhasnodatavalue);
END;
$$
LANGUAGE 'plpgsql';
@ -141,7 +186,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, expression text)
CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, band integer, expression text, nodatavalexpr text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, $2, $3, nodatavalexpr, NULL)'
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast raster, expression text, nodatavalexpr text, pixeltype text)
@ -303,7 +348,7 @@ CREATE OR REPLACE FUNCTION ST_HasNoBand(raster, int)
--------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_BandIsNoData(raster, int)
RETURNS boolean
AS 'SELECT ST_BandHasNodataValue($1, $2) AND false'
AS 'SELECT ST_HasNoBand($1, $2) AND ST_BandHasNodataValue($1, $2) AND false'
LANGUAGE 'SQL' IMMUTABLE;
--CREATE OR REPLACE FUNCTION toto()
@ -329,15 +374,9 @@ CREATE OR REPLACE FUNCTION ST_BandIsNoData(raster, int)
-- -SECOND: Same extent as the second) raster. Default.
-- -INTERSECTION: Intersection of extent of the two rasters.
-- -UNION: Union oof extent of the two rasters.
-- rast1nodatavalrepl float8
-- rast2nodatavalrepl float8 - Values used in replacement when the raster
-- pixel involved in the expression are nodata values or
-- do not exist. When both pixels values are nodata values
-- (or one is nodata and the other does not exist), there
-- is no replacement and the expression is evaluated to nodata.
-- When nodatavalrepl is NULL the missing value is replaced
-- with the 'NULL' string in the expression and can be evaluated
-- using 'IS NULL' in the expression. e.g. "CASE WHEN rast2 IS NULL OR rast1 < rast2 THEN rast1 WHEN rast1 IS NULL OR rast2 < rast1 THEN rast2 END"
-- nodata1expr text - Expression used when rast1 pixel value is nodata or absent.
-- nodata2expr text - Expression used when rast2 pixel value is nodata or absent.
-- nodatanodataexpr text - Expression used when both pixel values are nodata values or absent.
-- Further enhancements:
-- -Optimization for UNION & FIRST. We might not have to iterate over all the new raster.
@ -353,8 +392,9 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
expression text,
pixeltype text,
extentexpr text,
rast1nodatavalrepl float8,
rast2nodatavalrepl float8)
nodata1expr text,
nodata2expr text,
nodatanodatadaexpr text)
RETURNS raster AS
$$
DECLARE
@ -412,7 +452,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
newval float;
newexpr text;
initexpr text;
BEGIN
-- We have to deal with NULL, empty, hasnoband and hasnodatavalue band rasters...
-- These are respectively tested by "IS NULL", "ST_IsEmpty()", "ST_HasNoBand()" and "ST_BandIsNodata()"
@ -714,16 +754,9 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
newpixeltype := ST_BandPixelType(rast1, band1);
END IF;
END IF;
-- Determine new notadavalue
IF (upper(extentexpr) = 'SECOND' AND NOT ST_HasNoBand(rast2, band2)) OR ST_HasNoBand(rast1, band1) THEN
newnodatavalue := ST_BandNodataValue(rast2, band2);
ELSE
newnodatavalue := ST_BandNodataValue(rast1, band1);
END IF;
-- Get the nodata value for first raster and set newnodatavalue
IF NOT ST_HasNoBand(rast1, band1) AND ST_BandHasNodatavalue(rast1, band1) THEN
-- Get the nodata value for first raster
IF NOT ST_HasNoBand(rast1, band1) AND ST_BandHasNodataValue(rast1, band1) THEN
rast1nodataval := ST_BandNodatavalue(rast1, band1);
ELSE
rast1nodataval := NULL;
@ -735,6 +768,16 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
rast2nodataval := NULL;
END IF;
-- Determine new notadavalue
IF (upper(extentexpr) = 'SECOND' AND NOT rast2nodataval IS NULL) THEN
newnodatavalue := rast2nodataval;
ELSEIF NOT rast1nodataval IS NULL THEN
newnodatavalue := rast1nodataval;
ELSE
RAISE NOTICE 'ST_MapAlgebra: Both source rasters do not have a nodata value, nodata value for new raster set to the minimum value possible';
newnodatavalue := ST_MinPossibleVal(newrast);
END IF;
-------------------------------------------------------------------
--Create the raster receiving all the computed values. Initialize it to the new nodatavalue.
newrast := ST_AddBand(ST_MakeEmptyRaster(newwidth, newheight, newulx, newuly, newpixelsizex, newpixelsizey, newskewx, newskewy, newsrid), newpixeltype, newnodatavalue, newnodatavalue);
@ -744,52 +787,49 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
-- and there is no replacement value for those missing values
-- and this raster IS involved in the expression
-- return NOW with the nodata band raster.
--IF (rast1 IS NULL OR ST_IsEmpty(rast1) OR ST_HasNoBand(rast1, band1) OR ST_BandIsNoData(rast1, band1)) AND nodatavalrepl IS NULL AND position('rast1' in expression) != 0 THEN
--IF (rast1 IS NULL OR ST_IsEmpty(rast1) OR ST_HasNoBand(rast1, band1) OR ST_BandIsNoData(rast1, band1)) AND nodatavalrepl IS NULL AND position('RAST1' in upper(expression)) != 0 THEN
-- RETURN newrast;
--END IF;
--IF (rast2 IS NULL OR ST_IsEmpty(rast2) OR ST_HasNoBand(rast2, band2) OR ST_BandIsNoData(rast2, band2)) AND nodatavalrepl IS NULL AND position('rast2' in expression) != 0 THEN
--IF (rast2 IS NULL OR ST_IsEmpty(rast2) OR ST_HasNoBand(rast2, band2) OR ST_BandIsNoData(rast2, band2)) AND nodatavalrepl IS NULL AND position('RAST2' in upper(expression)) != 0 THEN
-- RETURN newrast;
--END IF;
initexpr := 'SELECT ' || expression;
-- There is place for optimization here when doing a UNION we don't want to iterate over the empty space.
FOR x IN 1..newwidth LOOP
FOR y IN 1..newheight LOOP
r1 := ST_Value(rast1, band1, x - rast1offsetx, y - rast1offsety);
r2 := ST_Value(rast2, band2, x - rast2offsetx1, y - rast2offsety1);
-- Check if both values are outside the extent or nodata values
IF (r1 IS NULL OR r1 = rast1nodataval) AND (r2 IS NULL OR r2 = rast2nodataval) AND rast1nodatavalrepl IS NULL AND rast2nodatavalrepl IS NULL THEN
newval := newnodatavalue;
ELSE
newexpr := initexpr;
-- Check if r1 or r2 are outside the extent or are nodata value
IF r1 IS NULL OR r1 = rast1nodataval THEN
IF rast1nodatavalrepl IS NULL THEN
newexpr := replace(newexpr, 'rast1', 'NULL');
ELSE
newexpr := replace(newexpr, 'rast1', rast1nodatavalrepl::text);
END IF;
IF (r1 IS NULL OR r1 = rast1nodataval) AND (r2 IS NULL OR r2 = rast1nodataval) THEN
IF nodatanodataexpr IS NULL THEN
newval := newnodatavalue;
ELSE
newexpr := replace(newexpr, 'rast1', r1::text);
EXECUTE 'SELECT ' || nodatanodataexpr INTO newval;
END IF;
IF r2 IS NULL OR r2 = rast2nodataval THEN
IF rast2nodatavalrepl IS NULL THEN
newexpr := replace(newexpr, 'rast2', 'NULL');
ELSE
newexpr := replace(newexpr, 'rast2', rast2nodatavalrepl::text);
END IF;
ELSIF r1 IS NULL OR r1 = rast1nodataval THEN
IF nodata1expr IS NULL THEN
newval := newnodatavalue;
ELSE
newexpr := replace(newexpr, 'rast2', r2::text);
newexpr := replace('SELECT ' || upper(nodata1expr), 'RAST2', r2);
EXECUTE newexpr INTO newval;
END IF;
-- Evaluate the epression
EXECUTE newexpr INTO newval;
IF newval IS NULL THEN
ELSIF r2 IS NULL OR r2 = rast2nodataval THEN
IF nodata2expr IS NULL THEN
newval := newnodatavalue;
ELSE
newexpr := replace('SELECT ' || upper(nodata2expr), 'RAST1', r1);
EXECUTE newexpr INTO newval;
END IF;
ELSE
newexpr := replace('SELECT ' || upper(expression), 'RAST1', r1);
newexpr := replace(newexpr, 'RAST2', r2);
EXECUTE newexpr INTO newval;
END IF;
IF newval IS NULL THEN
newval := newnodatavalue;
END IF;
--RAISE NOTICE 'x=%, y=%, r1=%, r2=%, newval=%',x,y,r1,r2,newval;
newrast = ST_SetValue(newrast, 1, x, y, newval);
END LOOP;
END LOOP;
@ -798,6 +838,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
$$
LANGUAGE 'plpgsql';
--------------------------------------------------------------------
-- ST_MapAlgebra (two raster version) variants
--------------------------------------------------------------------
@ -808,10 +849,11 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
expression text,
pixeltype text,
extentexpr text,
rast1nodatavalrepl float8,
rast2nodatavalrepl float8)
nodata1expr text,
nodata2expr text,
nodatanodataexpr text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, $5, $6, $7)'
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, $5, $6, $7, $8)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant 6
@ -821,17 +863,16 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
pixeltype text,
extentexpr text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, $5, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, $5, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant 7
CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
rast2 raster,
expression text,
pixeltype text,
extentexpr text)
pixeltype text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, $5, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, $4, NULL, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
@ -840,7 +881,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
rast2 raster,
expression text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, NULL, NULL, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, 1, $2, 1, $3, NULL, NULL, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE STRICT;
-- Variant 10
@ -852,9 +893,9 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
expression text,
pixeltype text,
extentexpr text,
rastnodatavalrepl float8)
nodataexpr text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, $7, $8, $8)'
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, $7, $8, $8, $8)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant 11
@ -866,7 +907,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
pixeltype text,
extentexpr text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, $7, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, $7, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant 12
@ -877,7 +918,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
expression text,
pixeltype text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, NULL, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, $6, NULL, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant 13
@ -887,7 +928,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster,
band2 integer,
expression text)
RETURNS raster
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, NULL, NULL, NULL, NULL)'
AS 'SELECT ST_MapAlgebra($1, $2, $3, $4, $5, NULL, NULL, NULL, NULL, NULL)'
LANGUAGE 'SQL' IMMUTABLE;
--test MapAlgebra with NULL

652
raster/scripts/plpgsql/st_mapalgebra_optimized.sql Normal file
View file

@ -0,0 +1,652 @@
----------------------------------------------------------------------
--
-- $Id: st_mapalgebra_optimized.sql 6127 2010-10-25 16:06:00Z jorgearevalo $
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- Note: The functions provided in this script are in developement. Do not use.
-- Note: this script is dependent on
-- _MapAlgebraParts(r1x int, r1y int, r1w int, r1h int, r2x int, r2y int, r2w int, r2h int)
-- ST_SameAlignment(rast1ulx float8, rast1uly float8, rast1pixsizex float8, rast1pixsizey float8, rast1skewx float8, rast1skewy float8, rast2ulx float8, rast2uly float8, rast2pixsizex float8, rast2pixsizey float8, rast2skewx float8, rast2skewy float8)
-- ST_IsEmpty(raster)
-- ST_HasNoBand(raster, int)
-- ST_BandIsNoData(raster, int)
-- to be found in the script/plpgsql folder
--------------------------------------------------------------------
-- ST_MapAlgebra - (two rasters version) Return a raster which
-- values are the result of an SQL expression involving
-- pixel values from input rasters bands.
-- Arguments
-- rast1 raster - First raster referred by rast1 in the expression.
-- band1 integer - Band number of the first raster. Default to 1.
-- rast2 raster - Second raster referred by rast2 in the expression.
-- band2 integer - Band number of the second raster. Default to 1.
-- expression text - SQL expression. Ex.: "rast1 + 2 * rast2"
-- pixeltype text - Pixeltype assigned to the resulting raster. Expression
-- results are truncated to this type. Default to the
-- pixeltype of the first raster.
-- extentexpr text - Raster extent of the result. Can be:
-- -FIRST: Same extent as the first raster. Default.
-- -SECOND: Same extent as the second) raster. Default.
-- -INTERSECTION: Intersection of extent of the two rasters.
-- -UNION: Union oof extent of the two rasters.
-- nodata1expr text - Expression used when only rast1 pixel value are nodata or absent, i.e. rast2 pixel value is with data.
-- nodata2expr text - Expression used when only rast2 pixel value are nodata or absent, i.e. rast1 pixel value is with data.
-- nodatanodataexpr text - Expression used when both pixel values are nodata values or absent.
-- Further enhancements:
-- -Optimization for UNION & FIRST. We might not have to iterate over all the new raster.
-- -Make the function being able to look for neighbour pixels. Ex. rast1[-1,1] or rast2[-3,2]
-- -Make the function to work with neighbour tiles pixels.
-- -Resample the second raster when necessary (Require ST_Resample)
-- -More test with rotated images
--------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_MapAlgebra2(rast1 raster,
band1 integer,
rast2 raster,
band2 integer,
expression text,
pixeltype text,
extentexpr text,
nodata1expr text,
nodata2expr text,
nodatanodataexpr text)
RETURNS raster AS
$$
DECLARE
x integer;
y integer;
r1 float;
r2 float;
rast1ulx float8;
rast1uly float8;
rast1width int;
rast1height int;
rast1pixsizex float8;
rast1pixsizey float8;
rast1skewx float8;
rast1skewy float8;
rast1nodataval float8;
rast1srid int;
rast1offsetx int;
rast1offsety int;
rast2ulx float8;
rast2uly float8;
rast2width int;
rast2height int;
rast2pixsizex float8;
rast2pixsizey float8;
rast2skewx float8;
rast2skewy float8;
rast2nodataval float8;
rast2srid int;
r1x int;
r1y int;
r1w int;
r1h int;
r2x int;
r2y int;
r2w int;
r2h int;
newrx int;
newry int;
newrast raster;
tmprast raster;
newsrid int;
newpixelsizex float8;
newpixelsizey float8;
newskewx float8;
newskewy float8;
newnodatavalue float8;
newpixeltype text;
newulx float8;
newuly float8;
newwidth int;
newheight int;
newoffsetx1 int;
newoffsety1 int;
newoffsetx2 int;
newoffsety2 int;
newval float;
newexpr text;
upnodatanodataexpr text;
upnodata1expr text;
upnodata2expr text;
upexpression text;
nodatanodataval float;
skipcomputation int;
zones int[];
z11x int;
z11y int;
z11w int;
z11h int;
z12x int;
z12y int;
z12w int;
z12h int;
z13x int;
z13y int;
z13w int;
z13h int;
z14x int;
z14y int;
z14w int;
z14h int;
z21x int;
z21y int;
z21w int;
z21h int;
z22x int;
z22y int;
z22w int;
z22h int;
z23x int;
z23y int;
z23w int;
z23h int;
z24x int;
z24y int;
z24w int;
z24h int;
zcx int;
zcy int;
zcw int;
zch int;
BEGIN
-- We have to deal with NULL, empty, hasnoband and hasnodatavalue band rasters...
-- These are respectively tested by "IS NULL", "ST_IsEmpty()", "ST_HasNoBand()" and "ST_BandIsNodata()"
-- If both raster are null, we return NULL. ok
-- If both raster do not have extent (are empty), we return an empty raster. ok
-- If both raster do not have the specified band,
-- we return a no band raster with the correct extent (according to the extent expression). ok
-- If both raster bands are nodatavalue and there is no replacement value, we return a nodata value band. ok
-- If only one raster is null or empty or has no band or hasnodata band we treat it as a nodata band raster.
-- If there is a replacement value we replace the missing raster values with this replacement value. ok
-- If there is no replacement value, we return a nodata value band. ok
-- What to do when only one raster is NULL or empty
-- If the extent expression is FIRST and the first raster is null we return NULL. ok
-- If the extent expression is FIRST and the first raster do not have extent (is empty), we return an empty raster. ok
-- If the extent expression is SECOND and the second raster is null we return NULL. ok
-- If the extent expression is SECOND and the second raster do not have extent (is empty), we return an empty raster. ok
-- If the extent expression is INTERSECTION and one raster is null or do not have extent (is empty), we return an empty raster. ok
-- If the extent expression is UNION and one raster is null or do not have extent (is empty),
-- we return a raster having the extent and the band characteristics of the other raster. ok
-- What to do when only one raster do not have the required band.
-- If the extent expression is FIRST and the first raster do not have the specified band,
-- we return a no band raster with the correct extent (according to the extent expression). ok
-- If the extent expression is SECOND and the second raster do not have the specified band,
-- we return a no band raster with the correct extent (according to the extent expression). ok
-- If the extent expression is INTERSECTION and one raster do not have the specified band,
-- we treat it as a nodata raster band. ok
-- If the extent expression is UNION and one raster do not have the specified band,
-- we treat it as a nodata raster band. ok
-- In all those cases, we make a warning.
-- Check if both rasters are NULL
RAISE NOTICE 'ST_MapAlgebra2 000';
IF rast1 IS NULL AND rast2 IS NULL THEN
RAISE NOTICE 'ST_MapAlgebra: Both raster are NULL. Returning NULL';
RETURN NULL;
END IF;
-- Check if both rasters are empty (width or height = 0)
IF ST_IsEmpty(rast1) AND ST_IsEmpty(rast2) THEN
RAISE NOTICE 'ST_MapAlgebra: Both raster are empty. Returning an empty raster';
RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, -1);
END IF;
rast1ulx := ST_UpperLeftX(rast1);
rast1uly := ST_UpperLeftY(rast1);
rast1width := ST_Width(rast1);
rast1height := ST_Height(rast1);
rast1pixsizex := ST_PixelSizeX(rast1);
rast1pixsizey := ST_PixelSizeY(rast1);
rast1skewx := ST_SkewX(rast1);
rast1skewy := ST_SkewY(rast1);
rast1srid := ST_SRID(rast1);
rast2ulx := ST_UpperLeftX(rast2);
rast2uly := ST_UpperLeftY(rast2);
rast2width := ST_Width(rast2);
rast2height := ST_Height(rast2);
rast2pixsizex := ST_PixelSizeX(rast2);
rast2pixsizey := ST_PixelSizeY(rast2);
rast2skewx := ST_SkewX(rast2);
rast2skewy := ST_SkewY(rast2);
rast2srid := ST_SRID(rast2);
-- Check if the first raster is NULL or empty
IF rast1 IS NULL OR ST_IsEmpty(rast1) THEN
rast1ulx := rast2ulx;
rast1uly := rast2uly;
rast1width := rast2width;
rast1height := rast2height;
rast1pixsizex := rast2pixsizex;
rast1pixsizey := rast2pixsizey;
rast1skewx := rast2skewx;
rast1skewy := rast2skewy;
rast1srid := rast2srid;
END IF;
-- Check if the second raster is NULL or empty
IF rast2 IS NULL OR ST_IsEmpty(rast2) THEN
rast2ulx := rast1ulx;
rast2uly := rast1uly;
rast2width := rast1width;
rast2height := rast1height;
rast2pixsizex := rast1pixsizex;
rast2pixsizey := rast1pixsizey;
rast2skewx := rast1skewx;
rast2skewy := rast1skewy;
rast2srid := rast1srid;
END IF;
-- Check for SRID
IF rast1srid != rast2srid THEN
RAISE EXCEPTION 'ST_MapAlgebra: Provided raster with different SRID. Aborting';
END IF;
newsrid := rast1srid;
-- Check for alignment. (Rotation problem here)
IF NOT ST_SameAlignment(rast1ulx, rast1uly, rast1pixsizex, rast1pixsizey, rast1skewx, rast1skewy, rast2ulx, rast2uly, rast2pixsizex, rast2pixsizey, rast2skewx, rast2skewy) THEN
-- For now print an error message, but a more robust implementation should resample the second raster to the alignment of the first raster.
RAISE EXCEPTION 'ST_MapAlgebra: Provided raster do not have the same alignment. Aborting';
END IF;
-- Set new pixel size and skew. We set it to the rast1 pixelsize and skew
-- since both rasters are aligned and thus have the same pixelsize and skew
newpixelsizex := rast1pixsizex;
newpixelsizey := rast1pixsizey;
newskewx := rast1skewx;
newskewy := rast1skewx;
--r1x & r2x are the offset of each rasters relatively to global extent
r1x := 0;
r2x := -st_world2rastercoordx(rast2, rast1ulx, rast1uly) + 1;
IF r2x < 0 THEN
r1x := -r2x;
r2x := 0;
END IF;
r1y := 0;
r2y := -st_world2rastercoordy(rast2, rast1ulx, rast1uly) + 1;
IF r2y < 0 THEN
r1y := -r2y;
r2y := 0;
END IF;
r1w := rast1width;
r1h := rast1height;
r2w := rast2width;
r2h := rast2height;
zones := _MapAlgebraParts(r1x + 1, r1y + 1, r1w, r1h, r2x + 1, r2y + 1, r2w, r2h);
z11x := zones[1];
z11y := zones[2];
z11w := zones[3];
z11h := zones[4];
z12x := zones[5];
z12y := zones[6];
z12w := zones[7];
z12h := zones[8];
z13x := zones[9];
z13y := zones[10];
z13w := zones[11];
z13h := zones[12];
z14x := zones[13];
z14y := zones[14];
z14w := zones[15];
z14h := zones[16];
z21x := zones[17];
z21y := zones[18];
z21w := zones[19];
z21h := zones[20];
z22x := zones[21];
z22y := zones[22];
z22w := zones[23];
z22h := zones[24];
z23x := zones[25];
z23y := zones[26];
z23w := zones[27];
z23h := zones[28];
z24x := zones[29];
z24y := zones[30];
z24w := zones[31];
z24h := zones[32];
zcx := zones[33];
zcy := zones[34];
zcw := zones[35];
zch := zones[36];
-- Compute x and y relative index of master and slave according to the extent expression (FIRST, SECOND, INTERSECTION or UNION)
IF extentexpr IS NULL OR upper(extentexpr) = 'FIRST' THEN
-- Check if rast1 is NULL
IF rast1 IS NULL THEN
RAISE NOTICE 'ST_MapAlgebra: FIRST raster is NULL. Returning NULL';
RETURN NULL;
END IF;
-- Check if rast1 is empty
IF ST_IsEmpty(rast1) THEN
RAISE NOTICE 'ST_MapAlgebra: FIRST raster is empty. Returning an empty raster';
RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid);
END IF;
-- Check if rast1 has the required band
IF ST_HasNoBand(rast1, band1) THEN
RAISE NOTICE 'ST_MapAlgebra: FIRST raster has no band. Returning a raster without band';
RETURN ST_MakeEmptyRaster(rast1width, rast1height, rast1ulx, rast1uly, rast1pixsizex, rast1pixsizey, rast1skewx, rast1skewy, rast1srid);
END IF;
newulx := rast1ulx;
newuly := rast1uly;
newwidth := rast1width;
newheight := rast1height;
newrx := r1x;
newry := r1y;
z21w := 0;
z22w := 0;
z23w := 0;
z24w := 0;
ELSIF upper(extentexpr) = 'SECOND' THEN
-- Check if rast2 is NULL
IF rast2 IS NULL THEN
RAISE NOTICE 'ST_MapAlgebra: SECOND raster is NULL. Returning NULL';
RETURN NULL;
END IF;
-- Check if rast2 is empty
IF ST_IsEmpty(rast2) THEN
RAISE NOTICE 'ST_MapAlgebra: SECOND raster is empty. Returning an empty raster';
RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid);
END IF;
-- Check if rast2 has the required band
IF ST_HasNoBand(rast2, band2) THEN
RAISE NOTICE 'ST_MapAlgebra: SECOND raster has no band. Returning an empty raster';
RETURN ST_MakeEmptyRaster(rast2width, rast2height, rast2ulx, rast2uly, rast2pixsizex, rast2pixsizey, rast2skewx, rast2skewy, rast2srid);
END IF;
newulx := rast2ulx;
newuly := rast2uly;
newwidth := rast2width;
newheight := rast2height;
newrx := r2x;
newry := r2y;
z11w := 0;
z12w := 0;
z13w := 0;
z14w := 0;
ELSIF upper(extentexpr) = 'INTERSECTION' THEN
-- Check if the intersection is empty.
IF zcw = 0 OR zch = 0 OR
rast1 IS NULL OR ST_IsEmpty(rast1) OR
rast2 IS NULL OR ST_IsEmpty(rast2) THEN
RAISE NOTICE 'ST_MapAlgebra: INTERSECTION of provided rasters is empty. Returning an empty raster';
RETURN ST_MakeEmptyRaster(0, 0, 0, 0, 0, 0, 0, 0, newsrid);
END IF;
-- Compute the new ulx and uly
newulx := st_raster2worldcoordx(rast1, zcx - r1x + 1, zcy - r1y + 1);
newuly := st_raster2worldcoordy(rast1, zcx - r1x + 1, zcy - r1y + 1);
newwidth := zcw;
newheight := zch;
newrx := zcx;
newry := zcy;
z11w := 0;
z12w := 0;
z13w := 0;
z14w := 0;
z21w := 0;
z22w := 0;
z23w := 0;
z24w := 0;
ELSIF upper(extentexpr) = 'UNION' THEN
IF rast1 IS NULL OR ST_IsEmpty(rast1) THEN
newulx := rast2ulx;
newuly := rast2uly;
newwidth := rast2width;
newheight := rast2height;
newrx := r2x;
newry := r2y;
z21w := 0;
z22w := 0;
z23w := 0;
z24w := 0;
ELSIF rast2 IS NULL OR ST_IsEmpty(rast2) THEN
newulx := rast1ulx;
newuly := rast1uly;
newwidth := rast1width;
newheight := rast1height;
newrx := r1x;
newry := r1y;
z11w := 0;
z12w := 0;
z13w := 0;
z14w := 0;
ELSE
newrx := 0;
newry := 0;
newulx := st_raster2worldcoordx(rast1, r1x + 1, r1y + 1);
newuly := st_raster2worldcoordy(rast1, r1x + 1, r1y + 1);
newwidth := max(r1x + r1w, r2x + r2w);
newheight := max(r1y + r1h, r2y + r2h);
END IF;
ELSE
RAISE EXCEPTION 'ST_MapAlgebra: Unhandled extent expression "%". Only MASTER, INTERSECTION and UNION are accepted. Aborting.', upper(extentexpr);
END IF;
-- Check if both rasters do not have the specified band.
IF ST_HasNoband(rast1, band1) AND ST_HasNoband(rast2, band2) THEN
RAISE NOTICE 'ST_MapAlgebra: Both raster do not have the specified band. Returning a no band raster with the correct extent';
RETURN ST_MakeEmptyRaster(newwidth, newheight, newulx, newuly, newpixelsizex, newpixelsizey, newskewx, newskewy, newsrid);
END IF;
-- Check newpixeltype
newpixeltype := pixeltype;
IF newpixeltype NOTNULL AND 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_MapAlgebra: Invalid pixeltype "%". Aborting.', newpixeltype;
END IF;
-- If no newpixeltype was provided, get it from the provided rasters.
IF newpixeltype IS NULL THEN
IF (upper(extentexpr) = 'SECOND' AND NOT ST_HasNoBand(rast2, band2)) OR ST_HasNoBand(rast1, band1) THEN
newpixeltype := ST_BandPixelType(rast2, band2);
ELSE
newpixeltype := ST_BandPixelType(rast1, band1);
END IF;
END IF;
-- Get the nodata value for first raster
IF NOT ST_HasNoBand(rast1, band1) AND ST_BandHasNodataValue(rast1, band1) THEN
rast1nodataval := ST_BandNodatavalue(rast1, band1);
ELSE
rast1nodataval := NULL;
END IF;
-- Get the nodata value for second raster
IF NOT ST_HasNoBand(rast2, band2) AND ST_BandHasNodatavalue(rast2, band2) THEN
rast2nodataval := ST_BandNodatavalue(rast2, band2);
ELSE
rast2nodataval := NULL;
END IF;
-- Determine new notadavalue
IF (upper(extentexpr) = 'SECOND' AND NOT rast2nodataval IS NULL) THEN
newnodatavalue := rast2nodataval;
ELSEIF NOT rast1nodataval IS NULL THEN
newnodatavalue := rast1nodataval;
ELSE
RAISE NOTICE 'ST_MapAlgebra: Both source rasters do not have a nodata value, nodata value for new raster set to the minimum value possible';
newnodatavalue := ST_MinPossibleVal(newrast);
END IF;
upnodatanodataexpr := upper(nodatanodataexpr);
upnodata1expr := upper(nodata1expr);
upnodata2expr := upper(nodata2expr);
upexpression := upper(expression);
-- Initialise newrast with nodata-nodata values. Then we don't have anymore to set values for nodata-nodata pixels.
IF upnodatanodataexpr IS NULL THEN
nodatanodataval := newnodatavalue;
ELSE
EXECUTE 'SELECT ' || upnodatanodataexpr INTO nodatanodataval;
IF nodatanodataval IS NULL THEN
nodatanodataval := newnodatavalue;
ELSE
newnodatavalue := nodatanodataval;
END IF;
END IF;
-------------------------------------------------------------------
--Create the raster receiving all the computed values. Initialize it to the new nodatavalue.
newrast := ST_AddBand(ST_MakeEmptyRaster(newwidth, newheight, newulx, newuly, newpixelsizex, newpixelsizey, newskewx, newskewy, newsrid), newpixeltype, newnodatavalue, newnodatavalue);
-------------------------------------------------------------------
RAISE NOTICE 'ST_MapAlgebra2 111 z11x=%, z11y=%, z11w=%, z11h=%', z11x, z11y, z11w, z11h;
-- First zone with only rast1 (z11)
IF z11w > 0 AND z11h > 0 AND NOT ST_BandIsNodata(rast1, band1) AND NOT nodata2expr IS NULL THEN
IF upnodata2expr = 'RAST' THEN
-- IF rast1nodataval != nodatanodataval THEN
RAISE NOTICE 'ST_MapAlgebra2 222';
-- newrast := ST_SetValues(newrast, 1, z11x, z11y, z11w, z11h, nodatanodataval);
-- END IF;
RAISE NOTICE 'ST_MapAlgebra2 333 z11x=%, z11y=%, z11w=%, z11h=%', z11x, z11y, z11w, z11h;
newrast := ST_SetValues(newrast, 1, z11x, z11y, z11w, z11h, rast1, band1, NULL, TRUE);
ELSE
RAISE NOTICE 'ST_MapAlgebra2 444';
tmprast := ST_Clip(rast1, z11x - r1x + 1, z11y - r1y + 1, z11w, z11h);
newrast := ST_Mapalgebra2(newrast, 1, tmprast, 1, replace(upnodata2expr, 'RAST', 'RAST2'), NULL, 'FIRST', NULL, 'RAST', NULL);
END IF;
END IF;
RAISE NOTICE 'ST_MapAlgebra2 555';
-- Common zone (zc)
skipcomputation = 0;
IF zcw > 0 AND zch > 0 AND (NOT ST_BandIsNodata(rast1, band1) OR NOT ST_BandIsNodata(rast2, band2)) THEN
RAISE NOTICE 'ST_MapAlgebra2 666';
-- Initialize the zone with nodatavalue. We will not further compute nodata nodata pixels
-- newrast := ST_SetValues(newrast, 1, zcx + 1, zcy + 1, zcw, zch, newnodatavalue);
-- If rast1 is only nodata values, apply nodata1expr
IF ST_BandIsNodata(rast1, band1) THEN
IF nodata1expr IS NULL THEN
-- Do nothing
skipcomputation = 0;
ELSEIF upnodata1expr = 'RAST' THEN
-- Copy rast2 into newrast
newrast := ST_SetValues(newrast, 1, zcx, zcy, zcw, zch, , rast2, band2, NULL, 'KEEP');
ELSE
-- If nodata1expr resume to a constant
IF position('RAST' in upnodata1expr) = 0 THEN
EXECUTE 'SELECT ' || upnodata1expr INTO newval;
IF newval IS NULL OR newval = newnodatavalue THEN
-- The constant is equal to nodata. We have nothing to compute since newrast was already initialized to nodata
skipcomputation := 2;
ELSEIF newnodatavalue IS NULL THEN
-- We can globally initialize to the constant only if there was no newnodatavalue.
newrast := ST_SetValues(newrast, 1, zcx, zcy, zcw, zch, newval);
skipcomputation := 2;
ELSE
-- We will assign the constant pixel by pixel.
skipcomputation := 1;
END IF;
END IF;
IF skipcomputation < 2 THEN
FOR x IN 1..zcw LOOP
FOR y IN 1..zch LOOP
r2 := ST_Value(rast2, band2, x + r2x, y + r2y);
IF (r2 IS NULL OR r2 = rast2nodataval) THEN
-- Do nothing since the raster have already been all set to this value
ELSE
IF skipcomputation < 1 THEN
newexpr := replace('SELECT ' || upnodata1expr, 'RAST', r2);
EXECUTE newexpr INTO newval;
IF newval IS NULL THEN
newval := newnodatavalue;
END IF;
END IF;
newrast = ST_SetValue(newrast, 1, x + zcx, y + zcy, newval);
END IF;
END LOOP;
END LOOP;
END IF;
END IF;
ELSEIF ST_BandIsNodata(rast2, band2) THEN
ELSE
END IF;
END IF;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
CREATE OR REPLACE FUNCTION ST_TestRaster(ulx float8, uly float8, val float8)
RETURNS raster AS
$$
DECLARE
BEGIN
RETURN ST_AddBand(ST_MakeEmptyRaster(5, 5, ulx, uly, 1, -1, 0, 0, -1), '32BF', val, -1);
END;
$$
LANGUAGE 'plpgsql';
SELECT asbinary((gv).geom), (gv).val
FROM st_pixelaspolygons(ST_TestRaster(-10, 2, 1)) gv;
SELECT asbinary(_MapAlgebraAllPartsGeom(0, 0, 2, 1, 1, 0, 2, 1))
SELECT asbinary(pix.geom) as geom, pix.val
FROM st_pixelaspolygons(ST_MapAlgebra2(ST_TestRaster(0, 1, 1), 1, ST_TestRaster(1, 0, 1), 1, '(rast1 + rast2) / 2', NULL, 'union', '2*rast', 'rast', NULL), 1) as pix

5
raster/scripts/plpgsql/st_pixelaspolygons.sql
View file

@ -1,6 +1,9 @@
-- $Id$
----------------------------------------------------------------------
--
-- $Id$
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-----------------------------------------------------------------------

18
raster/scripts/plpgsql/st_reclass.sql
View file

@ -1,8 +1,17 @@
----------------------------------------------------------------------
--
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- Note: this script is dependent on
-- ST_DeleteBand(rast raster, band int)
-- ST_AddBand(newrast, rast, b, NULL)
-- ST_MapAlgebra(rast raster, band integer, expression text, nodatavalueexpr text, pixeltype text)
-- to be found in the script/plpgsql folder
CREATE OR REPLACE FUNCTION ST_Reclass(rast raster,
band int,
reclassexpr text)
@ -47,17 +56,18 @@ CREATE OR REPLACE FUNCTION ST_Reclass(rast raster,
END IF;
-- Build the maexpr expression
IF fromstr[2] IS NULL OR fromstr[2] = '' THEN
maexpr := maexpr || ' WHEN ' || fromstr[1] || ' = rast1 THEN ' || reclassstr[2] || ' ';
maexpr := maexpr || ' WHEN ' || fromstr[1] || ' = rast THEN ' || reclassstr[2] || ' ';
ELSE
maexpr := maexpr || ' WHEN ' || fromstr[1] || ' <= rast1 AND rast1 < ' || fromstr[2] || ' THEN ' || reclassstr[2] || ' ';
maexpr := maexpr || ' WHEN ' || fromstr[1] || ' <= rast AND rast < ' || fromstr[2] || ' THEN ' || reclassstr[2] || ' ';
END IF;
END LOOP;
maexpr := maexpr || 'ELSE rast1 END';
maexpr := maexpr || 'ELSE rast END';
newrast := ST_AddBand(rast, ST_MapAlgebra(rast, band, maexpr), 1, band);
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
SELECT ST_Value(ST_TestRaster(1, 1, 4),1,1)
SELECT ST_Value(ST_Reclass(ST_TestRaster(1, 1, 4), 1, '1:2|2:2|3-5:10'), 1, 1);

419
raster/scripts/plpgsql/st_setvalues.sql Normal file
View file

@ -0,0 +1,419 @@
----------------------------------------------------------------------
--
-- $Id: st_setvalues.sql 6127 2010-10-25 16:06:00Z jorgearevalo $
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
--NOTE: Both ST_SetValues functions found in this file are ready to be being implemented in C
CREATE OR REPLACE FUNCTION max(a int, b int)
RETURNS int
AS 'SELECT CASE WHEN $1 < $2 THEN $2 ELSE $1 END'
LANGUAGE 'SQL' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION min(a int, b int)
RETURNS int
AS 'SELECT CASE WHEN $1 < $2 THEN $1 ELSE $2 END'
LANGUAGE 'SQL' IMMUTABLE STRICT;
CREATE OR REPLACE FUNCTION _Inter(x1 int, w1 int, x2 int, w2 int)
RETURNS int[] AS
$$
DECLARE
nx int;
nw int;
BEGIN
nx := x1 + min(max(0, x2 - x1), w1);
nw := max(min(x1 + w1, x2 + w2), x1) - nx;
RETURN ARRAY[nx, nw];
END;
$$
LANGUAGE 'plpgsql';
--Test
SELECT _Inter(1, 3, -2, 2);
SELECT _Inter(1, 3, 0, 2);
SELECT _Inter(1, 3, 0, 5);
SELECT _Inter(1, 3, 2, 1);
SELECT _Inter(1, 3, 2, 2);
SELECT _Inter(1, 3, 4, 2);
SELECT _Inter(-2, 2, 1, 3);
SELECT _Inter(0, 2, 1, 3);
SELECT _Inter(0, 5, 1, 3);
SELECT _Inter(2, 1, 1, 3);
SELECT _Inter(2, 2, 1, 3);
SELECT _Inter(4, 2, 1, 3);
SELECT _Inter(1, 3, 2, 0);
SELECT _Inter(4, 5, 2, 3);
--------------------------------------------------------------------
-- ST_SetValues - Set a range of raster pixels to a value.
--
-- Arguments
--
-- rast raster - Raster to be edited.
-- band integer - Band number of the raster to be edited. Default to 1.
-- x, y - Raster coordinates of the upper left corner of the range
-- of pixel to be edited.
-- width, height - Width and height of the range of pixel to be edited.
-- val - Value to set the range. If NULL, pixels are set to nodata.
-- keepdestnodata - Flag indicating not to change pixels set to nodata value.
-- Default to FALSE.
--
-- When x, y, width or height are out of the raster range, only the part
-- of the range intersecting with the raster is set.
--------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_SetValues(rast raster, band int, x int, y int, width int, height int, val float8, keepdestnodata boolean)
RETURNS raster AS
$$
DECLARE
newraster raster := rast;
cx int;
cy int;
newwidth int := width;
newheight int := height;
newband int := band;
oldx int := x;
oldy int := y;
newx int;
newy int;
newval float8 := val;
newkeepdestnodata boolean := keepdestnodata;
rastnodataval int := 0;
BEGIN
IF rast IS NULL THEN
RAISE NOTICE 'ST_SetValues: No raster provided. Returns NULL';
RETURN NULL;
END IF;
IF ST_IsEmpty(rast) OR ST_HasNoBand(rast, band) THEN
RAISE NOTICE 'ST_SetValues: Empty or no band raster provided. Returns rast';
RETURN rast;
END IF;
IF newband IS NULL THEN
newband := 1;
END IF;
IF newband < 1 THEN
RAISE NOTICE 'ST_SetValues: band out of range. Returns rast';
RETURN rast;
END IF;
IF width IS NULL OR width < 1 OR height IS NULL OR height < 1 THEN
RAISE NOTICE 'ST_SetValues: invalid width or height. Returns rast';
RETURN rast;
END IF;
IF x IS NULL THEN
oldx := 1;
END IF;
newx := 1 + min(max(0, oldx - 1), ST_Width(rast));
newwidth := max(min(1 + ST_Width(rast), oldx + newwidth), 1) - newx;
IF y IS NULL THEN
oldy := 1;
END IF;
newy := 1 + min(max(0, oldy - 1), ST_Height(rast));
newheight := max(min(1 + ST_Height(rast), oldy + newheight), 1) - newy;
IF newwidth < 1 OR newheight < 1 THEN
RETURN rast;
END IF;
IF newkeepdestnodata IS NULL THEN
newkeepdestnodata := FALSE;
END IF;
IF newkeepdestnodata THEN
IF ST_BandHasNodataValue(rast, newband) THEN
rastnodataval := ST_BandNoDataValue(rast, newband);
ELSE
RAISE NOTICE 'ST_SetValues: keepdestnodata was set to TRUE but rast1 does not have a nodata value. keepdestnodata reset to FALSE';
newkeepdestnodata := FALSE;
END IF;
IF ST_BandIsNodata(rast, band) THEN
RETURN rast;
END IF;
END IF;
IF newval IS NULL THEN
newval := ST_BandNoDataValue(rast, band);
END IF;
IF newval IS NULL THEN
RAISE NOTICE 'ST_SetValues: val is NULL and no nodata exist for rast. Returns rast';
RETURN rast;
END IF;
IF newkeepdestnodata THEN
FOR cx IN newx..newx + newwidth - 1 LOOP
FOR cy IN newy..newy + newheight - 1 LOOP
IF ST_Value(newraster, newband, cx, cy) != rastnodataval THEN
newraster := ST_SetValue(newraster, newband, cx, cy, newval);
END IF;
END LOOP;
END LOOP;
ELSE
FOR cx IN newx..newx + newwidth - 1 LOOP
FOR cy IN newy..newy + newheight - 1 LOOP
newraster := ST_SetValue(newraster, newband, cx, cy, newval);
END LOOP;
END LOOP;
END IF;
RETURN newraster;
END;
$$
LANGUAGE 'plpgsql';
-- Variant with band = 1
CREATE OR REPLACE FUNCTION ST_SetValues(rast raster, x int, y int, width int, height int, val float8, keepdestnodata boolean)
RETURNS raster
AS 'SELECT ST_SetValues($1, 1, $2, $3, $4, $5, $6, $7)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant with band = 1 & keepdestnodata = FALSE
CREATE OR REPLACE FUNCTION ST_SetValues(rast raster, x int, y int, width int, height int, val float8)
RETURNS raster
AS 'SELECT ST_SetValues($1, 1, $2, $3, $4, $5, $6, FALSE)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant with keepdestnodata = FALSE
CREATE OR REPLACE FUNCTION ST_SetValues(rast raster, band int, x int, y int, width int, height int, val float8)
RETURNS raster
AS 'SELECT ST_SetValues($1, $2, $3, $4, $5, $6, $7, FALSE)'
LANGUAGE 'SQL' IMMUTABLE;
-- Test
SELECT ST_SetValues(ST_TestRaster(0, 0, 1), 2, 2, 1, 1, 0)
SELECT (pix).geom, (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_TestRaster(0, 0, 1)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, 1), 2, 1, 1, 1, 0)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, 1), 2, 1, 1, 10, 0)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, 1), 1, 1, 4, 4, 0)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, 1), 0, 3, 4, 4, 0)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, -1), 2, 2, 2, 2, 0, TRUE)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_SetBandHasNoDataValue(ST_TestRaster(0, 0, -1), FALSE), 2, 2, 2, 2, 0, TRUE)) as pix) foo
--------------------------------------------------------------------
-- ST_SetValues - Set a range of raster pixels to values copied from
-- the corresponding pixels in another raster.
-- Arguments
--
-- rast1 raster - Raster to be edited.
-- band1 integer - Band number of the raster to be edited. Default to 1.
-- x, y - Raster coordinates of the upper left corner of the
-- range of pixel to be edited.
-- width, height - Width and height of the range of pixel to be edited.
-- rast2 - Raster values are copied from.
-- band2 - Band number of the raster values are copied from.
-- keepdestnodata - Flag indicating not to change pixels (in the edited
-- raster) set to nodata value. Default to FALSE.
-- keepsourcetnodata - Flag indicating not to copy pixels (from the source
-- raster) set to nodata value. Default to FALSE.
--
-- When x, y, width or height are out of the raster range, only the part
-- of the range intersecting with the raster is set.
--------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_SetValues(rast1 raster, band1 int, x int, y int, width int, height int, rast2 raster, band2 int, keepdestnodata boolean, keepsourcenodata boolean)
RETURNS raster AS
$$
DECLARE
newraster raster := rast1;
newwidth int := width;
newheight int := height;
oldx int := x;
oldy int := y;
newx int;
newy int;
newband1 int := band1;
newband2 int := band2;
newkeepdestnodata boolean := keepdestnodata;
newkeepsourcenodata boolean := keepsourcenodata;
r2val float;
rast1nodataval float;
rast2nodataval float;
x2 int;
y2 int;
cx int;
cy int;
BEGIN
IF rast1 IS NULL THEN
RAISE NOTICE 'ST_SetValues: No raster provided. Return NULL';
RETURN NULL;
END IF;
IF ST_IsEmpty(rast1) OR ST_HasNoBand(rast1, band1) THEN
RAISE NOTICE 'ST_SetValues: Empty or no band destination raster provided. Returns rast1';
RETURN rast1;
END IF;
IF ST_IsEmpty(rast2) OR ST_HasNoBand(rast2, band2) OR rast2 IS NULL THEN
RAISE NOTICE 'ST_SetValues: Empty or no band source raster provided. Returns rast1';
RETURN rast1;
END IF;
IF newband1 IS NULL THEN
newband1 := 1;
END IF;
IF newband1 < 1 THEN
RAISE NOTICE 'ST_SetValues: band1 out of range. Returns rast';
RETURN rast1;
END IF;
IF newband2 IS NULL THEN
newband2 := 1;
END IF;
IF newband2 < 1 THEN
RAISE NOTICE 'ST_SetValues: band2 out of range. Returns rast';
RETURN rast1;
END IF;
IF x IS NULL THEN
oldx := 1;
END IF;
newx := 1 + min(max(0, oldx - 1), ST_Width(rast1));
newwidth := max(min(1 + ST_Width(rast1), oldx + newwidth), 1) - newx;
oldx := newx;
IF y IS NULL THEN
oldy := 1;
END IF;
--RAISE NOTICE 'aaa oldy=%, newheight=%', oldy, newheight;
newy := 1 + min(max(0, oldy - 1), ST_Height(rast1));
newheight := max(min(1 + ST_Height(rast1), oldy + newheight), 1) - newy;
oldy := newy;
--RAISE NOTICE 'bbb newx=%, newy=%', newx, newy;
--RAISE NOTICE 'ccc newwidth=%, newheight=%', newwidth, newheight;
IF newwidth < 1 OR newheight < 1 THEN
RETURN rast1;
END IF;
x2 := ST_World2RasterCoordX(rast1, ST_Raster2WorldCoordX(rast2, 1, 1), ST_Raster2WorldCoordY(rast2, 1, 1));
y2 := ST_World2RasterCoordY(rast1, ST_Raster2WorldCoordY(rast2, 1, 1), ST_Raster2WorldCoordY(rast2, 1, 1));
--RAISE NOTICE '111 x2=%, y2=%', x2, y2;
newx := x2 + min(max(0, oldx - x2), ST_Width(rast2));
newwidth := max(min(x2 + ST_Width(rast2), oldx + newwidth), x2) - newx;
newy := y2 + min(max(0, oldy - y2), ST_Height(rast2));
newheight := max(min(y2 + ST_Height(rast2), oldy + newheight), y2) - newy;
IF newwidth < 1 OR newheight < 1 THEN
RETURN rast1;
END IF;
--RAISE NOTICE '222 newx=%, newy=%', newx, newy;
--RAISE NOTICE '333 newwidth=%, newheight=%', newwidth, newheight;
IF newkeepdestnodata IS NULL THEN
newkeepdestnodata := FALSE;
END IF;
IF newkeepdestnodata THEN
IF ST_BandHasNodataValue(rast1, newband1) THEN
rast1nodataval := ST_BandNoDataValue(rast1, newband1);
ELSE
RAISE NOTICE 'ST_SetValues: keepdestnodata was set to TRUE but rast1 does not have a nodata value. keepdestnodata reset to FALSE';
newkeepdestnodata := FALSE;
END IF;
IF ST_BandIsNodata(rast1, newband1) THEN
RETURN rast1;
END IF;
END IF;
IF newkeepsourcenodata IS NULL THEN
newkeepsourcenodata := FALSE;
END IF;
IF newkeepsourcenodata THEN
IF ST_BandHasNodataValue(rast2, newband2) THEN
rast2nodataval := ST_BandNoDataValue(rast2, newband2);
ELSE
RAISE NOTICE 'ST_SetValues: keepsourcenodata was set to true but rast2 does not have a nodata value. keepsourcenodata reset to false';
newkeepsourcenodata := FALSE;
END IF;
IF ST_BandIsNodata(rast2, newband2) THEN
RETURN rast1;
END IF;
END IF;
IF newkeepdestnodata THEN
FOR cx IN newx..newx + newwidth - 1 LOOP
FOR cy IN newy..newy + newheight - 1 LOOP
r2val := ST_Value(rast2, newband2, cx - x2 + 1, cy - y2 + 1);
--RAISE NOTICE '444 x=%, y=%', cx, cy;
IF (ST_Value(newraster, newband1, cx, cy) != rast1nodataval) AND (NOT newkeepsourcenodata OR r2val != rast2nodataval) THEN
newraster := ST_SetValue(newraster, newband1, cx, cy, r2val);
END IF;
END LOOP;
END LOOP;
ELSE
--RAISE NOTICE '555 newx + newwidth - 1=%, newy + newheight - 1=%', newx + newwidth - 1, newy + newheight - 1;
FOR cx IN newx..newx + newwidth - 1 LOOP
FOR cy IN newy..newy + newheight - 1 LOOP
r2val := ST_Value(rast2, newband2, cx - x2 + 1, cy - y2 + 1);
--RAISE NOTICE '666 x=%, y=%', cx, cy;
IF NOT newkeepsourcenodata OR r2val != rast2nodataval THEN
newraster := ST_SetValue(newraster, newband1, cx, cy, r2val);
END IF;
END LOOP;
END LOOP;
END IF;
RETURN newraster;
END;
$$
LANGUAGE 'plpgsql';
-- Variant with both band = 1
CREATE OR REPLACE FUNCTION ST_SetValues(rast1 raster, x int, y int, width int, height int, rast2 raster, keepdestnodata boolean, keepsourcenodata boolean)
RETURNS raster
AS 'SELECT ST_SetValues($1, 1, $2, $3, $4, $5, $6, 1, $7, $8)'
LANGUAGE 'SQL' IMMUTABLE;
-- Variant with both band = 1 & both keepnodata = FALSE
CREATE OR REPLACE FUNCTION ST_SetValues(rast1 raster, x int, y int, width int, height int, rast2 raster)
RETURNS raster
AS 'SELECT ST_SetValues($1, 1, $2, $3, $4, $5, $6, 1, FALSE, FALSE)'
LANGUAGE 'SQL' IMMUTABLE;
-- Test
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_TestRaster(0, 0, 1)) as pix) foo
SELECT asbinary((pix).geom), (pix).val
FROM (SELECT ST_PixelAsPolygons(ST_SetValues(ST_TestRaster(0, 0, 1), 2, 1, 3, 1, ST_TestRaster(3, 0, 3))) as pix) foo

27
raster/scripts/plpgsql/st_sevalues.sql
View file

@ -1,27 +0,0 @@
-- $Id$
----------------------------------------------------------------------
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
----------------------------------------------------------------------
CREATE OR REPLACE FUNCTION ST_SetValues(rast raster, band, x int, y int, width int, height int, val float8)
RETURNS raster AS
$$
DECLARE
newraster raster := rast;
cx int;
cy int;
newwidth int := CASE WHEN x + width > ST_Width(rast) THEN ST_Width(rast) - x ELSE width END;
newheight int := CASE WHEN y + height > ST_Height(rast) THEN ST_Height(rast) - y ELSE height END;
BEGIN
newrast
FOR cx IN 1..newwidth LOOP
FOR cy IN 1..newheight LOOP
newrast := ST_SetValue(newrast, band, cx, cy, val);
END LOOP;
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';

55
raster/scripts/plpgsql/st_union.sql
View file

@ -1,8 +1,13 @@
----------------------------------------------------------------------
--
-- $Id$
----------------------------------------------------------------------
--
-- Copyright (c) 2009-2010 Pierre Racine <pierre.racine@sbf.ulaval.ca>
--
----------------------------------------------------------------------
-- Note: The functions provided in this script are in developement. Do not use.
CREATE OR REPLACE FUNCTION ST_MultiBandMapAlgebra(rast1 raster,
rast2 raster,
expression text,
@ -29,54 +34,6 @@ CREATE OR REPLACE FUNCTION ST_MultiBandMapAlgebra(rast1 raster,
$$
LANGUAGE 'plpgsql';
CREATE OR REPLACE FUNCTION ST_AddBand(rast1 raster, rast2 raster, band int, index int)
RETURNS raster AS
$$
DECLARE
newraster raster;
newnodatavalue int;
newindex int := index;
newband int := band;
x int;
y int;
BEGIN
IF ST_Width(rast1) != ST_Width(rast2) OR ST_Height(rast1) != ST_Height(rast2) THEN
RAISE EXCEPTION 'ST_AddBand: Attempting to add a band with different width or height';
END IF;
IF newindex < 1 THEN
newindex := 1;
END IF;
IF newindex IS NULL OR newindex > ST_NumBands(rast1) THEN
newindex := ST_NumBands(rast1) + 1;
END IF;
IF newband < 1 THEN
newband := 1;
END IF;
IF newband IS NULL OR newband > ST_NumBands(rast2) THEN
newband := ST_NumBands(rast2);
END IF;
IF newband = 0 THEN
RETURN rast1;
END IF;
raise notice 'aaaa newband=% newindex=%', newband, newindex;
newnodatavalue := ST_BandNodataValue(rast2, newband);
raise notice 'bbbb';
newraster := ST_AddBand(rast1, newindex, ST_BandPixelType(rast2, newband), newnodatavalue, newnodatavalue);
raise notice 'cccc';
FOR x IN 1..ST_Width(rast2) LOOP
FOR y IN 1..ST_Height(rast2) LOOP
newraster := ST_SetValue(newraster, newindex, x, y, ST_Value(rast2, newband, x, y));
END LOOP;
END LOOP;
RETURN newraster;
END;
$$
LANGUAGE 'plpgsql';
--Test
SELECT ST_NumBands(ST_AddBand(ST_TestRaster(0, 0, 1), ST_TestRaster(0, 0, 1), 1, 2))
CREATE OR REPLACE FUNCTION MultiBandMapAlgebra4Union(rast1 raster, rast2 raster, expression text)
RETURNS raster
AS 'SELECT ST_MultiBandMapAlgebra(rast1, rast2, expression, NULL, 'UNION')'