From fff1c06741963c0ee721adf881f5721ed01d40ff Mon Sep 17 00:00:00 2001 From: Pierre Racine Date: Tue, 7 Feb 2012 18:40:08 +0000 Subject: [PATCH] Added ST_TileAsGeom() returning only the extent of the planned tiles as polygons git-svn-id: http://svn.osgeo.org/postgis/trunk@9069 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/scripts/plpgsql/st_tile.sql | 91 +++++++++++++++++++++++++++--- 1 file changed, 82 insertions(+), 9 deletions(-) diff --git a/raster/scripts/plpgsql/st_tile.sql b/raster/scripts/plpgsql/st_tile.sql index ee96a9ce9..d7eee98ed 100644 --- a/raster/scripts/plpgsql/st_tile.sql +++ b/raster/scripts/plpgsql/st_tile.sql @@ -1,4 +1,4 @@ ----------------------------------------------------------------------- +---------------------------------------------------------------------- -- -- $Id: st_tile.sql 8255 2011-11-29 16:34:48Z pracine $ -- @@ -6,7 +6,11 @@ -- ---------------------------------------------------------------------- -- ST_Tile --- Split a raster into a set of raster tiles, one tile per row returned. +-- Split a raster into a set of raster tiles, one tile per row returned. +-- Works on multiband rasters. There is no way to specify the upper left +-- corner of the new tiled grid. The grid start at the upperleft corner +-- of the provided raster. +-- -- rast - Raster to be tiled. -- width - Width of the tiles. -- height - Height of the tiles @@ -17,7 +21,7 @@ -- nodatavalue - nodata value to use to pad the outbound tiles when the provided -- raster do not have a nodata value defined. If not provided and -- the raster do not have a nodata value defined --- ST_MinPossibleValue(ST_BandPixelType(rast, band)) is used. +-- ST_MinPossibleValue(ST_BandPixelType(rast, band)) is used for each band. -- -- Example producing 120 x 120 pixel tiles -- @@ -25,7 +29,7 @@ -- SELECT ST_Tile(rast, 120, 120, true), generate_series(1, 3600) rid -- FROM srtm_22_03; ---------------------------------------------------------------------------------------------------------------------- -DROP FUNCTION ST_Tile(rast raster, width integer, height integer, padwithnodata boolean, nodatavalue double precision); +DROP FUNCTION IF EXISTS ST_Tile(rast raster, width integer, height integer, padwithnodata boolean, nodatavalue double precision); CREATE OR REPLACE FUNCTION ST_Tile(rast raster, width integer, height integer, padwithnodata boolean DEFAULT FALSE, nodatavalue double precision DEFAULT NULL) RETURNS SETOF raster AS $$ @@ -38,7 +42,7 @@ CREATE OR REPLACE FUNCTION ST_Tile(rast raster, width integer, height integer, p newnodata double precision; newpixtype text; nbband integer; - bnadi integer; + bandi integer; newrast raster; initvalue double precision; grid record; @@ -98,6 +102,75 @@ CREATE OR REPLACE FUNCTION ST_Tile(rast raster, width integer, height integer, p END; $$ LANGUAGE 'plpgsql'; + +---------------------------------------------------------------------- +-- ST_TileAsGeom +-- Split a raster into a set of raster tiles, returning only the geometry +-- corresponding to each tile. +-- There is no way to specify the upper left corner of the new tiled grid. +-- The grid start at the upperleft corner of the provided raster. +-- +-- rast - Raster to be tiled. +-- width - Width of the tiles. +-- height - Height of the tiles +-- padwithnodata - If TRUE, the produced tiles are strictly width x heigth pixels. +-- Pixels outside the extent of the passed raster are filled with +-- nodata value. When FALSE out of bound tiles are clipped to the +-- extent of the raster. Default to FALSE. +-- +-- Example producing 120 x 120 pixel tiles +-- +-- SELECT ST_TileAsGeom(rast, 130, 130, true) +-- FROM srtm_22_03; +---------------------------------------------------------------------------------------------------------------------- +DROP FUNCTION IF EXISTS ST_TileAsGeom(rast raster, width integer, height integer, padwithnodata boolean); +CREATE OR REPLACE FUNCTION ST_TileAsGeom(rast raster, width integer, height integer, padwithnodata boolean DEFAULT FALSE) + RETURNS SETOF geometry AS + $$ + DECLARE + gridrast raster; + rastwidth integer; + rastheight integer; + gridwidth integer; + gridheight integer; + nbband integer; + rastextent geometry; + BEGIN + IF rast IS NULL THEN + RETURN; + END IF; + + nbband := ST_Numbands(rast); + IF nbband < 1 THEN + RAISE NOTICE 'Raster do not have band %. Returning null', band; + RETURN; + END IF; + + rastwidth := ST_Width(rast); + IF width IS NULL THEN + width := rastwidth; + END IF; + + rastheight := ST_Height(rast); + IF height IS NULL THEN + height := rastheight; + END IF; + + gridwidth := (rastwidth / width) + CASE WHEN rastwidth % width > 0 THEN 1 ELSE 0 END; + gridheight := (rastheight / height) + CASE WHEN rastheight % height > 0 THEN 1 ELSE 0 END; + + gridrast := ST_AddBand(ST_MakeEmptyRaster(gridwidth, gridheight, ST_UpperLeftX(rast), ST_UpperLeftY(rast), ST_ScaleX(rast) * width, ST_ScaleY(rast) * height, ST_SkewX(rast), ST_SkewY(rast), ST_SRID(rast)), '8BUI'::text, 1, 0); + IF padwithnodata THEN + RETURN QUERY SELECT (ST_PixelAsPolygons(gridrast)).geom geom; + ELSE + rastextent := rast::geometry; + RETURN QUERY SELECT ST_Intersection(rastextent, (ST_PixelAsPolygons(gridrast)).geom) geom; + END IF; + RETURN; + END; + $$ + LANGUAGE 'plpgsql'; + -- Redefine ST_TestRaster() CREATE OR REPLACE FUNCTION ST_TestRaster(ulx float8, uly float8, val float8) @@ -118,10 +191,6 @@ FROM (SELECT ST_PixelAsPolygons(ST_TestRaster(0, 0, 1)) gv) foo; -- Tile it to 10x10 tiles SELECT ST_Tile(ST_TestRaster(0, 0, 1), 10, 10) --- Old test for the extent of the tile grid -SELECT ST_AsBinary((gv).geom) geom, (gv).val -FROM (SELECT ST_PixelAsPolygons(ST_Rescale(ST_TestRaster(0, 0, 1), 0.003, 0.004)) gv) foo; - -- Display the result of the tile function SELECT ST_AsBinary((gv).geom) geom, (gv).val FROM (SELECT ST_PixelAsPolygons(ST_Tile(ST_TestRaster(0, 0, 1), 3, 4)) gv) foo; @@ -146,6 +215,10 @@ SELECT ST_AsBinary((gv).geom) geom, (gv).val, rid FROM (SELECT ST_PixelAsPolygons(rast) gv, rid FROM (SELECT ST_Tile(ST_TestRaster(0, 0, 1), 10, 10, true) rast, generate_series(1, 35) rid) foo WHERE rid = 2) foo2; +-- Test ST_TileAsGeom +SELECT ST_TileAsGeom(ST_TestRaster(0, 0, 1), 10, 10); + + -- Other tests SELECT ST_Tile(ST_AddBand(ST_MakeEmptyRaster(48, 63, 0, 0, 0.001, -0.001, 0, 0, 4269), '8BSI'::text, 0, -2), 10, 10, true, -10); -- 2.40.0