Commit 9d07a9bd authored by Sandro Santilli's avatar Sandro Santilli

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
parent 4ae13a35
......@@ -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
......
......@@ -1323,6 +1323,108 @@ WHERE short_name = 'GTiff') As g;
<para>
<xref linkend="UpdateGeometrySRID"/>
</para>
</refsection>
</refentry>
<refentry id="RT_Retile">
<refnamediv>
<refname>ST_Retile</refname>
<refpurpose>
Return a set of configured tiles from an arbitrarily tiled raster coverage.
</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
<funcdef>SETOF raster <function>ST_Retile</function></funcdef>
<paramdef><type>regclass </type> <parameter>tab</parameter></paramdef>
<paramdef><type>name </type> <parameter>col</parameter></paramdef>
<paramdef><type>geometry </type> <parameter>ext</parameter></paramdef>
<paramdef><type>float8 </type> <parameter>sfx</parameter></paramdef>
<paramdef><type>float8 </type> <parameter>sfy</parameter></paramdef>
<paramdef><type>int </type> <parameter>tw</parameter></paramdef>
<paramdef><type>int </type> <parameter>th</parameter></paramdef>
<paramdef choice="opt"><type>text </type> <parameter>algo='NearestNeighbor'</parameter></paramdef>
</funcprototype>
</funcsynopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
<para>
Return a set of tiles having the specified scale (<varname>sfx</varname>,
<varname>sfy</varname>) and max size (<varname>tw</varname>,
<varname>th</varname>) and covering the specified extent
(<varname>ext</varname>) with data coming from the specified
raster coverage (<varname>tab</varname>, <varname>col</varname>).
</para>
<para>Algorithm options are: 'NearestNeighbor', 'Bilinear', 'Cubic', 'CubicSpline', and 'Lanczos'. Refer to: <ulink url="http://www.gdal.org/gdalwarp.html">GDAL Warp resampling methods</ulink> for more details.</para>
<para>Availability: 2.2.0</para>
</refsection>
<refsection>
<title>See Also</title>
<para>
<xref linkend="RT_CreateOverview"/>
</para>
</refsection>
</refentry>
<refentry id="RT_CreateOverview">
<refnamediv>
<refname>ST_CreateOverview</refname>
<refpurpose>
Create an reduced resolution version of a given raster coverage.
</refpurpose>
</refnamediv>
<refsynopsisdiv>
<funcsynopsis>
<funcprototype>
<funcdef>regclass <function>ST_CreateOverview</function></funcdef>
<paramdef><type>regclass </type> <parameter>tab</parameter></paramdef>
<paramdef><type>name </type> <parameter>col</parameter></paramdef>
<paramdef><type>int </type> <parameter>factor</parameter></paramdef>
<paramdef choice="opt"><type>text </type> <parameter>algo='NearestNeighbor'</parameter></paramdef>
</funcprototype>
</funcsynopsis>
</refsynopsisdiv>
<refsection>
<title>Description</title>
<para>
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/<varname>factor</varname> of the original in both directions).
</para>
<para>
The overview table will be made available in the
<varname>raster_overviews</varname> catalog and will have raster
constraints enforced.
</para>
<para>Algorithm options are: 'NearestNeighbor', 'Bilinear', 'Cubic', 'CubicSpline', and 'Lanczos'. Refer to: <ulink url="http://www.gdal.org/gdalwarp.html">GDAL Warp resampling methods</ulink> for more details.</para>
<para>Availability: 2.2.0</para>
</refsection>
<refsection>
<title>See Also</title>
<para>
<xref linkend="RT_Retile"/>,
<xref linkend="RT_AddOverviewConstraints"/>,
<xref linkend="RT_AddRasterConstraints"/>,
<xref linkend="RT_Raster_Overviews"/>
</para>
</refsection>
</refentry>
</sect1>
......
......@@ -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
-------------------------------------------------------------------
......
......@@ -109,7 +109,8 @@ TEST_UTILITY = \
rt_reclass \
rt_gdalwarp \
rt_asraster \
rt_dumpvalues
rt_dumpvalues \
rt_createoverview
TEST_MAPALGEBRA = \
rt_mapalgebraexpr \
......
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;
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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment