rr-sfterra

Why does terra in R find the buffer of a polygon to not fully contain the buffer of the centroid?

I am using the terra package in R to get a buffer around polygons. However, I have noticed that terra::buffer() is not behaving as expected. In particular, for some polygons and buffer widths the buffered polygon does not fully contain a buffer with the same width of its centroid. Is there a bug or am I doing something wrong? I provide reproducible code below and a picture of the nonintuitive area difference. Thanks for your help!

Buffer Differences

library(terra)

# Define the vertices of the triangle (right-angled at point (0, 0))
vertices <- rbind(
  c(-93.746, 46.99),     # Vertex 1 (right angle)
  c(-93.736, 46.99),    # Vertex 2 (height of 10 units)
  c(-93.746, 46.9886),     # Vertex 3 (base of 5 units)
  c(-93.746, 46.99)      # Close the polygon by returning to Vertex 1
)

# Create the triangle polygon - EPSG:4269 is lat/lon
triangle_polygon <- vect(vertices, type = "polygons", crs = "EPSG:4269")

# Plot the rectangular triangle
plot(triangle_polygon)

# Get centroids
centroid <- terra::centroids(triangle_polygon)

# Buffer centroid (5000m)
centroid_buffer <- terra::buffer(centroid, 5000)

# Buffer polygon (5000m)
buffer <- terra::buffer(triangle_polygon, 5000)

# Do buffered centroid minus buffered polygon
diff <- terra::erase(centroid_buffer, buffer)

# The picture shows area contained in the buffered centroid, but not in the buffered polygon
plot(diff)

Posted as an issue in the terra github: Here.

Solution

  • It has to do with the way terra appears to be handling geographic coordinate system (GCS) spheroid coordinates, and the fact that some functions are not optimised for GCS data. As such, it is better to project your data, run your functions, then project back to the original GCS.

    Out of the box, the sf package is arguably easier in this use case. Then, if you need your data as a SpatVector or as a SpatialPolygons object, it is straightforward to convert them with terra::vect() or sf::as_Spatial() respectively.

    First, an example of why GCS data are problematic in this use case:

    library(terra)
library(sf)
library(ggplot2)

# Define polygon vertices
vertices <- rbind(
  c(-93.746, 46.99),
  c(-93.736, 46.99),
  c(-93.746, 46.9886),
  c(-93.746, 46.99)
)

# Create polygon sf with NAD83/EPSG:4269 CRS
triangle_polygon <- st_sfc(st_polygon(list(vertices)), crs = 4269) |>
  st_as_sf()

# Get polygon centroid
centroid <- st_centroid(triangle_polygon)

# Buffer centroid and polygon
centroid_buffer <- st_buffer(centroid, 5000) |>
  st_make_valid()

buffer <- st_buffer(triangle_polygon, 5000) |>
  st_make_valid()

# Return difference
diff <- st_difference(buffer, centroid_buffer)

ggplot() +
  geom_sf(data = diff) +
  theme(panel.background = element_blank())

    1

    Note the polygons have jagged edges. This is a result of applying the functions to GCS data. The st_buffer() function does have a nQuadSegs parameter, which increases the number of segments used to form a buffer, but this had no effect.

    Further, it was necessary to use st_make_valid(), otherwise an error was returned (this was specific to the 5000m buffer size and did not occur when a lower value was set).

    Now see the difference when using projected data. This example uses the UTM Zone CRS appropriate for your example data. I have provided a function to calculate the correct UTM Zone EPSG code. If your data spans more than one UTM Zone, and/or you want to use an alternative CRS, search online for a more appropriate EPSG code.

    # Get coordinates, convert to WGS84/EPSG:4326 for calculating UTM Zone EPSG
xy <- st_transform(triangle_polygon, 4326) |>
  st_coordinates()

# Function to return UTM Zone
utm_epsg <- function(lon, lat) {
  # Calculate UTM zone number
  utm_zone <- floor((lon + 180) / 6) + 1
  # Determine if Northern or Southern Hemisphere and assign EPSG code
  if (lat >= 0) {
    epsg <- 32600 + utm_zone  # Northern Hemisphere EPSG
  } else {
    epsg <- 32700 + utm_zone  # Southern Hemisphere EPSG
  }
  
}

# Get UTM Zone EPSG code
unique(mapply(utm_epsg, xy[,"X"], xy[,"Y"]))
# [1] 32615

# Project to UTM Zone
triangle_32615 <- st_transform(triangle_polygon, 32615)

# Get centroid
centroid <- st_centroid(triangle_32615)

# Buffer centroid and polygon
centroid_buffer <- st_buffer(centroid, 5000)

buffer <- st_buffer(triangle_32615, 5000)

# Return difference, project back to NAD83/EPSG:4269
diff <- st_difference(buffer, centroid_buffer) |>
  st_transform(4269)

# Convert to SpatVector or SpatialPolygons object
sv_diff <- vect(diff)
sp_diff <- as_Spatial(diff)

ggplot() +
  geom_sf(data = diff) +
  theme(panel.background = element_blank())

    2