How to Fillet (round) Interior Angles of SF line intersections

Question

I am working with OSM data to create vector street maps. For the roads, I use line geometry provided by OSM and add a buffer to convert the line to geometry that looks like a road.

My question is related to geometry, not OSM, so I will use basic lines for simplicity.

library(dplyr)
library(ggplot2)
library(sf)

# two lines representing an intersection of roads
l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))

# union the two lines (so they "intersect" instead of overlap once buffered
l2 <- l1 %>% st_union()

ggplot()+
  geom_sf(data = st_buffer(l2, dist = 0.75))

This will yield something like this:

Ultimately, I am trying to render roads similarly to map APIs like Google, OSM, Bing, etc. These fillet (round) the inner angles like this:

I've searched through the st_ series methods in the sf package but haven't found a solution. The closest I have found is some arguments for st_buffer that control the shape of endcaps and the number of angles used for outer angles (but not the inner angles).

Is there a simple/practical solution for this, or am I better off getting used to un-rounded intersections?

Thanks

UPDATE:

Lovalery and mrhellmann present two solutions below.

Ultimately, I selected mrhellmann's response because it fit my particular need; however, I did compare and explore the merits of each which I share below. I leave this discussion here, because someone else may have a slightly different use case and the differences are useful to understand.

Let's define the different methods used:

l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))
l2 <- l1 %>% st_union()

l2.base <- l2 %>% st_buffer(dist = 0.75, endCapStyle = "ROUND") # OP
l2.neg <- l2.base %>% st_buffer(dist = -0.25) # mrhellmann (st_buffer)
l2.smth <- l2.base %>% smooth(method = "ksmooth", smoothness = 3, n= 50L)# Lovalery (smoothr)

Lovalery raises a good point, the negative buffer approach also smooths the ends/extremities.

Comparing the three methods shows that the negative buffer approach affects more than the endcaps, it shrinks the geometry by the distance used in st_buffer on all sides.

The l2.neg issue can be solved in a workaround by adding the "smoothing" value (value used in the negative buffer) to the original buffering distance (0.75 + 0.25).

l3 <- l1 %>% st_union %>% st_buffer(dist = 0.75+0.25, endCapStyle = "ROUND")
l3.neg <- st_buffer(l3, dist = -0.25)

With this adjustment, the methods have similar outcomes. In fact, l3 may be more similar to l2.base than l2.smth. At this point I might express a preference for the l3.neg method for the following reasons:

it uses the same package (not a big deal),
if appears to be slightly more true to the l2 geometry,
the distance value for st_buffer (using the same units as the map's crs) might be easier for users to understand than the smoothness setting for smoothr (a factor of the mean distance between vertices -- more on this below).

HOWEVER, Lovalery's point about the accuracy of line ends raised another good point. The default (ROUND) endCapStyle used by st_buffer extends the length of the line by the buffer distance setting.

If the plot limits crops all street ends, this is irrelevant. If the ends of streets are to be rendered in the plot, they should be the correct length.

A more appropriate buffer setting is endCapStyle = FLAT. Note that the smoothness setting for l2.smth changes from smoothness = 3 to smoothness = 0.2 when not fitting a rounded endcap.

The line length renders as intended for the l2 and l2.smth methods. The l2.neg and l3.neg preserve the 90 deg. angles on the ends (a "nice to have" but not critical feature), but shortens the line length by the negative distance used in the smoothing buffer. Point for smoothr.

Finally, to be nitpicky, the curve provided by smoothr is not symmetrical. The greater the smoothness value, and the greater the difference in vertices lengths, the more obvious the asymmetry is. Here is the smoothr method using a smoothness = 1.

Note the elongation of the curves on the longer vertices. This may not be noticeable at most scales.

As Lovalery mentioned, this comes down to what's best for the individual use case. Also, will the average viewer notice if the street is one metre shorter than OSM says, or that the curve is ever so slightly asymmetrical? Precision is a funny thing like that with graphics...perfect precision probably doesn't matter, but it's math and numbers so it feels like it should.

All code from above:

library(dplyr)
library(ggplot2)
library(sf)
library(smoothr)


l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))#,"LINESTRING(-5 7,5 7)"))
l2 <- l1 %>% st_union() 

#### endCapStyle = ROUND ####

l2.base <- l2 %>% st_buffer(dist = 0.75, endCapStyle = "ROUND")
l2.neg <- l2 %>% st_buffer(dist = -0.25)
l2.smth <- l2 %>% smooth(method = "ksmooth", smoothness = 3, n= 50L)

l3 <- l1 %>% st_union %>% st_buffer(dist = 0.75+0.25, endCapStyle = "ROUND")
l3.neg <- st_buffer(l3, dist = -0.25)

ggplot()+
  geom_sf(data = l2.base, colour = "darkgreen", fill = "palegreen")+
  geom_sf(data = l2.neg, colour = "blue", fill = NA)+
  geom_sf(data = l2.smth, colour = "red", fill = NA)+
  geom_sf(data = l3.neg, colour = "cyan", fill = NA)+
  ggtitle("Method Comparison", 
        subtitle = "EndCapStyle ROUND: l2 (green) v l2.smth (red), l2.neg (blue), l3.neg (cyan)") +
  scale_y_continuous(breaks = c(-1:6), limits = c(-1,6))+
  scale_x_continuous(breaks = c(-6:6), limits = c(-6,6))

#### endCapStyle = FLAT ####   

l2.base <- l1 %>% st_union() %>% st_buffer(dist = 0.75, endCapStyle = "FLAT")
l2.neg <- l2.base %>% st_buffer(dist = -0.25)
l2.smth <- l2.base %>% smooth(method = "ksmooth", smoothness = 0.2, n= 50L)

l3 <- l1 %>% st_union %>% st_buffer(dist = 0.75+0.25, endCapStyle = "FLAT")
l3.neg <- st_buffer(l3, dist = -0.25)

ggplot()+
  geom_sf(data = l2.base, colour = "darkgreen", fill = "palegreen")+
  geom_sf(data = l2.neg, colour = "blue", fill = NA)+
  geom_sf(data = l2.smth, colour = "red", fill = NA)+
  geom_sf(data = l3.neg, colour = "cyan", fill = NA)+
  ggtitle("Method Comparison", 
          subtitle = "EndCapStyle FLAT: l2 (green) v l2.smth (red), l2.neg (blue), l3.neg (cyan)") +
  scale_y_continuous(breaks = c(-1:6), limits = c(-1,6))+
  scale_x_continuous(breaks = c(-6:6), limits = c(-6,6))

#### illustrate asymmetry of smoothr method as l4.smth ####

l4.smth <- l2.base %>% smooth(method = "ksmooth", smoothness = 1, n= 50L)

ggplot()+
  geom_sf(data = l2.base, colour = "darkgreen", fill = "palegreen")+
  geom_sf(data = l4.smth, colour = "red", fill = NA)+
  geom_sf(data = l3.neg, colour = "blue", fill = NA)+
  geom_abline(slope = 1, intercept = 3, colour = "gray", linetype = "dashed")+
  ggtitle("Method Comparison", 
          subtitle = "EndCapStyle FLAT: l3.neg (blue), smoothr smoothness = 1 (red)") +
  coord_sf(xlim = c(0,5), ylim = c(3,5))

user_7547327 · Answer

You can buffer the lines and then negative buffer that result:

ggplot()+
  geom_sf(data = st_buffer(st_buffer(l2, dist = 1), dist = -0.5))

More verbose, less nested:

library(dplyr)
library(ggplot2)
library(sf)

# two lines representing an intersection of roads
l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))

# union the two lines (so they "intersect" instead of overlap once buffered
l2 <- l1 %>% st_union()

# buffer the lines with a positive dist argument
l2_buffered <- st_buffer(l2, dist = 1)

# un-buffer (negative buffer?) the lines for rounded inner corners
#  with a negative dist argument
l2_unbuffered <- st_buffer(l2_buffered, dist = -0.5)

Both together, the red de-buffered polygon is overplotted by the black except in the curved interior angles:

ggplot()+
  geom_sf(data = st_buffer(st_buffer(l2, dist = 1.5), dist = -0.75), fill = NA, col = 'red') + 
  geom_sf(data = st_buffer(l2, dist = .75), fill = NA, col = 'black')

Zoomed in on one corner:

ggplot()+
  geom_sf(data = st_buffer(st_buffer(l2, dist = 1.5), dist = -0.75), fill = NA, col = 'red') + 
  geom_sf(data = st_buffer(l2, dist = .75), fill = NA, col = 'black') + 
  coord_sf(xlim = c(-2,0), ylim = c(1,3))

user_13661061 · Answer

I suggest an alternative approach to the one proposed by @mrhellmann. So, you will be able to choose the one that suits you best in relation to your goal.

I suggest you to smooth your buffered polygon using the smoothr package. The only difference with @mrhellmann's solution is that the ends/extremities of the polygon are also smoothed so they don't exactly overlap the original polygon boundary at that point, as you'll see in the reprex below (though I'm not sure if this difference is important for your purposes)

Reprex

library(sf)
library(smoothr) 

# two lines representing an intersection of roads
l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))

# union the two lines (so they "intersect" instead of overlap once buffered
l2 <- l1 %>% st_union()

Smooth of your buffer

You can, of course, modify the algorithm and parameters to control the type and degree of smoothing to suit your needs, cf. ?smooth.

l2_smooth <- smooth(st_buffer(l2, dist = 0.75), 
                    method = "ksmooth", smoothness = 3, n = 50L)

Overview of the smoothed polygon

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l2_smooth)

General comparison of the smoothed polygon (red) with the original one (black)

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l2_smooth, col = "red", fill = NA) + 
  ggplot2::geom_sf(data = st_buffer(l2, dist = .75),  col = 'black', fill = NA)

Zoom in one corner

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l2_smooth, col = "red", fill = NA) + 
  ggplot2::geom_sf(data = st_buffer(l2, dist = .75),  col = 'black', fill = NA) +
  ggplot2::coord_sf(xlim = c(-1.2,-0.6), ylim = c(1.9,2.4))

Zoom in one extremity

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l2_smooth, col = "red", fill = NA) + 
  ggplot2::geom_sf(data = st_buffer(l2, dist = .75),  col = 'black', fill = NA) +
  ggplot2::coord_sf(xlim = c(-6,-4), ylim = c(2,4))

^{Created on 2021-10-16 by the reprex package (v2.0.1)}

EDIT/UPDATE

First of all, thank you very much @AWaddington for this in-depth report of the different proposed solutions: yours, @mrhellmann's answer and mine. It's always nice to have such constructive exchanges.

So, as you invited me to do, the purpose of this edit/update is to clarify and improve the approach I proposed and to give you my comparative analysis of the different methods.

I agree with your argument, and that of mrhellmann on the following points:

The solution proposed by mrhellmann does not indeed imply any additional package (which can be an advantage in some cases). On this point I can't "fight" and I gladly bow ;-)

The difference you point out concerning the end of the segments between the argument endCapStyle = "ROUND" and endCapStyle = "FLAT" does exist concerning the preservation (or not) of the exact length of the segments. In terms of cartographic representation "we are however in the thickness of the line" (except, of course, at a very large – and unlikely – scale). That said, let's be conservative and keep the actual length: it's always better if for some reason (geocomputation) you need the exact length to the nearest meter. Therefore in the reprex below, I have only used the argument endCapStyle = FLAT considering that this option is suitable for all cases.

However, it is possible to correct the smoothing method in such a way that points 2 and 3 that you raised become invalid:

Regarding point 2, as you will see in the reprex below, one can even bring the rounding of the internal angle even closer to the right angle than with the solution proposed by mrhellman. And it is easy to choose the shape of the rounding by varying the smoothness argument of the smooth function: the closer it is to 0, the closer the rounding is to the right angle. On the other hand, concerning the symmetry problem that you rightly raised, you just need to increase the number of nodes/points before smoothing. Indeed, the smoothing methods rely on these nodes/points and if they are not evenly/symmetrically distributed on each side of the angles, the smoothing itself cannot be symmetrical. To fix this, simply increase the number of points (in the reprex below, I have parametrized the function to insert a node/point every meter).

Finally, since it is possible to insert additional nodes/points in map units (here, in the reprex, every meter) the difference highlighted in point 3 no longer seems valid.

To conclude, the remaining differences between mrhellmann's approach and mine seem extremely tenuous and I think either solution may be appropriate.

mrhellman's approach:

Pros: no additional package, conservation of right angles at the ends of the segments.

Cons: less easy to adjust the rounding of the internal angle (to my opinion!)

lovalery' approach:

Pros: easier to adjust the rounding of the internal angle.

Cons: an additional package, slightly rounded corners at the ends of the segments (but not visible at a normal map scale)

Reprex

NB: only endCapStyle = "FLAT" option has been retained for this reprex (cf. explanation above)

Use of solutions l2.base and l3.neg for comparisons

library(sf)
library(smoothr)

l1 <- st_as_sfc(c("LINESTRING(0 0,0 5)","LINESTRING(-5 3,5 3)"))
l2 <- l1 %>% st_union()

l2.base <- l2 %>% st_buffer(dist = 0.75, endCapStyle = "FLAT") # OP FLAT version

l3 <- l1 %>% st_union %>% st_buffer(dist = 0.75+0.25, endCapStyle = "FLAT") # mrhellmann (st_buffer) modified, FLAT version
l3.neg <- st_buffer(l3, dist = -0.25)

Smooth of your buffer (i.e. l2.base) while respecting the symmetry of the internal angles = lovalery (smoothr) modified, FLAT version

# 1. Densify the number of nodes/points (here, one node every meters, cf. max_distance argument)
l3.smth.dens <- smooth(st_buffer(l2, dist = 0.75,endCapStyle = "FLAT"), method = "densify", max_distance = 1)

# 2. Smooth with kernel method (as in my first post). Test of two smoothness indices (i.e. 0.2 and 1)
l3.smth.dens1 <- smooth(l3.smth.dens, method = "ksmooth", smoothness = 0.2, n = 50L)

l3.smth.dens2 <- smooth(l3.smth.dens, method = "ksmooth", smoothness = 1, n = 50L)

General comparison of one of the two smoothed polygon (red, the one with smoothness index = 1) with the original one (black) and the mrhellman's one (blue)

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l3.smth.dens2, col = "red", fill = NA) + 
  ggplot2::geom_sf(data = l3.neg, col = "blue", fill = NA) +
  ggplot2::geom_sf(data = l2.base,  col = 'black', fill = NA) +
  ggplot2::geom_abline(intercept = 3, slope = 1, color="gray",linetype="dashed", size=1)

Zoom in one corner: comparison of the two smoothed polygons (red - smoothness index = 1 - and orange -smoothness index = 0.2) with the original one (black) and the mrhellman's one (blue)

Please, note that with a smoothness index = 0.2, the rounding of the inner angle is closer to the original geometry than the rounding of the mrhellman solution (if you want to have about the same rounding, set the smoothness index to ca. 0.6)

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l3.smth.dens2, col = "red", fill = NA) + 
  ggplot2::geom_sf(data = l3.smth.dens1, col = "orange", fill = NA) +
  ggplot2::geom_sf(data = l3.neg, col = "blue", fill = NA) +
  ggplot2::geom_sf(data = l2.base,  col = 'black', fill = NA) +
  ggplot2::geom_abline(intercept = 3, slope = 1, color="gray",linetype="dashed", size=1) +
  ggplot2::coord_sf(xlim = c(-1.55,-0.6), ylim = c(1.4,2.4))

Zoom in one extremity: comparison of the two smoothed polygons (red - smoothness index = 1 - and orange -smoothness index = 0.2) with the original one (black) and the mrhellman's one (blue)

Please, note that with a smoothness index = 0.2, the rounding of the ends is insignificant despite a disproportionately large rendering scale

ggplot2::ggplot()+
  ggplot2::geom_sf(data = l3.smth.dens2, col = "red", fill = NA) +
  ggplot2::geom_sf(data = l3.smth.dens1, col = "orange", fill = NA) +
  ggplot2::geom_sf(data = l3.neg, col = "blue", fill = NA) +
  ggplot2::geom_sf(data = l2.base,  col = 'black', fill = NA) +
  ggplot2::geom_abline(intercept = 3, slope = 1, color="gray",linetype="dashed", size=1) +
  ggplot2::coord_sf(xlim = c(-6,-4), ylim = c(2,4))

^{Created on 2021-10-17 by the reprex package (v2.0.1)}

How to Fillet (round) Interior Angles of SF line intersections

UPDATE:

2 Answers

EDIT/UPDATE