Search code examples
rlinepolygonterra

Unordered lines to closed polygon

Is there a way to convert a set of unordered lines in R (e.g. using terra) into a polygon?

Here is an example:

library(terra)

geom <- as.matrix(data.frame(
  "id" = rep(1:6,each = 2),
  "x" = c(1,1,3,3,4,3,3,1,4,3,1,3),
  "y" = c(2,3,1,2,2,2,4,3,2,4,2,1)
))

l <- vect(geom, type = "lines")
plot(l)

It does not need to be spatial polygon created in some package, I just want the coordinates to be arranged in a way that if I connect each point to the next point I get the polygon.

Here the required output for my example:

coords <- data.frame(
  "x" = c(1,1,3,4,3,3,1),
  "y" = c(2,3,4,2,2,1,2))

plot(coords)
lines(coords)

Edit: Thanks to the great answers, this is my solution when not using any additional packages.

geom$x <- geom$x - mean(geom$x)
geom$y <- geom$y - mean(geom$y)
geom$angle <- atan2(geom$x, geom$y)
geom <- unique(geom[order(geom$angle),])
geom <- rbind(geom, geom[1,])
Solution

  • Ok, thanks for the edit - let's start all over then. Now I get why you asked about the coordinates being re-sorted. Plotting your polygon with distinct coordinates results in borders crossing each other and does not meet your expected result, got it. In order to "connect each point to the next point" - usually one would do this (anti-)clockwise - let's try to work with some angles (adjusted from RodrigoAgronomia/PAR):

    library(terra)
#> terra 1.5.34
library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

# load data
geom <- data.frame(
  "id" = rep(1:6, each = 2),
  "x" = c(1, 1, 3, 3, 4, 3, 3, 1, 4, 3, 1, 3),
  "y" = c(2, 3, 1, 2, 2, 2, 4, 3, 2, 4, 2, 1)
) |> st_as_sf(coords = c("x", "y"))

# define a little helper function to calculate angles
calc_angle_to_mean <- function(sf) {
  
  x <- st_coordinates(sf) |> as.data.frame()
  y <- st_coordinates(sf) |> colMeans()
  
  diff <- x - y
  
  mapply(atan2, diff[["X"]], diff[["Y"]]) * 180 / pi
}

geom[["angle_deg"]] <- calc_angle_to_mean(geom)

# inspect result
geom
#> Simple feature collection with 12 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 4
#> CRS:           NA
#> First 10 features:
#>    id    geometry  angle_deg
#> 1   1 POINT (1 2) -108.43495
#> 2   1 POINT (1 3)  -63.43495
#> 3   2 POINT (3 1)  161.56505
#> 4   2 POINT (3 2)  116.56505
#> 5   3 POINT (4 2)  108.43495
#> 6   3 POINT (3 2)  116.56505
#> 7   4 POINT (3 4)   18.43495
#> 8   4 POINT (1 3)  -63.43495
#> 9   5 POINT (4 2)  108.43495
#> 10  5 POINT (3 4)   21.80141

# sort by angle, extract coordinates, keep only distinct ones
p <- geom |> 
  dplyr::arrange(angle_deg) |> 
  st_coordinates() |> 
  tibble::as_tibble() |> 
  dplyr::distinct()

# close the polygon
p <- rbind(p, p[1,])

# create SpatVector polygon
p <- as.matrix(p) |> vect(type = "polygons")

# inspect again
p
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 1, 0  (geometries, attributes)
#>  extent      : 1, 4, 1, 4  (xmin, xmax, ymin, ymax)
#>  coord. ref. :

plot(p)