[postgis-users] Cleaning polygons, function documentation and a buffering problem (bug?) (longish)

Sat Nov 24 03:27:51 PST 2007

Greetings list,

Martin Davis/Kevin Neufeld provided a few weeks ago a technique to clean 
non-hole polygons from
topological errors.  In my PostGIS function list their technique now 
features as the following function:

CREATE OR REPLACE FUNCTION cleanpoly_nohole(geometry)
  RETURNS geometry AS
$BODY$
select st_buildarea(uniongeom) as geom
from ( select st_union(ringgeom, st_startpoint(ringgeom)) as uniongeom
       from ( select st_exteriorring($1::geometry) as ringgeom
             ) as ring
      ) as fixed_ring
$BODY$
  LANGUAGE 'sql' STABLE STRICT;

See their original contributions for further explanation.  This 
technique has also made it into the wiki.  There,
it is written that "this may not create the expected topology if the 
polygon contains holes."  I believe this should
read "this will not . . . "  So, it begs the question what to do with 
polygons that have holes.  I believe the following
works nicely, but I would appreciate informed comments.  The philosophy 
is to repair the outer polygon, and the
holes separately, and then determine their difference:

CREATE OR REPLACE FUNCTION cleanpoly(geometry)
  RETURNS geometry AS
$BODY$
select st_difference(
          cleanpoly_nohole($1::geometry),
          st_collect(
            cleanpoly_nohole(
              st_makepolygon(st_interiorringn($1,n)))))
from generate_series(1,st_numinteriorrings($1)) n ;
$BODY$
  LANGUAGE 'sql' STABLE STRICT;

If this works, the next challenge is cleaning multipolygons.

In general, I find it difficult to assess what some of the more advanced 
postgis functions actually do, and I find myself
experimenting a lot before making progress.  I would appreciate any 
references for deepening understanding of functions
beyond the online manual.  I have especially battled with collect, union 
and buffer recently.  But similar questions can be
posed for instance on polygonize. See one example problem below.

I would, possibly naively, expect that a buffer operation on a polygon 
with a positive second argument always returns
another polygon.  Not so.  See below real life example of a polygon that 
when buffered with value 15.0 gives us a
multipolygon.

"POLYGON((714.103 691.448,718.28 680.791,723.529 670.828,735.946
650.734,741.534 643.564,743.456 642.132,755.886 628.959,768.711
614.556,781.448 599.061,790.854 583.975,794.955 575.398,796.672
570.676,798.261 564.386,799.47 551.874,798.977 542.831,798.641
502.569,798.126 496.618,797.128 479.322,797.694 475.348,797.338
462.208,796.228 458.223,794.151 455.166,791.173 453.476,740.495
430.063,683.268 403.624,654.936 387.98,651.863 386.284,633.721
374.736,553.295 323.328,529.519 307.93,477.499 274.242,330.223
178.866,323.818 189.765,299.791 175.088,269.819 214.939,298.751
233.579,289.75 253.256,340.273 277.111,343.51 281.428,337.242
298.22,309.767 368.698,323.552 373.714,329.658 361.494,319.098
357.065,350.654 280.095,361.081 263.681,357.098 260.717,372.841
235.538,436.141 276.203,410.421 315.518,374.397 292.474,360.3
326.201,380.311 334.391,386.153 319.782,414.786 330.706,390.006
392.687,395.293 394.583,381.808 426.063,368.122 458.655,311.366
435.146,303.302 425.776,288.84 421.292,286.146 428.064,283.88
433.761,269.744 469.302,254.538 508.521,254.274 510.133,310.026
533.366,325.919 539.932,332.466 542.638,382.588 563.345,387.912
565.545,392.954 567.628,439.578 585.259,462.115 525.738,462.625
520.627,461.605 516.539,397.755 490.986,379.111 484.349,383.426
478.412,474.885 512.961,447.108 588.097,476.048 599.005,524.184
614.068,531.969 616.761,563.932 629.03,605.771 647.025,638.199
661.907,649.533 622.053,615.53 612.78,628.819 564.054,636.005
566.014,633.376 575.654,640.598 577.623,667.347 479.54,672.605
480.974,674.804 472.912,689.913 477.033,695.103 457.998,682.204
454.481,678.532 467.94,663.389 463.81,668.59 444.741,695.231
452.007,701.732 428.168,719.365 432.978,714.575 450.54,729.42
454.589,713.916 514.387,687.436 507.164,665.386 585.132,696.139
593.26,680.908 650.895,650.154 642.767,644.063 664.598,699.23
689.915,709.353 693.956,711.519 693.775,712.692 693.118,714.103
691.448),(608.479 579.933,580.072 572.982,600.562 489.249,654.485
502.444,644.038 545.137,625.727 540.657,615.685 581.696,608.479
579.933),(556.661 510.789,584.304 443.146,604.636 448.842,611.27
450.7,599.405 485.234,579.45 550.524,544.394 540.273,522.945
531.348,525.518 525.164,508.635 518.141,574.869 358.435,578.787
360.065,584.617 346.569,593.37 350.211,577.827 387.567,580.973
388.876,534.147 501.422,556.661 510.789),(301.911 505.413,297.714
510.076,268.586 498.749,283.823 461.743,279.374 459.912,289.081
436.335,293.938 438.335,292.329 442.242,327.382 456.675,317.295
481.174,321.988 483.106,301.911 505.413))"

What we get is a multipolygon with three nested shells.  Is this a bug 
of the
buffer function?  I am using: PostGIS version "1.3 USE_GEOS=1
USE_PROJ=1 USE_STATS=1" and GEOS "3.0.0rc4-CAPI-1.3.3"

regards,

Rolf de By