GEOSBuffer on a valid Polygon creates MultiPolygon with sliver geometry

[DahnJ]

This was originally reported by @bocalml in geopandas/geopandas#3330

Problem statement

For a valid Polygon, GEOSBuffer creates a MultiPolygon that contains the original buffered geometry and a small sliver polygon.

I reproduced it first in PyGEOS and then GEOS, posting both examples:

PyGEOS

import pygeos

# Original geometry in EPSG:32651
p = pygeos.polygons([
 (217166.2392768794, 2965585.426031656),
 (217157.0169161774, 2965590.9332145103),
 (217152.36481728853, 2965585.858676655),
 (217148.8168604683, 2965582.0901328623),
 (217142.12103207916, 2965581.3085645237),
 (217135.51007477118, 2965590.095169524),
 (217131.64161110792, 2965596.0147710186),
 (217135.46261837851, 2965605.754331463),
 (217141.33031409973, 2965611.3339293986),
 (217146.9890091528, 2965608.0245872694),
 (217151.19246039074, 2965613.0648114844),
 (217154.5256337439, 2965620.1316876863),
 (217162.73226182582, 2965620.823985329),
 (217167.91836137482, 2965624.4116914785),
 (217179.40574377665, 2965625.3853994287),
 (217182.00630587857, 2965615.7349281944),
 (217180.09617763676, 2965608.658363588),
 (217181.07185720128, 2965602.448688585),
 (217177.533621619, 2965593.3348772344),
 (217178.31749983027, 2965588.3382303184),
 (217173.1534413562, 2965581.290163225),
 (217166.2392768794, 2965585.426031656)
])
buffered_polygon = pygeos.buffer(p, 100_000)

print("Original Polygon WKT:", pygeos.to_wkt(p))
print("Buffered Polygon WKT:", pygeos.to_wkt(buffered_polygon))

# Output
Original Polygon WKT: POLYGON ((217166.239277 2965585.426032, ...
Buffered Polygon WKT: MULTIPOLYGON (((165831.979558 2879767.049334,...

GEOS

#include <stdio.h>  
#include <stdarg.h> 

#include <geos_c.h>

static void
geos_msg_handler(const char* fmt, ...)
{
    va_list ap;
    va_start(ap, fmt);
    vprintf (fmt, ap);
    va_end(ap);
}

int main()
{
    initGEOS(geos_msg_handler, geos_msg_handler);

    GEOSWKTReader* reader = GEOSWKTReader_create();
    GEOSGeometry* geom = GEOSWKTReader_read(
        reader, "\
        POLYGON ((\
            217166.239277 2965585.426032,\
            217157.016916 2965590.933215,\
            217152.364817 2965585.858677,\
            217148.81686 2965582.090133,\
            217142.121032 2965581.308565,\
            217135.510075 2965590.09517,\
            217131.641611 2965596.014771,\
            217135.462618 2965605.754331,\
            217141.330314 2965611.333929,\
            217146.989009 2965608.024587,\
            217151.19246 2965613.064811,\
            217154.525634 2965620.131688,\
            217162.732262 2965620.823985,\
            217167.918361 2965624.411691,\
            217179.405744 2965625.385399,\
            217182.006306 2965615.734928,\
            217180.096178 2965608.658364,\
            217181.071857 2965602.448689,\
            217177.533622 2965593.334877,\
            217178.3175 2965588.33823,\
            217173.153441 2965581.290163,\
            217166.239277 2965585.426032))"
    );

    GEOSGeometry* buffered = GEOSBuffer(geom, 100000, 30);

    GEOSWKTWriter* writer = GEOSWKTWriter_create();
    char* wkt = GEOSWKTWriter_write(writer, buffered);
    printf("Geometry: %s\n", wkt);

    /* Clean up allocated objects */
    GEOSWKTReader_destroy(reader);
    GEOSWKTWriter_destroy(writer);
    GEOSGeom_destroy(geom);
    GEOSGeom_destroy(buffered);
    GEOSFree(wkt);

    /* Clean up the global context */
    finishGEOS();
    return 0;
}

Returns

Geometry: MULTIPOLYGON (((165831.97174934018 2879767.0540048014 ...

Discussion

With sufficiently high quadsegs this goes away. The combination of large radius and low quadsegs is likely the culprit here.

>>> buffered_polygon = pygeos.buffer(p, 100_000, quadsegs=100)

Original Polygon WKT: POLYGON ((217166.239277 2965585.426032....
Buffered Polygon WKT: POLYGON ((208819.005842 2865931.168169, 208250.19331...

Environment

GEOS 3.12.1

[mwtoews]

Here is a simple reproducer:

geosop -a "POLYGON ((217166.239 2965585.426, 217157.017 2965590.933, 217152.365 2965585.859, 217148.817 2965582.09, 217142.121 2965581.309, 217135.51 2965590.095, 217131.642 2965596.015, 217135.463 2965605.754, 217141.33 2965611.334, 217146.989 2965608.025, 217151.192 2965613.065, 217154.526 2965620.132, 217162.732 2965620.824, 217167.918 2965624.412, 217179.406 2965625.385, 217182.006 2965615.735, 217180.096 2965608.658, 217181.072 2965602.449, 217177.534 2965593.335, 217178.317 2965588.338, 217173.153 2965581.29, 217166.239 2965585.426))" buffer 50000

This does not appear to be an issue with JTS, but here is what the input looks like:

and the MULTIPOLYGON output from geosop:

The bottom of the output geometry has a denser pattern of vertices, which is unusual.

A good workaround for such a large buffer is to simplify the input geometry.