[postgis-users] Splitting large polygons into smaller polygons

Thu May 27 18:42:33 PDT 2010

Hi All,

I have a requirement to split large polygons into smaller pieces which
can't have any interior rings.

I started writing a plpgsql function for my solution with limited
spatial experience and less postgis experience and want some tips on how
it can be improved.

I pass the function two parameters, the geometry and the size I want the
smaller polygons to be. It then returns a setof geometries with the
results.

If the split of polygon has an rings in it, the function recursively
calls itself with pdp_size / 2 until all objects have no interior
rings..

The function works well with larger values in pdp_size (ie 1.0) but is
quite slow when using smaller values (ie 0.002).

To give you an idea of the data I am working with, it is Oceans around
the Australian mainland and Lakes/Rivers/Catchments in Australia.
Ocean's I split up using a size of 1.0 and Lakes/Rives/Catchments using
a size of 0.002.

My postgres/Postgis versions are

SELECT version(), postgis_full_version();

PostgreSQL 8.3.6 on x86_64-redhat-linux-gnu, compiled by GCC gcc (GCC)
4.1.2 20071124 (Red Hat 4.1.2-42),POSTGIS="1.3.5"
GEOS="2.2.3-CAPI-1.1.1" PROJ="Rel. 4.6.0, 21 Dec 2007" USE_STATS

Regards,

Josh

-- Function: principal_gis.denonut_primitive(geometry, double precision)

-- DROP FUNCTION principal_gis.denonut_primitive(geometry, double
precision);

CREATE OR REPLACE FUNCTION principal_gis.denonut_primitive(pgeom_obj
geometry, pdp_size double precision)

  RETURNS SETOF geometry AS

$BODY$

DECLARE

  vrec_ring record;

  -- Bounding box extents

  vdp_xmin double precision;

  vdp_xmax double precision;

  vdp_ymin double precision;

  vdp_ymax double precision;

  -- Cursors

  vdp_xpos double precision;

  vdp_ypos double precision;

  -- Stepping

  vdp_xstep double precision;

  vdp_ystep double precision;

  -- Variables to hold different versions of the geometry

  vgeom_orig geometry;

  vgeom_test geometry;

  vgeom_cut  geometry;

  vgeom_coll geometry;

  vgeom_dump geometry;

  -- Variable to hold SRID

  vi_srid integer;

BEGIN

  /******************************************************

    Do split

  ******************************************************/

  -- Get Bounding box extents

  vdp_xmin := ST_XMin(ST_Envelope(pgeom_obj));

  vdp_xmax := ST_XMax(ST_Envelope(pgeom_obj));

  vdp_ymin := ST_YMin(ST_Envelope(pgeom_obj));

  vdp_ymax := ST_YMax(ST_Envelope(pgeom_obj));

  vi_srid := ST_SRID(pgeom_obj);

  -- Setup xmin, xmax as xstep so they go east to west.

  IF vdp_xmax < 0 THEN

    vdp_xpos := vdp_xmax;

    vdp_xmax := vdp_xmin;

    vdp_xmin := vdp_xpos;

    vdp_xstep := 0 - pdp_size;

  ELSE

--    vdp_xpos := vdp_xmin;

    vdp_xstep := pdp_size;

  END IF;

  -- Setup ymin, ymax as ystep so they go north to south.

  IF vdp_ymax < 0 THEN

    vdp_ypos := vdp_ymax;

    vdp_ymax := vdp_ymin;

    vdp_ymin := vdp_ypos;

    vdp_ystep := 0 - pdp_size;

  ELSE 

--    vdp_ypos := vdp_ymin;

    vdp_ystep := pdp_size;

  END IF;

  -- Make Test Polgon for Intersection

  vgeom_orig := ST_GeomFromText('POLYGON((' || vdp_xmin            || '
' || vdp_ymin || ',' 

                                           || vdp_xmin + vdp_xstep || '
' || vdp_ymin || ',' 

                                           || vdp_xmin + vdp_xstep || '
' || vdp_ymin + vdp_ystep || ','

                                           || vdp_xmin             || '
' || vdp_ymin + vdp_ystep || ','

                                           || vdp_xmin             || '
' || vdp_ymin || '))', vi_srid);

  -- Check geometry 

  IF ST_Area(pgeom_obj) < (ST_Area(vgeom_orig) / 2) AND

     ST_Area(ST_Envelope(pgeom_obj)) < (ST_Area(ST_Envelope(vgeom_orig))
* 2) AND

     NumInteriorRings(pgeom_obj) = 0

   THEN

     -- Regardless of the shape, the object is small enough to pass.

     --   That is, the actual area of the geometry is less than half the
split area AND

     --   the bounding box of the geometry is not double the split area
AND

     --   the geometry has no inner rings

     RETURN NEXT pgeom_obj;

  ELSE

    -- While we have not gone past the original polygon (To the east)

    WHILE NOT pgeom_obj << vgeom_orig LOOP

      -- Initalise Test Polygon (Y axis)

      vgeom_test := vgeom_orig;

      --While we have not gone below the original polygon (To the south)

      WHILE NOT pgeom_obj |>> vgeom_test LOOP

        -- Get the intersection

        vgeom_cut := SetSRID(ST_Intersection(AsText(pgeom_obj),
SetSRID(vgeom_test,-1)), vi_srid);

/*

EXPLANATION OF ABOVE QUERY

 -- In our version of postgis (or more precisely GEOS), there is a bug
with precision of geometeries.

 -- The AsText function has some precision reduction which makes the
error go way. This bug has been

 -- fixed in later versions.

 -- Info at

 --
http://geos.refractions.net/pipermail/geos-devel/2005-May/001441.html

*/

        -- Did anything intersect?

        IF ST_IsEmpty(vgeom_cut) THEN

          -- No, do nothing

          NULL;

        ELSE

          -- Yes, check what turned out in the results

          IF GeometryType(vgeom_cut) = 'POLYGON' THEN

            -- Was there rings?

            IF NumInteriorRings(vgeom_cut) > 0 THEN

              -- Yes Rings, Dump each polygon through this function
again

              RETURN QUERY SELECT denonut_primitive FROM
principal_gis.denonut_primitive(vgeom_cut, pdp_size / 2);

            ELSE

              RETURN NEXT vgeom_cut;

            END IF;

          ELSE

            -- Not a polygon (linestring, point, multi*something* or
geometrycollection)

            FOR vgeom_dump IN SELECT geom FROM ST_Dump(vgeom_cut) LOOP

              IF ST_IsValid(vgeom_dump) THEN

                vgeom_coll := ST_BuildArea(vgeom_dump);

                IF NumInteriorRings(vgeom_coll) > 0 THEN

                  RETURN QUERY SELECT denonut_primitive FROM
principal_gis.denonut_primitive(vgeom_coll, pdp_size / 2);

                ELSE

                  RETURN NEXT vgeom_coll;

                END IF;

              ELSE

                -- Yes, skip it

                CONTINUE;

              END IF;

            END LOOP;

          END IF;

        END IF; /* ST_IsEmpty */

        -- Move the Test Polygon (Y Axis)

        vgeom_test := ST_Translate(vgeom_test, 0, vdp_ystep);

      END LOOP;

      -- Move the Test Polygon (X Axis)

      vgeom_orig := ST_Translate(vgeom_orig, vdp_xstep, 0);

    END LOOP;

  END IF;

  -- We are finished, return the result set now

  RETURN;

END;

$BODY$

  LANGUAGE 'plpgsql' STABLE STRICT

  COST 100

  ROWS 1000;

ALTER FUNCTION principal_gis.denonut_primitive(geometry, double
precision) OWNER TO postgres;

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20100528/e0ae1f2b/attachment.html>