Unordered lines to closed polygon

Question

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,])

Does this question help? https://stackoverflow.com/questions/47147242/convert-spatial-lines-to-spatial-polygons — Peter, Oct 13 '22 at 12:55

Chris · Answer 1 · 2022-10-15T21:09:52.880

If you are agnostic about the shape of the resulting polygon

library(terra)
coords_spV <- vect(coords, geom = c('x','y'))
plot(convHull(coords_spV))

Presented as you may be with unordered coordinates, and in search of a concave polygon for further purposes, you may choose to assemble and test resulting polygons for their characteristics, 'complex', 'convex', 'concave', depending on the order of assembly. Taking your geom and a borrowed function perm(), using only terra:

# a borrowed permutation function 
perm <- function(v) {
  n <- length(v)
  if (n == 1) v
  else {
    X <- NULL
    for (i in 1:n) X <- rbind(X, cbind(v[i], perm(v[-i])))
    X
  }
}
factorial(6)
[1] 720 # how many possible polygons
new_order <- perm(1:6)
c(rbind(new_order[1, ], new_order[1, ]))
 [1] 1 1 2 2 3 3 4 4 5 5 6 6 #original order
c(rbind(new_order[720, ], new_order[720, ]))
 [1] 6 6 5 5 4 4 3 3 2 2 1 1 

geoms <- list()
for (k in 1:nrow(new_order)) {
geoms[[k]] <- matrix(c(sample(geom[,2], 12, replace = TRUE), 
 sample(geom[,3], 12, replace = TRUE))[order(c(rbind(new_order[k, ],
 new_order[k, ])))], nrow=12,ncol=2, byrow = TRUE)
}
geoms_v <- list()
for (j in 1:length(geoms)) {
    geoms_v[[j]] <- vect(geoms[[j]], type = 'polygons')
}

geoms_v_valid <- list()
for (i in 1:length(geoms_v)) {
 geoms_v_valid[[i]] <- is.valid(geoms_v[[i]])
 }
length(which(unlist(geoms_v_valid) == TRUE))
[1] 0 # all complex?

# Okay, make them valid
geoms_v_mV <- list()
 for (m in 1:length(geoms_v)) {
  geoms_v_mV[[m]] <- makeValid(geoms_v[[m]])
  }

# get centroids
geoms_v_centr <- list()
for (r in 1:length(geoms_v_mV)) {
  geoms_v_centr[[r]] <- centroids(geoms_v_mV[[r]])
  }

# see if centroids 'within' poly(s):convex, else concave
geoms_v_within <- list()
_mVfor (s in 1:length(geoms_v_centr)) {
  geoms_v_within[[s]] <- is.related(geoms_v_centr[[s]], geoms_v_mV[[s]], relation = 'within') 
}
> length(which(unlist(geoms_v_within) == TRUE)) #convex
[1] 544
> length(which(unlist(geoms_v_within) == FALSE)) #concave
[1] 176

plot(geoms_v_mV[[25]], col = 'blue') # changes as no set.seed above
plot(geoms_v_centr[[25]], col = 'red', add = TRUE)

This leaves the researcher to ponder which of these 'concave' are desired, that likely depends on the goal of the process that returned these points.

Perhaps 'area' of concave and convex is useful expanse:

expanse(geoms_v_mV[[which(unlist(geoms_v_within) == TRUE)[1]]])
[1] 3
Warning message:
[expanse] unknown CRS. Results can be wrong

With caution...

This only works for convex, not concave polygons though? My current solution is to use `concaveman`, but this will be part of a package and I want to avoid collecting many dependencies... — Zoe, Oct 13 '22 at 13:52

dimfalk · Accepted Answer · 2022-10-14T08:01:13.497

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)

You sorted the points in the correct order though? That's the output I require as I have unordered points or lines. :( If I try this with my unordered points, I get some hexagon-looking shape. — Zoe, Oct 13 '22 at 13:55
Actually, the order shoud be as-is - haven't touched anything in this regard. Have you applied this to `geom` or to `coords`? — dimfalk, Oct 13 '22 at 14:06
I applied it to `geom` because that's the input I have. `coords` is the reqired output. (I wrote it more explicit in the question now, sorry.) — Zoe, Oct 13 '22 at 14:08
@Zoe: Ok, probably I got it now. Please have a look at the revised answer. — dimfalk, Oct 13 '22 at 19:43
It's a bit like magic to me. Do you calculate the angles between the center and the outer points and just order them according to their angle so it's either clock-wise or anti-clockwise depending on the sorting direction (ascending/descending)? Genius :D I can even get this to work without adding dependencies, Thanks! — Zoe, Oct 14 '22 at 05:21
@Zoe: Yeah, that's pretty much exactly what I did - with some filtering for distinct nodes in-between. No magic, just some good old trigonometry.. ;-) I'm not sure how robust this approach is when it comes to more complex geometries, but it seems to work well in this case. You're welcome! — dimfalk, Oct 14 '22 at 08:06

Robert Hijmans · Answer 3 · 2022-10-14T18:08:33.390

Here is a simplified version of falk-env's approach

library(terra)
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))
  

poly_from_seqments <- function(s, crs="") {
    d <- s - colMeans(s)
    angle <- atan2(d[, 1], d[, 2]) 
    s <- s[order(angle), ]  |> unique() 
    vect(as.matrix(s), type="polygons", crs=crs)
}

p <- poly_from_seqments(geom[, c("x", "y")], "+proj=lonlat")

This approach would need some additional work if you have duplicate angles; and it would fail, I think, if the mean location of the coordinates is not inside the polygon.

Unordered lines to closed polygon

3 Answers3