prepr: Automatic Repair of Spatial Polygons

Introduction

The prepr R R package is used to repair broken spatial polygons. To achieve this, the package uses the constrained triangulation approach implemented in the prepair C++ library (Ledoux et al. 2014). Specifically, the main st_prepair() function is used to repair polygons.

Tutorial

Here we provide a short tutorial showing how the prepr package can be used to repair broken polygon geometries. To start with, we will load the package. We will also load the sf package for working with spatial data.

# load packages
library(sf)    # package for working with spatial data
library(prepr) # package for repairing spatial geometries

We will now show how the prepr package can be used to repair broken geometries. To achieve this, we will examine a range of different mechanisms by which polygons can be broken. For each of these mechanisms, we will manually create a broken polygon geometry based on the mechanism, repair the geometry (using st_prepair()), and then compare the repaired geometry with the original broken geometry. We will also confirm that the broken geometries have indeed been repaired (using sf::st_is_valid()).

A ‘bowtie’ polygon

# create broken polygon
x <- st_as_sfc("POLYGON((0 0, 0 10, 10 0, 10 10, 0 0))")

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Self-intersection[5 5]"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

Square with wrong orientation

# create broken polygon
x <- st_as_sfc("POLYGON((0 0, 0 10, 10 10, 10 0, 0 0))")

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Valid Geometry"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

# create broken polygon
x <- st_as_sfc(
  "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 2, 5 7, 10 7, 10 2, 5 2))"
)

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Self-intersection[10 2]"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

Dangling edge

# create broken polygon
x <- st_as_sfc("POLYGON((0 0, 10 0, 15 5, 10 0, 10 10, 0 10, 0 0))")

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Self-intersection[15 5]"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

Two adjacent inner rings

# create broken polygon
x <- st_as_sfc(
  paste(
    "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),",
    "(1 1, 1 8, 3 8, 3 1, 1 1), (3 1, 3 8, 5 8, 5 1, 3 1))"
  )
)

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Self-intersection[3 8]"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

Polygon with an inner ring inside another inner ring

# create broken polygon
x <- st_as_sfc(
  paste(
    "POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),",
    "(2 8, 5 8, 5 2, 2 2, 2 8), (3 3, 4 3, 3 4, 3 3))"
  )
)

# verify that polygon is broken, and display reason why it's broken
st_is_valid(x, reason = TRUE)

## [1] "Holes are nested[3 3]"

# repair the polygon
y <- st_prepair(x)

# verify that repaired polygon is NOT broken
st_is_valid(y)

## [1] TRUE

# visualize the polygons
par(mfrow = c(1, 2), mar = c(0, 0, 2.5, 0))
plot(x, main = "original", col = "#999999")
points(st_coordinates(st_cast(x, "POINT")), col = "#e41a1c", pch = 16, cex = 2)
plot(y, main = "repaired", col = "#999999")
points(st_coordinates(st_cast(y, "POINT")), col = "#e41a1c", pch = 16, cex = 2)

Citation

We recommend citing the methodology that underpins the prepr R package. To cite this work, please use:

Ledoux H, Arroyo Ohori K, and Meijers M (2014) A triangulation-based approach to automatically repair GIS polygons. Computers & Geosciences 66:121–131.

Getting help

If you have any questions about using the prepr R package or suggestions for improving it, please file an issue at the package’s online code repository.