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