2

How can I convert an open linestring into a closed polygon, but when the points are not in order?

Here is the easy case, when the points are in order (from an sf vignette):

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
mat <- matrix(c(0,0,
                0,1,
                1,1,
                0.5,0.5
),,
2,
byrow=TRUE)
ls = st_linestring(mat)
plot(ls) # an open curve

st_polygonize(ls) # fails
#> GEOMETRYCOLLECTION EMPTY

It seems that st_polygonize lacks the ability to automatically close the polygon.

In this simple case, we can easily take the first point and add it to the end of the matrix, then use st_polygonize to convert to a polygon

mat <- rbind(mat, mat[1,])

ls = st_linestring(mat)
plot(ls) # an open curve

pol <- st_polygonize(ls) # works
plot(pol)

But what about a more complex situation where the points are not ordered, and we cannot easily find the points at each side of the opening? Here's an example similar to my problem:

eg <-  structure(c(77, 1967, 656.689484908856, 1417.41686984243, 2022.29766310616, 
                   2115.48512309123, 255.539657579976, 712.790838796637, 1551.34291678482, 
                   256.301574078794, 2254.86461653393, 219.398084314749, 308.336271812538, 
                   210.641580119875, 124.755913242167, 65.3263108593827, 2244, 23, 
                   2087.68279235023, 1616.80814193697, 2099.5290296304, 1544.41763616998, 
                   1966.27425582228, 1409.51461937173, 1626.3609089262, 27, 2256, 
                   1166.68169373894, 28, 2016, 1316.00439547048, 1856.79371972606, 
                   159.310515091144, 40.4168698424283, 1937.29766310616, 1482.51487690877, 
                   475.460342420024, 106.209161203363, 93.3429167848175, 2004.69842592121, 
                   938.864616533932, 1589.39808431475, 1891.66372818746, 2062.64158011987, 
                   1456.75591324217, 1339.32631085938, 1093.34356771967, 857.481418293397, 
                   1495.31720764977, 127.808141936968, 554.529029630397, 90, 370.274255822284, 
                   38, 133, 1214.67021042315, 978.302571822642, 20, 1122.11042488164, 
                   1970.91320285213), .Dim = c(30L, 2L), .Dimnames = list(NULL, 
                                                                          c("X", "Y")))
plot(eg, asp = 1)

We can clearly see this is a circular shape with a flared opening at the top. I want to turn it into a closed polygon by drawing a line across the opening, then use st_sample to sample some points from that line so I have a set of points that describes a closed shape.

If we look at the coords for this example, we see that the first and last rows do not correspond to the points at the opening

eg[1,] # way off 
#>        X        Y 
#>   77.000 1316.004
eg[nrow(eg),] # close, but can we be sure?
#>        X        Y 
#> 2016.000 1970.913

We can see more clearly that these points are not in an order that reflects the shape by plotting:

plot(eg, type = 'l')

I would like to know how can we put the points in order that reflects the shape, so we can get the first point of the line (i.e. upper left of the opening) and copy it to the end of the matrix, so that st_polygonize can close the polygon?

devtools::session_info()
#> Session info -------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.5.0 (2018-04-23)
#>  system   x86_64, darwin15.6.0        
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  tz       America/Los_Angeles         
#>  date     2018-07-20
#> Packages -----------------------------------------------------------------
#>  package   * version date       source        
#>  backports   1.1.2   2017-12-13 CRAN (R 3.5.0)
#>  base      * 3.5.0   2018-04-24 local         
#>  class       7.3-14  2015-08-30 CRAN (R 3.5.0)
#>  classInt    0.2-3   2018-04-16 CRAN (R 3.5.0)
#>  compiler    3.5.0   2018-04-24 local         
#>  curl        3.2     2018-03-28 CRAN (R 3.5.0)
#>  datasets  * 3.5.0   2018-04-24 local         
#>  DBI         1.0.0   2018-05-02 CRAN (R 3.5.0)
#>  devtools    1.13.6  2018-06-27 cran (@1.13.6)
#>  digest      0.6.15  2018-01-28 CRAN (R 3.5.0)
#>  e1071       1.6-8   2017-02-02 CRAN (R 3.5.0)
#>  evaluate    0.10.1  2017-06-24 CRAN (R 3.5.0)
#>  graphics  * 3.5.0   2018-04-24 local         
#>  grDevices * 3.5.0   2018-04-24 local         
#>  grid        3.5.0   2018-04-24 local         
#>  htmltools   0.3.6   2017-04-28 CRAN (R 3.5.0)
#>  httr        1.3.1   2017-08-20 CRAN (R 3.5.0)
#>  knitr       1.20    2018-02-20 CRAN (R 3.5.0)
#>  magrittr    1.5     2014-11-22 CRAN (R 3.5.0)
#>  memoise     1.1.0   2017-04-21 CRAN (R 3.5.0)
#>  methods   * 3.5.0   2018-04-24 local         
#>  mime        0.5     2016-07-07 CRAN (R 3.5.0)
#>  R6          2.2.2   2017-06-17 CRAN (R 3.5.0)
#>  Rcpp        0.12.17 2018-05-18 CRAN (R 3.5.0)
#>  rmarkdown   1.10    2018-06-11 CRAN (R 3.5.0)
#>  rprojroot   1.3-2   2018-01-03 CRAN (R 3.5.0)
#>  sf        * 0.6-3   2018-05-17 CRAN (R 3.5.0)
#>  spData      0.2.9.0 2018-06-17 CRAN (R 3.5.0)
#>  stats     * 3.5.0   2018-04-24 local         
#>  stringi     1.2.2   2018-05-02 CRAN (R 3.5.0)
#>  stringr     1.3.1   2018-05-10 CRAN (R 3.5.0)
#>  tools       3.5.0   2018-04-24 local         
#>  units       0.6-0   2018-06-09 CRAN (R 3.5.0)
#>  utils     * 3.5.0   2018-04-24 local         
#>  withr       2.1.2   2018-03-15 CRAN (R 3.5.0)
#>  xml2        1.2.0   2018-01-24 CRAN (R 3.5.0)
#>  yaml        2.1.19  2018-05-01 CRAN (R 3.5.0)
  • By definition, a polygon is bounded by an ordered set of vertices. Without order, you've got an NP-Hard problem on your hands, possibly NP-Complete. If you can make some assumptions, you might be able to use a convex or concave hull on the point set, but it might not be accurate. – Vince Jul 21 '18 at 1:18
  • Thanks, that suggests I'm not missing a simple trick here... I was thinking something like this might work gis.stackexchange.com/q/138635/9803 – Ben Jul 21 '18 at 3:21
  • What might help is articulating the algorithm you would use to join the dots yourself - imagine you’re using pen and paper. How to choose a starting point? How to then choose the next point, and the next? The latter could simply be "choose the closest point and connect to that" for your sample data. – Simbamangu Jul 22 '18 at 4:29
  • 1
    How about alpha hull/alpha shapes/concave hull gis.stackexchange.com/questions/1200/…? – user30184 Aug 7 '18 at 8:19
3

Following off the concave hull comment, here is a solution using the concaveman package:

library(sf)
library(concaveman)
library(ggplot2)

eg <- st_as_sf(data.frame(eg), coords = c("X", "Y"))
poly <- concaveman(eg)

ggplot() + 
  geom_sf(data = poly) +
  geom_sf(data = eg)

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.