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;
_______________________________________________ postgis-users mailing list postgis-users@postgis.refractions.net http://postgis.refractions.net/mailman/listinfo/postgis-users