'Postgresql: Calculate the number of polygon islands or overlapping groups'


This post is about some work I did on how to calculate the number of number of groups of intersecting polygons in an array of polygons in Postgresql


First lets be a little more strict about what the definition of an island is. For the purposes of this algorithm an island consists of a set polygons where at least one polygon in the set overlaps another in the set.


I created a function to calculate the number of islands attached is the plpgsql code for the function.


CREATE OR REPLACE FUNCTION getNumIslands(polys polygon[]) RETURNS INT AS
\(BODY\)
DECLARE
islandList polygon[]; -- holds the set of obstacles that make up the current island.
islandCount INT; -- holds the count of the number of islands.
poly polygon; -- holds the polygon to start an island search
polys_copy polygon[]; -- A copy of the input list
inddex INT; -- keeps track of the current index in the polys_copy search

BEGIN
   
islandCount := 0;

islandList := polys;
polys_copy := polys;
inddex := 1;

WHILE array_length(polys_copy,1) > 0
LOOP
    -- RAISE NOTICE 'Poly_list length is: %', array_length(polys_copy, 1);
    -- grab the first poly
    poly := polys_copy[1];
    islandList := array_append(ARRAY[]::polygon[], poly);
    -- RAISE NOTICE 'Current island is %', islandList;
    -- remove the first poly because we dont care about it overlaping itself.
    polys_copy := polys_copy[2:array_length(polys_copy, 1)];
   
    -- Make a list of all the polys left that overlap with poly
    WHILE inddex <= array_length(polys_copy, 1)
    LOOP
        -- If the poly_copy overlaps the poly then add its index to the list
        -- of polys to be removed
        FOREACH poly IN ARRAY islandList
        LOOP
            IF polys_copy[inddex] && poly
            THEN
                islandList := array_append(islandList, polys_copy[inddex]);
                --
                IF array_length(polys_copy,1) = 1
                THEN
                    polys_copy := null;
                ELSE
                    polys_copy := polys_copy[1:inddex-1] || polys_copy[inddex+1:array_length(polys_copy, 1)];
                END IF;
                -- RAISE NOTICE 'Intersects with poly index: %', inddex;
                -- When a new intersecting poly is found it needs to be compared
                -- to every poly not just the polys left in the list
                inddex := 0;
                EXIT;
            END IF;
        END LOOP;
        -- RAISE NOTICE 'Current island is %', islandList;
        inddex := inddex + 1;
    END LOOP;
   
    islandCount := islandCount + 1;
    inddex := 1;
END LOOP;

RETURN islandCount;
               
END
               
\(BODY\) LANGUAGE plpgsql;
© 2013 Glen Berseth | www.fracturedplane.com