1

Given the 2D coordinates of a number of vertices, I would like to calculate the internal angles of the according polygon.

Here's a dataframe with a few polygons and the coordinates of their vertices. The problem is that the vertices aren't ordered.

structure(list(x = c(173L, 173L, 185L, 190L, 231L, 267L, 185L, 
190L, 233L, 260L, 190L, 231L, 260L, 230L, 230L, 172L, 233L, 230L, 
231L, 267L, 185L, 172L, 233L, 231L, 231L), y = c(299L, 299L, 
321L, 360L, 361L, 377L, 321L, 360L, 363L, 309L, 360L, 361L, 309L, 
322L, 322L, 378L, 363L, 322L, 391L, 377L, 321L, 378L, 363L, 391L, 
361L), polygonID = c(2L, 4L, 5L, 3L, 6L, 7L, 2L, 5L, 7L, 6L, 2L, 
3L, 4L, 5L, 6L, 2L, 3L, 4L, 3L, 6L, 4L, 3L, 6L, 7L, 5L)), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 9L, 10L, 12L, 13L, 14L, 15L, 16L, 17L, 
19L, 20L, 21L, 24L, 25L, 27L, 29L, 30L, 31L, 34L), class = "data.frame")

Essentially, I would go from vertex to vertex to calculate the edges/vectors and use those to calculate the internal angles. But it's not clear to me how I would do this in an algorithmic fashion. Any hints much appreciated.

  • Well, if you look at three points at a time you have this: https://stackoverflow.com/questions/1211212/how-to-calculate-an-angle-from-three-points – MrFlick Feb 13 '19 at 16:36
  • It looks like you may have removed the "row.names" component from the `dput` output, that's not optional with `class="data.frame"`. Try pasting this `structure` to your console and see what happens. – r2evans Feb 13 '19 at 16:45
  • 1
    @r2evans - sorry, corrected –  Feb 13 '19 at 16:47
  • @MrFlick So I could calculate the centroid and use that as the third point? Still there's the problem that I need to pick adjacent vertices.. –  Feb 13 '19 at 16:48
  • If you want the internal angles, then you just basically just need to do a rolling window to calculate the angle between consecutive points. I don't see a need to calculate a centroid unless i've misunderstood exactly what you are after. Providing the expected output for a simple reproducible example is helpful so we can test and verify possible answers. – MrFlick Feb 13 '19 at 16:50
  • 2
    It might help to use `chull(x,y)` to get a clockwise ordering of points within a particular `polygonID`. From that and @MrFlick's link, you should be able to iterate over each polygon. – r2evans Feb 13 '19 at 16:59
  • @r2evans - you mean sth like that? `ID = seq(min(DF$polygonID),max(DF$polygonID))` `DForder <- c()` `for (i in ID) {DForder[[i]] <- chull(xy.coords(DF[DF$polygonID==i,]))}` `DF <- DF[order(DF$polygonID),]` `DF$order = unlist(DForder)` `DF <- DF[order(DF$polygonID,DF$order),]` –  Feb 13 '19 at 17:24
  • 1
    I would actually use something like `by(DF, DF$polygonID, function(a) { o <- chull(a[,1:2]); atan2(...)-atan2(...); })` – r2evans Feb 13 '19 at 17:29
  • @r2evans - you mean like this? `by(DF, DF$polygonID, function(a) { o <- chull(a[,1:2]); atan2(a[o+2,2]-a[o,2],a[o+2,1]-a[o,1])-atan2(a[o+1,2]-a[o,2],a[o+1,1]-a[o,1]); })` . I think I do not understand yet how to implement a rolling window. I do not get angles at all vertices.. –  Feb 13 '19 at 19:20

1 Answers1

1

Here's a shot at one:

shp <- structure(list(x = c(173L, 173L, 185L, 190L, 231L, 267L, 185L, 
190L, 233L, 260L, 190L, 231L, 260L, 230L, 230L, 172L, 233L, 230L, 
231L, 267L, 185L, 172L, 233L, 231L, 231L), y = c(299L, 299L, 
321L, 360L, 361L, 377L, 321L, 360L, 363L, 309L, 360L, 361L, 309L, 
322L, 322L, 378L, 363L, 322L, 391L, 377L, 321L, 378L, 363L, 391L, 
361L), polygonID = c(2L, 4L, 5L, 3L, 6L, 7L, 2L, 5L, 7L, 6L, 
                     2L, 3L, 4L, 5L, 6L, 2L, 3L, 4L, 3L, 6L, 4L, 3L, 6L, 7L, 5L)), class = "data.frame",
row.names = c(NA, -25L))

aa <- shp[ shp$polygonID == 2, ]
aa <- aa[ chull(aa[,1:2]), ]
aa
#      x   y polygonID
# 7  185 321         2
# 1  173 299         2
# 16 172 378         2
# 11 190 360         2

Now aa is a 4-point polygon, ordered clockwise. Let's hard-code some calculations for now:

indm1 <- c(ind[-1], ind[1])
indp1 <- c(ind[length(ind)], ind[-length(ind)])
angles <- ((atan2(aa$y[indm1] - aa$y[ind], aa$x[indm1] - aa$x[ind]) -
              atan2(aa$y[indp1] - aa$y[ind], aa$x[indp1] - aa$x[ind])) * 180 / pi) %% 360
cbind(indm1,ind,indp1)
#      indm1 ind indp1
# [1,]     2   1     4
# [2,]     3   2     1
# [3,]     4   3     2
# [4,]     1   4     3
angles
# [1] 158.69530  29.33568  44.27478 127.69424

And let's see (I was initially perplexed that it was not visually-correlating until I realized that the aspect-ratio was off, ergo asp=1).

plot(y~x, data=aa, type='l', asp=1)
with(aa, text(x-5, y, seq_len(nrow(aa)), col="red"))
with(aa, text(x+5, y, round(angles, 0)))

first polygon

Okay, let's try to formalize this a little:

getangles <- function(aa) {
  aa <- aa[chull(aa[,1:2]),]
  ind <- seq_len(nrow(aa))
  indm1 <- c(ind[-1], ind[1])
  indp1 <- c(ind[length(ind)], ind[-length(ind)])
  ((atan2(aa$y[indm1] - aa$y[ind], aa$x[indm1] - aa$x[ind]) -
      atan2(aa$y[indp1] - aa$y[ind], aa$x[indp1] - aa$x[ind])) * 180 / pi) %% 360
}
by(shp, shp$polygonID, getangles)
# shp$polygonID: 2
# [1] 158.69530  29.33568  44.27478 127.69424
# ------------------------------------------------------------ 
# shp$polygonID: 3
# [1] 130.91438 136.39718 133.60282  57.42594  81.65967
# ------------------------------------------------------------ 
# shp$polygonID: 4
# [1]  29.98564  54.83259 119.88349 155.29828
# ------------------------------------------------------------ 
# shp$polygonID: 5
# [1] 92.74183 81.42121 98.70294 87.13402
# ------------------------------------------------------------ 
# shp$polygonID: 6
# [1]  72.44870 111.95989 136.46880 157.38014  61.74247
# ------------------------------------------------------------ 
# shp$polygonID: 7
# [1] 71.70548 64.66388 43.63064

(There might be some issues with rounding/modulus, I'll leave it to you to beautify and verify the others.)

r2evans
  • 141,215
  • 6
  • 77
  • 149
  • 1
    Much obliged. Disappointing to not find a solution by myself this time, but well compensated by being taught about a few useful things! –  Feb 13 '19 at 20:39
  • 1
    This is "one" solution, I suspect there are alternatives. I'm confident that `getangles` could be cleaned up a bit, though I don't know if it can be much faster. I don't use `chull` often but found that its "guarantees" with convex polygons and its performance make it appealing when it works. You may note that it specifically drops points that are inside the convex hull (e.g., `chull(c(4,4,2,1,1), c(4,1,2,1,4))` only uses four points), so use it with caution if you have more-varied data. – r2evans Feb 13 '19 at 21:23
  • you mean it doesn't work for concave polygons, ergo the name 'convex hull', right? :) I will go through it with a fresh brain but I tried it already on more diverse data and it works pretty nicely - population average follows expected trend. Will add a solution later if of additional value here. –  Feb 13 '19 at 21:35
  • One can only know it's for convex hulls if you read the docs (or heed all notes/warnings/cautions), which appears to be an exception here on SO :-) – r2evans Feb 13 '19 at 21:57