[postgis-users] Cut extent box2d into small regular boxes, and iterate through each box of features
Shaozhong SHI
shishaozhong at gmail.com
Tue Mar 29 02:33:01 PDT 2022
Thanks. Simon,
I will have a close look and study.
Regards,
David
On Tue, 29 Mar 2022 at 10:17, Simon SPDBA Greener <simon at 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 at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20220329/cd651481/attachment.html>
More information about the postgis-users
mailing list