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

plpgsql implementation for st_clip.sql(raster, geom)

Browse source 
git-svn-id: http://svn.osgeo.org/postgis/trunk@8267 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Pierre Racine 2011-11-30 22:42:13 +00:00
parent dab1e5887f
commit a90bc98922
1 changed files with 84 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -1,4 +1,4 @@
----------------------------------------------------------------------
----------------------------------------------------------------------
--
-- $Id$
--
@ -34,4 +34,87 @@ CREATE OR REPLACE FUNCTION ST_Clip(rast raster, x int, y int, width int, height
$$
LANGUAGE 'plpgsql';
CREATE OR REPLACE FUNCTION ST_Clip(rast raster, band int, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false)
RETURNS raster AS
$$
DECLARE
sourceraster raster := rast;
newrast raster;
geomrast raster;
numband int := ST_Numbands(rast);
bandstart int;
bandend int;
newextent text;
newnodata float8;
bandi int;
BEGIN
IF rast IS NULL THEN
RETURN null;
END IF;
IF geom IS NULL THEN
RETURN rast;
END IF;
IF band IS NULL THEN
bandstart := 1;
bandend := numband;
ELSEIF ST_HasNoBand(rast, band) THEN
bandstart := numband;
bandend := numband;
ELSE
bandstart := band;
bandend := band;
END IF;
IF nodata IS NULL THEN
newnodata := coalesce(ST_BandNodataValue(rast, bandstart), ST_MinPossibleVal(ST_BandPixelType(rast, bandstart)));
ELSE
newnodata := nodata;
END IF;
newextent := CASE WHEN trimraster THEN 'INTERSECTION' ELSE 'FIRST' END;
RAISE NOTICE 'bandstart=%, bandend=%', bandstart, bandend;
-- Convert the geometry to a raster
geomrast := ST_AsRaster(geom, rast, ST_BandPixelType(rast, band), 1, newnodata);
-- Compute the first raster band
sourceraster := ST_SetBandNodataValue(sourceraster, bandstart, newnodata);
newrast := ST_MapAlgebraExpr(sourceraster, bandstart, geomrast, 1, 'rast1');
FOR bandi IN bandstart+1..bandend LOOP
RAISE NOTICE 'bandi=%', bandi;
-- for each band we must determine the nodata value
IF nodata IS NULL THEN
newnodata := coalesce(ST_BandNodataValue(sourceraster, bandi), ST_MinPossibleVal(ST_BandPixelType(rast, bandi)));
ELSE
newnodata := nodata;
END IF;
sourceraster := ST_SetBandNodataValue(sourceraster, bandi, newnodata);
newrast := ST_AddBand(newrast, ST_MapAlgebraExpr(sourceraster, bandi, geomrast, 1, 'rast1'));
END LOOP;
RETURN newrast;
END;
$$
LANGUAGE 'plpgsql';
-- Test
CREATE OR REPLACE FUNCTION ST_TestRaster(h integer, w integer, val float8)
RETURNS raster AS
$$
DECLARE
BEGIN
RETURN ST_AddBand(ST_MakeEmptyRaster(h, w, 0, 0, 1, 1, 0, 0, 0), '32BF', val, -1);
END;
$$
LANGUAGE 'plpgsql';
SELECT ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4)) rast
-- Test one band raster
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_TestRaster(10, 10, 2), 1, ST_Buffer(ST_MakePoint(8, 5), 4))) gv
-- Test two bands raster
SELECT ST_AsBinary((gv).geom), (gv).val
FROM ST_PixelAsPolygons(ST_Clip(ST_AddBand(ST_TestRaster(10, 10, 2), '16BUI'::text, 4, 0), null, ST_Buffer(ST_MakePoint(8, 5), 4)), 2) gv