From 37b702b7a7697a0e890300c9cb8e26f391c19d97 Mon Sep 17 00:00:00 2001 From: Pierre Racine Date: Thu, 1 Dec 2011 15:52:39 +0000 Subject: [PATCH] Ticket 1342. Integrate ST_Clip script into rtpostgis.sql git-svn-id: http://svn.osgeo.org/postgis/trunk@8283 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/rt_pg/rtpostgis.sql.in.c | 122 ++++++++++++++++++++++++++++++++ 1 file changed, 122 insertions(+) diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 079c23c70..cf022490c 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -2712,6 +2712,35 @@ CREATE OR REPLACE FUNCTION st_raster2worldcoordy(rast raster, yr int) $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT; +----------------------------------------------------------------------- +-- ST_MinPossibleVal(pixeltype text) +-- Return the smallest value for a given pixeltyp. +-- Should be called like this: +-- SELECT ST_MinPossibleVal(ST_BandPixelType(rast, band)) +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION ST_MinPossibleVal(pixeltype text) + RETURNS float8 AS + $$ + DECLARE + newval int := 0; + BEGIN + newval := CASE + WHEN pixeltype = '1BB' THEN 0 + WHEN pixeltype = '2BUI' THEN 0 + WHEN pixeltype = '4BUI' THEN 0 + WHEN pixeltype = '8BUI' THEN 0 + WHEN pixeltype = '8BSI' THEN -128 + WHEN pixeltype = '16BUI' THEN 0 + WHEN pixeltype = '16BSI' THEN -32768 + WHEN pixeltype = '32BUI' THEN 0 + WHEN pixeltype = '32BSI' THEN -2147483648 + WHEN pixeltype = '32BF' THEN -2147483648 -- Could not find a function returning the smallest real yet + WHEN pixeltype = '64BF' THEN -2147483648 -- Could not find a function returning the smallest float8 yet + END; + RETURN newval; + END; + $$ + LANGUAGE 'plpgsql'; ----------------------------------------------------------------------- -- Raster Outputs @@ -3387,6 +3416,99 @@ CREATE AGGREGATE ST_Union(raster) ( FINALFUNC = MapAlgebra4UnionFinal1 ); +------------------------------------------------------------------- +-- ST_Clip(rast raster, band int, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false) +-- Clip the values of a raster to the shape of a polygon. +-- +-- rast - raster to be clipped +-- band - limit the result to only one band +-- geom - geometry defining the she to clip the raster +-- nodata - define (if there is none defined) or replace the raster nodata value with this value +-- trimraster - limit the extent of the result to the extent of the geometry +-- Todo: +-- test point +-- test line +-- test polygon smaller than pixel +-- test and optimize raster totally included in polygon +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; + newpixtype text; + 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 + RAISE NOTICE 'Raster do not have band %. Returning null', band; + RETURN null; + ELSE + bandstart := band; + bandend := band; + END IF; + newpixtype := ST_BandPixelType(rast, bandstart); + newnodata := coalesce(nodata, ST_BandNodataValue(rast, bandstart), ST_MinPossibleVal(newpixtype)); + newextent := CASE WHEN trimraster THEN 'INTERSECTION' ELSE 'FIRST' END; + +--RAISE NOTICE 'newextent=%', newextent; + -- Convert the geometry to a raster + geomrast := ST_AsRaster(geom, rast, ST_BandPixelType(rast, band), 1, newnodata); + + -- Set the newnodata + sourceraster := ST_SetBandNodataValue(sourceraster, bandstart, newnodata); + + -- Compute the first raster band + newrast := ST_MapAlgebraExpr(sourceraster, bandstart, geomrast, 1, 'rast1', newpixtype, newextent); + + FOR bandi IN bandstart+1..bandend LOOP +--RAISE NOTICE 'bandi=%', bandi; + -- for each band we must determine the nodata value + newpixtype := ST_BandPixelType(rast, bandi); + newnodata := coalesce(nodata, ST_BandNodataValue(sourceraster, bandi), ST_MinPossibleVal(newpixtype)); + sourceraster := ST_SetBandNodataValue(sourceraster, bandi, newnodata); + newrast := ST_AddBand(newrast, ST_MapAlgebraExpr(sourceraster, bandi, geomrast, 1, 'rast1', newpixtype, newextent)); + END LOOP; + RETURN newrast; + END; + $$ + LANGUAGE 'plpgsql'; + +-- Variant defaulting to band 1 +CREATE OR REPLACE FUNCTION ST_Clip(rast raster, geom geometry, nodata float8 DEFAULT null, trimraster boolean DEFAULT false) + RETURNS raster + AS $$ + SELECT ST_Clip($1, 1, $2, $3, $4); + $$ LANGUAGE 'SQL'; + +-- Variant defaulting nodata to the one of the raster or the min possible value +CREATE OR REPLACE FUNCTION ST_Clip(rast raster, band int, geom geometry, trimraster boolean) + RETURNS raster + AS $$ + SELECT ST_Clip($1, $2, $3, null, $4); + $$ LANGUAGE 'SQL'; + +-- Variant defaulting nodata to the one of the raster or the min possible value +CREATE OR REPLACE FUNCTION ST_Clip(rast raster, geom geometry, trimraster boolean) + RETURNS raster + AS $$ + SELECT ST_Clip($1, 1, $2, null, $3); + $$ LANGUAGE 'SQL'; + ------------------------------------------------------------------------------ -- RASTER_COLUMNS -- -- 2.40.0