[postgis-users] Cut extent box2d into small regular boxes, and iterate through each box of features
Simon SPDBA Greener
simon at spdba.com.au
Tue Mar 29 02:16:34 PDT 2022
I think the following procedure and function will do what you want.
There is some documentation in the driving procedure.
-- Example
CALL spdba.QuadTree(
p_SearchOwner := 'data',
p_SearchTable := 'valves',
p_SearchColumn := 'geom',
p_LL :=
ST_GeomFromEWKT('SRID=28356;POINT(515698.10890000034 6960213.1757)'),
p_UR :=
ST_GeomFromEWKT('SRID=28356;POINT(519045.1911000004 6965208.943 )'),
p_TargetOwner := 'data',
p_TargetTable := 'valves_q',
p_TargetColumn := 'geom',
p_MaxQuadLevel := 8,
p_MaxCount := 100
);
regards
Simon
----------------------------------------------------------------------------------------------------------
DROP procedure if exists
spdba.QuadTree(Varchar(250),VARCHAR(250),VARCHAR(250),geometry,geometry,Varchar(250),varchar(250),VARCHAR(250),integer,integer);
Create or Replace Procedure spdba.QuadTree(
p_SearchOwner IN Varchar(250),
p_SearchTable IN VARCHAR(250),
p_SearchColumn IN VARCHAR(250),
p_LL IN geometry,
p_UR IN geometry,
p_TargetOwner IN Varchar(250),
p_TargetTable IN varchar(250),
p_TargetColumn IN VARCHAR(250),
p_MaxQuadLevel IN integer,
p_MaxCount IN integer
)
LANGUAGE 'plpgsql'
SECURITY DEFINER
As
/****f* TESSELATION/ST_QuadTree
* NAME
* ST_QuadTree - Tesselates a two-dimensional space using a simple
recursive quad tree grid.
* SYNOPSIS
* Procedure spdba.QuadTree(
* p_SearchOwner IN Varchar(250),
* p_SearchTable IN VARCHAR(250),
* p_SearchColumn IN VARCHAR(250),
* p_LL IN geometry,
* p_UR IN geometry,
* p_TargetOwner IN Varchar(250),
* p_TargetTable IN varchar(250),
* p_TargetColumn IN VARCHAR(250),
* p_MaxQuadLevel IN integer,
* p_MaxCount IN integer
* )
* DESCRIPTION
* Recursively tesselates the two-dimensional space defined by
p_SearchTable using a
* quad tree algorithm based on a set of criteria:
* 1. Depth of the Quad Tree;
* 2. Max number of features per quad
* (If number in a quad is > max number, quad is divided into
four and each quad feature count is recomputed)
* The ouput polygons representing the quads that contain the data
* are written to the p_TagetTable with some specific fields
* PARAMETERS
* p_SearchOwner - Varchar(250) - Schema owner of p_SearchTable table
* p_SearchTable - VARCHAR(250) - Name of table in p_SchemaOwner
that is to be quadded
* p_SearchColumn - VARCHAR(250) - Geometry column in p_SearchTable
containing spatial data to be quadded.
* p_LL - geometry - Lower Left corner of extent of
data in p_SearchColumn to be quadded (normally LL of extent of all data
in p_searchColumn)
* p_UR - geometry - Upper Right corner of extent of
data in p_SearchColumn to be quadded (normally UR of extent of all data
in p_searchColumn)
* p_TargetOwner - Varchar(250) - Schema owner of p_TargetTable
* p_TargetTable - varchar(250) - Name of table that will be
created and will hold the quad tree rectangles.
* p_TargetColumn - VARCHAR(250) - Name of geometry column in
p_TargetTable that will hold resultant quad rectangles.
* p_MaxQuadLevel - integer - Maximum depth to recurse.
* p_MaxCount - integer - Max number of features per quad
tree rectangle.
* EXAMPLE
* CALL spdba.QuadTree(
* p_SearchOwner := 'data',
* p_SearchTable := 'valves',
* p_SearchColumn := 'geom',
* p_LL :=
ST_GeomFromEWKT('SRID=28356;POINT(515698.10890000034 6960213.1757)'),
* p_UR :=
ST_GeomFromEWKT('SRID=28356;POINT(519045.1911000004 6965208.943 )'),
* p_TargetOwner := 'data',
* p_TargetTable := 'valves_q',
* p_TargetColumn := 'geom',
* p_MaxQuadLevel := 8,
* p_MaxCount := 200
* );
* NOTE
* Ignores Z and only supports geometry objects not geography
(separate procedure)
* HISTORY
* Simon Greener - March 2022 - Converted from Oracle PL/SQL
* COPYRIGHT
* Licensed under a Creative Commons Attribution-Share Alike 2.5
Australia License. (http://creativecommons.org/licenses/by-sa/2.5/au/)
******/
$OUTER$
DECLARE
v_sql text;
v_insert_sql text;
v_select_sql text;
v_result integer;
v_searchOwnerTable text := case when p_searchOwner is null then ''
else p_searchOwner || '.' end || p_SearchTable;
v_TargetOwnerTable text := case when p_TargetOwner is null then ''
else p_TargetOwner || '.' end || p_TargetTable;
BEGIN
/***** Inner subroutine*/
Drop Function IF EXISTS spdba.QuadTree(integer, geometry, geometry,
integer, integer, integer, text, text);
Create Function spdba.QuadTree(p_QuadId IN INTEGER,
p_LL IN geometry,
p_UR IN geometry,
p_QuadLevel IN INTEGER,
p_maxQuadLevel in integer,
p_MaxCount in integer,
p_insert_sql in text,
p_search_sql in text)
RETURNS INTEGER
LANGUAGE 'plpgsql'
AS
$INNER$
DECLARE
v_count INTEGER;
v_temp double precision;
v_QuadLevel INTEGER;
v_QuadID INTEGER;
v_LL geometry;
v_UR geometry;
BEGIN
v_QuadId := p_QuadId;
v_QuadLevel := p_QuadLevel;
v_LL := p_LL;
v_UR := p_UR;
-- Performance optimisation
-- When at level 0 don't execute a search but start tesselating
If ( v_QuadLevel <> 0 ) Then
EXECUTE p_search_sql
INTO v_count
USING
ST_X(p_LL)::float,ST_Y(p_LL)::float,ST_X(p_UR)::float,ST_Y(p_UR)::float,ST_Srid(p_LL)::integer;
IF ( v_count = 0 ) THEN
RETURN v_QuadId;
END IF;
IF ( v_count <= p_MaxCount ) THEN
EXECUTE p_insert_sql
USING p_QuadId,
v_QuadLevel,
v_count,
ST_X(v_LL), ST_Y(v_LL),
ST_X(v_UR), ST_Y(v_UR),
ST_MakeEnvelope(ST_X(v_LL)::float,ST_Y(v_LL)::float,ST_X(v_UR)::float,ST_Y(v_UR)::float,ST_Srid(p_LL)::integer);
RETURN v_QuadId + 1;
END IF;
End If;
v_QuadLevel := p_QuadLevel + 1;
IF ( v_QuadLevel > p_MaxQuadLevel ) THEN
RETURN v_QuadId;
END IF;
-- Still need to tessellate the space!
--
-- initialize tmp node with corner data
-- +---+---R
-- | | |
-- +---+---+
-- | x | |
-- L---+---+
v_UR := ST_SetSrid(
ST_MakePoint(ST_X(v_LL) + ( ST_X(v_UR) -
ST_X(v_LL) ) / 2.0,
ST_Y(v_LL) + ( ST_Y(v_UR) - ST_Y(v_LL) ) / 2.0),
ST_SRID(v_LL)
);
v_QuadId := spdba.QuadTree(v_QuadId,v_LL,v_UR,v_QuadLevel,
p_MaxQuadLevel,p_MaxCount,p_insert_sql,p_search_sql);
-- +---+---+
-- | x | |
-- +---R---+
-- | 1 | |
-- L---+---+
v_temp := ( ST_Y(v_UR) - ST_Y(v_LL) );
v_LL :=
ST_SetSrid(ST_MakePoint(ST_X(v_LL),ST_Y(v_UR)),ST_Srid(v_UR));
v_UR := ST_SetSrid(ST_MakePoint(ST_X(v_UR),ST_Y(v_LL) +
v_temp),ST_Srid(v_LL));
v_QuadId := spdba.QuadTree(v_QuadId,v_LL,v_UR,v_QuadLevel,
p_MaxQuadLevel,p_MaxCount,p_insert_sql,p_search_sql);
-- +---R---+
-- | 2 | |
-- L---+---+
-- | | x |
-- +---+---+
v_temp := ST_X(v_UR) - ST_X(v_LL);
v_LL := ST_SetSrid(ST_MakePoint(ST_X(v_UR),
ST_Y(v_LL)),ST_SRID(v_LL));
v_UR := ST_SetSrid(ST_MakePoint(ST_X(v_LL) +
v_temp,ST_Y(v_UR)),ST_SRID(v_LL));
v_temp := ST_Y(v_UR) - ST_Y(v_LL);
v_UR := ST_SetSRID(ST_MakePoint(ST_X(v_UR),ST_Y(v_LL)),
ST_Srid(v_LL));
v_LL :=
ST_SetSRID(ST_MakePoint(ST_X(v_LL),ST_Y(v_UR)-v_temp),ST_Srid(v_LL));
v_QuadId :=
spdba.QuadTree(v_QuadId,v_LL,v_UR,v_QuadLevel,p_MaxQuadLevel,p_MaxCount,p_insert_sql,p_search_sql);
-- +---+---+
-- | | x |
-- +---+---U
-- | | 3 |
-- +---L---+
v_temp := ST_Y(v_UR) - ST_Y(v_LL);
v_LL := ST_SetSrid(ST_MakePoint(ST_X(v_LL),ST_Y(v_UR)),
ST_Srid(v_UR));
v_UR :=
ST_SetSrid(ST_MakePoint(ST_X(v_UR),ST_Y(v_LL)+v_temp),ST_Srid(v_ll));
RETURN spdba.QuadTree(v_QuadId,v_LL,v_UR,v_QuadLevel,
p_MaxQuadLevel,p_MaxCount,p_insert_sql,p_search_sql);
END;
$INNER$;
/***************************/
/***** End Inner subroutine*/
/***************************/
IF ( p_LL is null or p_UR is null ) THEN
RAISE EXCEPTION 'Starting LL and UR must not be null';
END IF;
v_insert_sql := 'INSERT INTO ' || v_TargetOwnerTable ||
' (quad_id,quad_level,feature_count,xlo,ylo,xhi,yhi,'
|| p_TargetColumn || ') ' ||
' VALUES ($1,$2,$3,$4,$5,$6,$7,$8)';
v_select_sql := 'SELECT count(*)' ||
' FROM '|| v_SearchOwnerTable || ' as a ' ||
' WHERE ST_Intersects(a.' || p_SearchColumn ||
',ST_MakeEnvelope($1, $2, $3, $4, $5))';
-- Drop target table.
EXECUTE 'DROP TABLE IF EXISTS ' || v_TargetOwnerTable;
-- CreateTargetTable
v_sql := 'CREATE TABLE ' || v_TargetOwnerTable || '(
quad_id bigint NOT NULL,
quad_level integer,
feature_Count integer,
xlo float,
ylo float,
xhi float,
yhi float,
'|| p_TargetColumn || ' geometry)';
EXECUTE v_sql;
v_result := spdba.QuadTree( 0, p_LL, p_UR, 0, p_MaxQuadLevel,
p_MaxCount, v_insert_sql, v_select_sql );
return;
END;
$OUTER$;
More information about the postgis-users
mailing list