smoothr: spatial feature smoothing in R

For a project at work, one of my colleagues is generating polygons from raster data, which he then needs to smooth out to turn the sharp corners into smooth, natural looking curves. Although polygon smoothing seems like it should be a fairly commonly used GIS tool, we’ve been unable to find a good open source solution to this problem. ArcGIS has the Smooth Polygon tool that works nicely; however, given that smoothing is the final step in a large, automated, R-based workflow on Linux, it’s frustrating to have to use a commercial Windows program for the final step. Although I know almost nothing about smoothing algorithms, I did a little Googling and started putting together an R package, So, I introduce smoothr, now on GitHub and CRAN.

In this post, I’ll introduce the package with the hopes of stimulating other, more knowledgeable, folks to help test the package and implement some more advanced smoothing algorithms.

The final dataset that comes with this package, jagged_raster, is a simulated occurrence probability for a species, consisting of a spatially auto-correlated raster layer with values between 0 and 1. This raster can be used to experiment with smoothing polygons generated from rasters.

Smoothing methods

Thus far, I’ve implemented two simple smoothing methods: Chaikin’s corner cutting algorithm and spline interpolation. Both are accessed with the smooth() function, and all methods work on spatial lines and polygons in sf and sp format.

Chaikin’s corner cutting algorithm

Chaikin’s corner cutting algorithm smooths by iteratively replacing every point by two new points: one 1/4 of the way to the next point and one 1/4 of the way to the previous point. Consult the references below for details, but essentially the idea is to iteratively cut off corners until the curve is smooth. I’ve found this method to produce fairly natural looking smooth curves, although they’re a little more “boxy” than I’d like, and the algorithm has the benefit of only requiring a single parameter: the number of smoothing iterations.

This method can be applied with smooth(x, method = "chaikin"). Here’s what this looks like for the polygons:

Spline interpolation

This method applies a spline interpolation to the x and y coordinates independently using the built-in spline() function. For polygons (and closed lines), method = "periodic" is used to avoid getting a kink at the start/end of the curve defining the boundary. Unlike the corner cutting algorithm, this method results in a curve that passes through the vertices of the original curve, which may be a desirable feature. Unfortunately, this results in an unnaturally wiggly curve. Spline interpolation requires a parameter specifying the number of points to interpolate at, which can either be an absolute number or a relative increase in the number of vertices.

This method can be applied with smooth(x, method = "spline"). Here’s what this looks like for the polygons:

Raster-to-polygon conversion

The whole point of this smoothr business was to smooth out polygons generated from rasters, so let’s work through a quick example of that. Treating jagged_raster as the occurrence probability for a species, imagine we want to produce a range map for this species, showing where it occurs with at least 50% probability. We can convert the raster to a binary presence/absence map, then polygonize.

Not perfect, it still clearly looks like this range map came from a raster, but those slightly smoother corners are certainly easier on the eyes!

Future development

As it stands, smoothr has two decent methods for smoothing out sharp corners on polygons; however, my hope is to build out this package with a little more functionality:

Densification: add a method for densification, i.e. adding extra points to curves. While not really a smoothing algorithm, this could be a precursor to other smoothing algorithms, and it is another GIS tool that doesn’t yet exist in R.

Kernel smoothing: kernel smoothing, with the built-in function ksmooth() or the KernSmooth package, could be an alternative way of smoothing curves, likely in conjunction with densification.

Local regression: I have seen some suggestion that local regression (something I’m not even remotely familiar with) could be used for smoothing. There’s the built-in loess() function as well as the locfit package.

PAEK: ArcGIS refers to their smoothing algorithm as PAEK (Polynomial Approximation with Exponential Kernel), which is apparently described in the appendix of this paper. However, having read the appendix multiple times, I’ve found their description extremely vague and confusing, possibly intentionally so to keep their algorithm proprietary.

Geographically aware algorithms: none of the methods I’ve looked at are geographically aware, they just treat coordinates as Cartesian points in the plane. Is this an issue? If so, is there a better way? I don’t know, but hopefully there’s someone out there smarter than me that does.