From: Sandro Santilli Date: Wed, 3 Sep 2014 14:52:33 +0000 (+0000) Subject: Implement ST_Retile and ST_CreateOverview (#2247) X-Git-Tag: 2.2.0rc1~863 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=9d07a9bde19223d88459a4c4a6babe9c43b00d0e;p=postgis Implement ST_Retile and ST_CreateOverview (#2247) Includes testcases and documentation git-svn-id: http://svn.osgeo.org/postgis/trunk@12941 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/NEWS b/NEWS index a6ba0998d..f1ddc0d31 100644 --- a/NEWS +++ b/NEWS @@ -17,6 +17,8 @@ PostGIS 2.2.0 * New Features * + - #2247, ST_Retile and ST_CreateOverview: in-db raster overviews creation + (Sandro Santilli / Vizzuality) - #899, -m shp2pgsql attribute names mapping -m switch (Regina Obe / Sandro Santilli) - #1678, Added GUC postgis.gdal_datapath to specify GDAL config diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index 263baf631..7b57052e0 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -1323,6 +1323,108 @@ WHERE short_name = 'GTiff') As g; + + + + + + ST_Retile + + Return a set of configured tiles from an arbitrarily tiled raster coverage. + + + + + + + + SETOF raster ST_Retile + regclass tab + name col + geometry ext + float8 sfx + float8 sfy + int tw + int th + text algo='NearestNeighbor' + + + + + + + Description + + +Return a set of tiles having the specified scale (sfx, +sfy) and max size (tw, +th) and covering the specified extent +(ext) with data coming from the specified +raster coverage (tab, col). + + + Algorithm options are: 'NearestNeighbor', 'Bilinear', 'Cubic', 'CubicSpline', and 'Lanczos'. Refer to: GDAL Warp resampling methods for more details. + + Availability: 2.2.0 + + + See Also + + + + + + + + + ST_CreateOverview + +Create an reduced resolution version of a given raster coverage. + + + + + + + + regclass ST_CreateOverview + regclass tab + name col + int factor + text algo='NearestNeighbor' + + + + + + + Description + + +Create an overview table with resampled tiles from the source table. +Output tiles will have the same size of input tiles and cover the same +spatial extent with a lower resolution (pixel size will be +1/factor of the original in both directions). + + + +The overview table will be made available in the +raster_overviews catalog and will have raster +constraints enforced. + + + Algorithm options are: 'NearestNeighbor', 'Bilinear', 'Cubic', 'CubicSpline', and 'Lanczos'. Refer to: GDAL Warp resampling methods for more details. + + Availability: 2.2.0 + + + See Also + + , + , + , + + diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in index 3a9255981..bbba976e7 100644 --- a/raster/rt_pg/rtpostgis.sql.in +++ b/raster/rt_pg/rtpostgis.sql.in @@ -8581,6 +8581,144 @@ CREATE OR REPLACE FUNCTION UpdateRasterSRID( AS $$ SELECT _UpdateRasterSRID('', $1, $2, $3) $$ LANGUAGE 'sql' VOLATILE STRICT; +------------------------------------------------------------------------------ +-- ST_Retile +------------------------------------------------------------------------------ + +-- Availability: 2.2.0 +-- @param ext extent to create overviews for, also used for grid origin +-- SRID must match source tile srid. +-- @param sfx scale factor x (pixel width) +-- @param sfy scale factor y (pixel height, usually negative) +-- @param tw max tile width +-- @param th max tile height +-- +CREATE OR REPLACE FUNCTION ST_Retile(tab regclass, col name, ext geometry, sfx float8, sfy float8, tw int, th int, algo text DEFAULT 'NearestNeighbour') +RETURNS SETOF raster AS $$ +DECLARE + rec RECORD; + ipx FLOAT8; + ipy FLOAT8; + tx int; + ty int; + te GEOMETRY; -- tile extent + ncols int; + nlins int; + srid int; + sql TEXT; +BEGIN + + RAISE DEBUG 'Target coverage will have sfx=%, sfy=%', sfx, sfy; + + -- 2. Loop over each target tile and build it from source tiles + ipx := st_xmin(ext); + ncols := ceil((st_xmax(ext)-ipx)/sfx/tw); + IF sfy < 0 THEN + ipy := st_ymax(ext); + nlins := ceil((st_ymin(ext)-ipy)/sfy/th); + ELSE + ipy := st_ymin(ext); + nlins := ceil((st_ymax(ext)-ipy)/sfy/th); + END IF; + + srid := ST_Srid(ext); + + RAISE DEBUG 'Target coverage will have % x % tiles, each of approx size % x %', ncols, nlins, tw, th; + RAISE DEBUG 'Target coverage will cover extent %', ext::box2d; + + FOR tx IN 0..ncols-1 LOOP + FOR ty IN 0..nlins-1 LOOP + te := ST_MakeEnvelope(ipx + tx * tw * sfx, + ipy + ty * th * sfy, + ipx + (tx+1) * tw * sfx, + ipy + (ty+1) * th * sfy, + srid); + --RAISE DEBUG 'sfx/sfy: %, %', sfx, sfy; + --RAISE DEBUG 'tile extent %', te; + sql := 'SELECT count(*), ST_Clip(ST_Union(ST_SnapToGrid(ST_Rescale(ST_Clip(' || quote_ident(col) + || ', st_expand($3, greatest($1,$2))),$1, $2, $6), $4, $5, $1, $2)), $3) g FROM ' || tab::text + || ' WHERE ST_Intersects(' || quote_ident(col) || ', $3)'; + --RAISE DEBUG 'SQL: %', sql; + FOR rec IN EXECUTE sql USING sfx, sfy, te, ipx, ipy, algo LOOP + --RAISE DEBUG '% source tiles intersect target tile %,% with extent %', rec.count, tx, ty, te::box2d; + IF rec.g IS NULL THEN + RAISE WARNING 'No source tiles cover target tile %,% with extent %', + tx, ty, te::box2d; + ELSE + --RAISE DEBUG 'Tile for extent % has size % x %', te::box2d, st_width(rec.g), st_height(rec.g); + RETURN NEXT rec.g; + END IF; + END LOOP; + END LOOP; + END LOOP; + + RETURN; +END; +$$ LANGUAGE 'plpgsql' STABLE STRICT; + +------------------------------------------------------------------------------ +-- ST_CreateOverview +------------------------------------------------------------------------------ + +-- Availability: 2.2.0 +CREATE OR REPLACE FUNCTION ST_CreateOverview(tab regclass, col name, factor int, algo text DEFAULT 'NearestNeighbour') +RETURNS regclass AS $$ +DECLARE + sinfo RECORD; -- source info + sql TEXT; + ttab TEXT; +BEGIN + + -- 0. Check arguments, we need to ensure: + -- a. Source table has a raster column with given name + -- b. Source table has a fixed scale (or "factor" would have no meaning) + -- c. Source table has a known extent ? (we could actually compute it) + -- d. Source table has a fixed tile size (or "factor" would have no meaning?) + -- # all of the above can be checked with a query to raster_columns + sql := 'SELECT r.r_table_schema sch, r.r_table_name tab, ' + || 'r.scale_x sfx, r.scale_y sfy, r.blocksize_x tw, ' + || 'r.blocksize_y th, r.extent ext, r.srid FROM raster_columns r, ' + || 'pg_class c, pg_namespace n WHERE r.r_table_schema = n.nspname ' + || 'AND r.r_table_name = c.relname AND r_raster_column = $2 AND ' + || ' c.relnamespace = n.oid AND c.oid = $1' + ; + EXECUTE sql INTO sinfo USING tab, col; + IF sinfo IS NULL THEN + RAISE EXCEPTION '%.% raster column does not exist', tab::text, col; + END IF; + IF sinfo.sfx IS NULL or sinfo.sfy IS NULL THEN + RAISE EXCEPTION 'cannot create overview without scale constraint, try select AddRasterConstraints(''%'', ''%'');', tab::text, col; + END IF; + IF sinfo.tw IS NULL or sinfo.tw IS NULL THEN + RAISE EXCEPTION 'cannot create overview without tilesize constraint, try select AddRasterConstraints(''%'', ''%'');', tab::text, col; + END IF; + IF sinfo.ext IS NULL THEN + RAISE EXCEPTION 'cannot create overview without extent constraint, try select AddRasterConstraints(''%'', ''%'');', tab::text, col; + END IF; + + -- TODO: lookup in raster_overviews to see if there's any + -- lower-resolution table to start from + + ttab := 'o_' || factor || '_' || sinfo.tab; + sql := 'CREATE TABLE ' || quote_ident(sinfo.sch) + || '.' || quote_ident(ttab) + || ' AS SELECT ST_Retile($1, $2, $3, $4, $5, $6, $7) ' + || quote_ident(col); + EXECUTE sql USING tab, col, sinfo.ext, + sinfo.sfx * factor, sinfo.sfy * factor, + sinfo.tw, sinfo.th, algo; + + -- TODO: optimize this using knowledge we have about + -- the characteristics of the target column ? + PERFORM AddRasterConstraints(sinfo.sch, ttab, col); + + PERFORM AddOverviewConstraints(sinfo.sch, ttab, col, + sinfo.sch, sinfo.tab, col, factor); + + RETURN ttab; +END; +$$ LANGUAGE 'plpgsql' VOLATILE STRICT; + ------------------------------------------------------------------- -- Debugging ------------------------------------------------------------------- diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index fe9892f10..6b38234fb 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -109,7 +109,8 @@ TEST_UTILITY = \ rt_reclass \ rt_gdalwarp \ rt_asraster \ - rt_dumpvalues + rt_dumpvalues \ + rt_createoverview TEST_MAPALGEBRA = \ rt_mapalgebraexpr \ diff --git a/raster/test/regress/rt_createoverview.sql b/raster/test/regress/rt_createoverview.sql new file mode 100644 index 000000000..3cc172ca0 --- /dev/null +++ b/raster/test/regress/rt_createoverview.sql @@ -0,0 +1,54 @@ +SET client_min_messages TO warning; +CREATE TABLE res1 AS SELECT + ST_AddBand( + ST_MakeEmptyRaster(10, 10, x, y, 1, -1, 0, 0, 0) + , 1, '8BUI', 0, 0 + ) r +FROM generate_series(-170, 160, 10) x, + generate_series(80, -70, -10) y; +SELECT addrasterconstraints('res1', 'r'); + +SELECT ST_CreateOverview('res1', 'r', 2)::text = 'o_2_res1'; + +SELECT ST_CreateOverview('res1', 'r', 4)::text = 'o_4_res1'; + +SELECT ST_CreateOverview('res1', 'r', 8)::text = 'o_8_res1'; + +SELECT ST_CreateOverview('res1', 'r', 16)::text = 'o_16_res1'; + +SELECT r_table_name tab, r_raster_column c, srid s, + scale_x sx, scale_y sy, + blocksize_x w, blocksize_y h, same_alignment a, + -- regular_blocking (why not regular?) + --extent::box2d e, + st_covers(extent::box2d, 'BOX(-170 -80,170 80)'::box2d) ec, + st_xmin(extent::box2d) = -170 as eix, + st_ymax(extent::box2d) = 80 as eiy, + (st_xmax(extent::box2d) - 170) <= scale_x as eox, + --(st_xmax(extent::box2d) - 170) eoxd, + abs(st_ymin(extent::box2d) + 80) <= abs(scale_y) as eoy + --,abs(st_ymin(extent::box2d) + 80) eoyd + FROM raster_columns +WHERE r_table_name like '%res1' +ORDER BY scale_x, r_table_name; + +SELECT o_table_name, o_raster_column, + r_table_name, r_raster_column, overview_factor +FROM raster_overviews +WHERE r_table_name = 'res1' +ORDER BY overview_factor; + +SELECT 'count', +(SELECT count(*) r1 from res1), +(SELECT count(*) r2 from o_2_res1), +(SELECT count(*) r4 from o_4_res1), +(SELECT count(*) r8 from o_8_res1), +(SELECT count(*) r16 from o_16_res1) +; + + +DROP TABLE o_16_res1; +DROP TABLE o_8_res1; +DROP TABLE o_4_res1; +DROP TABLE o_2_res1; +DROP TABLE res1; diff --git a/raster/test/regress/rt_createoverview_expected b/raster/test/regress/rt_createoverview_expected new file mode 100644 index 000000000..7173bcaa5 --- /dev/null +++ b/raster/test/regress/rt_createoverview_expected @@ -0,0 +1,15 @@ +t +t +t +t +t +res1|r|0|1|-1|10|10|t|t|t|t|t|t +o_2_res1|r|0|2|-2|10|10|t|t|t|t|t|t +o_4_res1|r|0|4|-4|10|10|t|t|t|t|t|t +o_8_res1|r|0|8|-8|10|10|t|t|t|t|t|t +o_16_res1|r|0|16|-16|10|10|t|t|t|t|t|t +o_2_res1|r|res1|r|2 +o_4_res1|r|res1|r|4 +o_8_res1|r|res1|r|8 +o_16_res1|r|res1|r|16 +count|544|136|36|10|3