Buffer intersects line it should buffer?

[zoeschindler]

I think I found abug in GEOS, but I am not 100% sure. I use the GEOS library via the R programming language using the terra, sf or geos packages which all rely on GEOS. With all of these three, I had issues with buffers around some lines. Starting from a certain buffer size, the buffer intersects the line it was supposed to buffer. These issues don't arise when the buffer is smaller, e.g. 1m instead of 2m.

I already opend two other issues for this in the respective packages: https://github.com/r-spatial/sf/issues/1986, https://github.com/rspatial/terra/issues/765. Someone suggested the issue would be at the GEOS package, so I wanted to let someone know although I use GEOS only via the R API.

Here is some example code and output in R:

# input data
ch_crd <- structure(
  c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3.50850510439515, 3.34156954142671, 
    2.87898194958676, 2.72374630132234, 2.2922345396554, 2.29312419890067, 
    2.32261681738335, 2.49720251348322, 3.19886481773371, 3.50850510439515, 
    1.62424111363924, 1.64507013529592, 1.71055495856562, 1.8867890234822, 
    2.38801074041958, 2.91005134572971, 3.22939872710867, 3.25684154033292, 
    2.97375702861173, 1.62424111363924), dim = c(10L, 3L), dimnames = list(
      NULL, c("id", "x", "y"))) 

# printed data frame with coordinates
> ch_crd
      id        x        y
 [1,]  1 3.508505 1.624241
 [2,]  1 3.341570 1.645070
 [3,]  1 2.878982 1.710555
 [4,]  1 2.723746 1.886789
 [5,]  1 2.292235 2.388011
 [6,]  1 2.293124 2.910051
 [7,]  1 2.322617 3.229399
 [8,]  1 2.497203 3.256842
 [9,]  1 3.198865 2.973757
[10,]  1 3.508505 1.624241

library(terra)

# create line object
ch_lin <- vect(ch_crd, "lines", crs = CRS("EPSG:25832")@projargs)

# create 2m buffer around object
ch_buff <- buffer(ch_lin, 2)

# plot buffer per section
plot(ch_buff, col = "lightblue")
plot(ch_lin, add = T)

library(sf)

# create line & buffer
ch_lin = st_linestring(ch_crd[, 2:3])
ch_lin = st_sfc(ch_lin, crs = "epsg:25832")
ch_buff = st_buffer(ch_lin, 2)

# plot line & buffer
plot(ch_buff, col = "lightblue")
plot(ch_lin, add = T)

[pramsey]

I'll wait for @dr-jts to weigh in on this, but to my mind this seems more like "do something weird, get a weird answer" than "answer is wrong". You're taking a linear loop and buffering it much further than its radius so the buffer lines on one side extend way past those on the other. Maybe the right answer is supposed to be "no hole please", it's certainly an odd hole, but I can see how it gets there. Does it go away if the start/end points are not the same?

Anyways, if you don't want the hole, converting the coordinates to a polygon instead of a closed line will surely give you the same shape but without the hole.

[zoeschindler]

Yeah, I mean if you say I accidentally misused the functions, then I can understand that. I just thought it was weird that it did work for some similar polygons but not for all. I actually want the hole in the middle, but either the correct one or the information that there is no hole. I think I will change my strategy and create negative buffers around polygons.

Changing the orders of the points (ch_lin <- vect(ch_crd[c(4:10, 2:4),], "lines", crs = CRS("EPSG:25832")@projargs)) does not change anything. Removing some points in-between changes the buffer but it's still intersecting the line. After opening the line by removing a point, no hole is created.

[pramsey]

The buffer and offset routines definitely have some historical issues with closed lines, and you've just exposed another one. If/when the buffering gets reworked, this may be a useful example. One way to get deterministic results on line buffers is to break them up into N little two-point lines, buffer all those, then union them together.

[changebetter1]

when i use geos in c++ project,same problem is happened. (geos.3.9.3 , geos3.11.0) double dDegree = dDistance / 100000.0;//NORTH_CHINA_METER_PER_DEGREE; auto pBuffer = pGeom->buffer(dDegree, iQuadrantSegments, iEndCapStyle); if dDistance is 1.0, pBuffer is right，but when dDistance is 100.0, pBuffer is error.

LINESTRING(113.612257740769 34.736520258499525,113.612239571148 34.736520283305,113.61223403854402 34.73651953299803)

[pramsey]

Add a picture, your explanation isn’t clear.

[changebetter1]

sorry, more infomations std::string str = pGeom->buffer(dDegree, iQuadrantSegments, iEndCapStyle)->toString();

pGeom is "LINESTRING(113.612257740769 34.736520258499525,113.612239571148 34.736520283305,113.61223403854402 34.73651953299803)"

if dDegree=0.00001 && iQuadrantSegments=8 && iEndCapStyle=2 the buffer is "POLYGON((113.61224023934092 34.73651028238346,113.61223538239776 34.73650962370657,113.61223269469023 34.73652944228945,113.61223822729423 34.73653019259646,113.61223958480016 34.73653028329568,113.61225775442118 34.73653025849016,113.61225772711683 34.7365102585088,113.61224023934092 34.73651028238346))".

if dDegree=0.0001 && iQuadrantSegments=8 && iEndCapStyle=2 the buffer is "POLYGON((113.61223943462633 34.736420283398175,113.61223943631178 34.73642151796493,113.61224726001439 34.73642204069104,113.61222060000638 34.73661862591251,113.61222613261037 34.73661937621951,113.61223970766967 34.73662028321179,113.61225787729067 34.736620258406305,113.61225760424732 34.73642025859267,113.61223943462633 34.736420283398175))"

if dDegree=0.0002 && iQuadrantSegments=8 && iEndCapStyle=2 the buffer is "MULTIPOLYGON(((113.61223929810465 34.73632028349138,113.61223930147555 34.73632275262487,113.6122574727879 34.73632396670726,113.61225746772564 34.736320258685865,113.61223929810465 34.73632028349138)),((113.61225750356027 34.73634650696535,113.61220716146876 34.73671771882704,113.61221269407275 34.73671846913399,113.61223984419135 34.736720283118615,113.61225801381234 34.73672025831311,113.61225750356027 34.73634650696535)),((113.61226053264417 34.73632417114587,113.61226611640751 34.73632454421457,113.61226644822325 34.73632209747596,113.61226091561927 34.73632134716898,113.61226053264417 34.73632417114587)))"

[changebetter1]

Maybe I found the reason when i try to use the function with QGis.

[pramsey]

I think this has to wait on a fix in https://github.com/locationtech/jts/issues/876