rgeos::gBuffer shrink without loss of points

Question

I need to be able to shrink a polygon of lat/lon data without loss of points; more the point, I need the points to be effectively "smooshed" in the correct direction. Typically, gBuffer works fine, but there is no assurance on the number of points nor the relative spacing of them. Ultimately, there are properties with each point that I need to preserve, and splines, smoothing, and other "nice efficiencies" with gBuffer and polygon growing/shrinking do not allow me to preserve those properties with sufficient confidence of a 1-to-1 mapping.

Example:

library(rgeos)   # gBuffer

dat <- structure(list(x = c(6, 5.98, 5.94, 5.86, 5.75, 5.62, 5.47, 5.31, 5.13, -4.87, -5.04, -5.22, -5.39, -5.55, -5.69, -5.81, -5.9, -5.96, -6, -6, -6, -5.96, -5.9, -5.81, -5.69, -5.55, -5.39, -5.22, -5.04, -3.04, -2.87, -2.69, -2.53, -2.38, -2.25, -2.14, -2.06, -2.02, -2, -2, -1.96, -1.9, -1.81, -1.69, -1.55, -1.39, -1.22, -1.04, -0.87, 1.13, 1.31, 1.47, 1.62, 1.75, 1.86, 1.94, 1.98, 2, 2, 2, 2.04, 2.1, 2.19, 2.31, 2.45, 2.61, 2.78, 2.96, 4.96, 5.13, 5.31, 5.47, 5.62, 5.75, 5.86, 5.94, 5.98, 6), y = c(5, 5.18, 5.35, 5.51, 5.66, 5.78, 5.88, 5.95, 5.99, 5.99, 6, 5.97, 5.92, 5.83, 5.72, 5.59, 5.43, 5.27, 5.09, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, -5, -1.91, -1.73, -1.57, -1.41, -1.28, -1.17, -1.08, -1.03, -1, -1.01, -1.01, -1.05, -1.12, -1.22, -1.34, -1.49, -1.65, -1.82, -2, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, 5)), row.names = c(NA, -78L), class = "data.frame")

# "shrink-wrap"
sp <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon( as.matrix(dat) )), "dat")))
sp2 <- gBuffer(sp, width = -0.5)
dat2 <- as.data.frame(sp2@polygons[[1]]@Polygons[[1]]@coords)

c(nrow(dat), nrow(dat2))
# [1] 78 97

Immediately we see the change in the number of points. I recognize that most of the time, this is a desired trait of gBuffer, so perhaps rgeos is not the best tool for this transformation.

library(ggplot2) # just for vis here
ggplot(dat, aes(x, y)) +
  geom_path() + geom_point() +
  geom_path(data = dat2, color = "red") + geom_point(data = dat2, color = "red")  

This image has the effect on the overall shape I want, but it has increased the number of points, which means that I can no longer rely on a 1-to-1 relationship to the original points.

In general, the polygons are not symmetric, and many have in-croppings like this where many methods of "pulling" points in a particular direction will be biased or in the wrong direction.

I can find no option in gBuffer nor other functions in rgeos to be able to preserve the number and basic spatial-relationship of the points. I don't need "perfect" shrinking, if that changes things, but it should not deviate significantly.

Answer 1

This might work if you're satisfied with how you're currently shrinking the polygons. It builds on that to get a 1:1 point mapping from the old (large) points to the new (smaller) polygon.

library(rgeos)   # gBuffer
library(sf)
library(tidyverse)
dat <- structure(list(x = c(6, 5.98, 5.94, 5.86, 5.75, 5.62, 5.47, 5.31, 5.13, -4.87, -5.04, -5.22, -5.39, -5.55, -5.69, -5.81, -5.9, -5.96, -6, -6, -6, -5.96, -5.9, -5.81, -5.69, -5.55, -5.39, -5.22, -5.04, -3.04, -2.87, -2.69, -2.53, -2.38, -2.25, -2.14, -2.06, -2.02, -2, -2, -1.96, -1.9, -1.81, -1.69, -1.55, -1.39, -1.22, -1.04, -0.87, 1.13, 1.31, 1.47, 1.62, 1.75, 1.86, 1.94, 1.98, 2, 2, 2, 2.04, 2.1, 2.19, 2.31, 2.45, 2.61, 2.78, 2.96, 4.96, 5.13, 5.31, 5.47, 5.62, 5.75, 5.86, 5.94, 5.98, 6), y = c(5, 5.18, 5.35, 5.51, 5.66, 5.78, 5.88, 5.95, 5.99, 5.99, 6, 5.97, 5.92, 5.83, 5.72, 5.59, 5.43, 5.27, 5.09, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, -5, -1.91, -1.73, -1.57, -1.41, -1.28, -1.17, -1.08, -1.03, -1, -1.01, -1.01, -1.05, -1.12, -1.22, -1.34, -1.49, -1.65, -1.82, -2, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, 5)), row.names = c(NA, -78L), class = "data.frame")

# "shrink-wrap"
sp <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon( as.matrix(dat) )), "dat")))
sp2 <- gBuffer(sp, width = -0.5)
dat2 <- as.data.frame(sp2@polygons[[1]]@Polygons[[1]]@coords)

## New methods begin here
# change objects to type `sf`
sp_sf <- st_as_sf(sp)
sp2_sf <- st_as_sf(sp2)

dat_sf <- dat %>% st_as_sf(coords = c('x', 'y'))
dat2_sf <- dat2 %>% st_as_sf(coords = c('x', 'y'))

# The plot so far, saved for building on further down
p <- ggplot() + 
  geom_sf(data = sp_sf, color = 'blue', fill = NA) + 
  geom_sf(data = dat_sf, color = 'blue') +
  geom_sf(data = sp2_sf, color = 'red', fill = NA) + 
  geom_sf(data = dat2_sf, color = 'red')




# Using st_nearest_points original points to new small polygon
#  results in (perpendicular?) lines from old points to new small polygon
near_lines <- st_nearest_points(dat_sf, sp2_sf)

#plotted together:
p + geom_sf(data = near_lines, color = 'black')

## Zooming in on a problem area
p + geom_sf(data = near_lines, color = 'black') + 
  coord_sf(xlim = c(-3, 0), ylim = c(-2,0))

# Get only 1:1 points for shrunken polygon
# a small buffer had to be added, as some points were not showing up
# you may need to adjust the buffer, depending on your data & projection
new_points <- st_intersection(st_buffer(near_lines, .001), sp2_sf)

# All together now:
p + geom_sf(data = near_lines, color = 'black') + 
  geom_sf(data = new_points, color = 'green', size = 4) +
  coord_sf(xlim = c(-3, 0), ylim = c(-2,0))

^{Created on 2020-12-20 by the reprex package (v0.3.0)}

Answer 2

I don't know of a way to do what you're asking using a package, but I've put together a little example that I think might help. This approach does assume a center around (0, 0).

The generic approach converts your cartesian points to polar, and then using some factor REDUCTION_FACTOR, you can scale the distance of your points from the origin. However, for the points that define the concave section of your polygon, you'll notice very little change. So what I've done is added a factor that shifts the points that are closer to the origin differently (according to some arbitrary cutoff CONVEX_SHIFT_CUTOFF that's a function of the extremes of the original polygon) by a factor defined by the REDUCTION_FACTOR. The CONVEX_SHIFT_CUTOFF could likely be set by considering the distribution of points of a given geometry.

I did have to fudge the multiples a bit, but I think you can likely adapt this for your dataset if the scale of your geometry asymmetry is not too great. There are probably also things you could do to orient each geometry in a similar fashion or ifelse conditions to properly account for differences in concavity but it's somewhat difficult to say without more information.

library(rgeos)   # gBuffer

dat <- structure(list(x = c(6, 5.98, 5.94, 5.86, 5.75, 5.62, 5.47, 5.31, 5.13, -4.87, -5.04, -5.22, -5.39, -5.55, -5.69, -5.81, -5.9, -5.96, -6, -6, -6, -5.96, -5.9, -5.81, -5.69, -5.55, -5.39, -5.22, -5.04, -3.04, -2.87, -2.69, -2.53, -2.38, -2.25, -2.14, -2.06, -2.02, -2, -2, -1.96, -1.9, -1.81, -1.69, -1.55, -1.39, -1.22, -1.04, -0.87, 1.13, 1.31, 1.47, 1.62, 1.75, 1.86, 1.94, 1.98, 2, 2, 2, 2.04, 2.1, 2.19, 2.31, 2.45, 2.61, 2.78, 2.96, 4.96, 5.13, 5.31, 5.47, 5.62, 5.75, 5.86, 5.94, 5.98, 6), y = c(5, 5.18, 5.35, 5.51, 5.66, 5.78, 5.88, 5.95, 5.99, 5.99, 6, 5.97, 5.92, 5.83, 5.72, 5.59, 5.43, 5.27, 5.09, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, -5, -1.91, -1.73, -1.57, -1.41, -1.28, -1.17, -1.08, -1.03, -1, -1.01, -1.01, -1.05, -1.12, -1.22, -1.34, -1.49, -1.65, -1.82, -2, -4.91, -5.09, -5.27, -5.43, -5.59, -5.72, -5.83, -5.92, -5.97, -6, -6, -5.99, -5.95, -5.88, -5.78, -5.66, -5.51, -5.35, -5.18, 5)), row.names = c(NA, -78L), class = "data.frame")

# "shrink-wrap"
sp <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon( as.matrix(dat) )), "dat")))
sp2 <- gBuffer(sp, width = -0.5)
dat2 <- as.data.frame(sp2@polygons[[1]]@Polygons[[1]]@coords)


REDUCTION_FACTOR <- 0.92
CONVEX_SHIFT_CUTOFF <- 0.55

x <- c()
y <- c()
max_x <- max(dat$x)
max_y <- max(dat$y)
for ( i in 1:length(dat$x)) {
  x_p <- dat$x[i]
  y_p <- dat$y[i]
  r <- sqrt(x_p^2 + y_p^2)
  theta <- atan2(y_p, x_p)
  r <- REDUCTION_FACTOR * r
  if(abs(x_p) < CONVEX_SHIFT_CUTOFF * max_x) {
    x <- c(x, r*cos(theta)*(1 + 4*(1-REDUCTION_FACTOR))  )
  } else {
    x <- c(x, r*cos(theta))
  }
  
  if(abs(y_p) < CONVEX_SHIFT_CUTOFF*max_y ) {
    y <- c(y, r*sin(theta)*(1 - 4*(1-REDUCTION_FACTOR) ))
  } else {
    y <- c(y, r*sin(theta)) 
  }
  
}

dat3 <- data.frame(x = x, y = y)

library(ggplot2) # just for vis here
ggplot(dat, aes(x, y)) +
  geom_path() + geom_point() +
  geom_path(data = dat2, color = "red") + geom_point(data = dat2, color = "red")  +
  geom_path(data = dat3, color = "blue") + geom_point(data = dat3, color = "blue")  

And this approach produces the following scaled result.