Thanks. Simon, I will have a close look and study.
Regards, David On Tue, 29 Mar 2022 at 10:17, Simon SPDBA Greener <si...@spdba.com.au> wrote: > 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$; > > _______________________________________________ > postgis-users mailing list > postgis-users@lists.osgeo.org > https://lists.osgeo.org/mailman/listinfo/postgis-users >
_______________________________________________ postgis-users mailing list postgis-users@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/postgis-users