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
