Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Row ordering for polygons

Tags:

r

My question is simple. Is there an automatic way to order you data so that it makes "clean" polygons? I have functions that are generating rings (specifically the ahull function), and I would like a way to cleanly produce polygons using such functions. Here is an example.

x <- c(1:3, 3:1, 1)
y <- c(1,1,1,3,3,2, 1)
xy <- cbind(x,y)

Sr1 <- Polygon(xy)
Srs1 = Polygons(list(Sr1), "s1")
SpP = SpatialPolygons(list(Srs1)) 
plot(SpP) 

z <- runif(7)
xyz <- cbind(x,y,z)
xyz <- xyz[order(z),]

xy <- xyz[,-3]
xy <- rbind(xy, xy[1,])

Sr1 <- Polygon(xy)
Srs1 = Polygons(list(Sr1), "s1")
SpP = SpatialPolygons(list(Srs1)) 
SpP = SpatialPolygons(list(Srs1)) 
plot(SpP)

Here is my real data: https://drive.google.com/file/d/0B8QG4cbDqH0UOUlobnlWaDgwOWs/edit?usp=sharing

like image 583
mgslocum Avatar asked Jun 10 '14 14:06

mgslocum


1 Answers

In a sense, you have answered your own question.

Assuming you have a set of points, and you use ahull(...) in the alphahull package to generate the convex hull, you can extract the points on the boundary, in the correct order, directly from the ahull object. Here is an example:

library(sp)
library(alphahull)
set.seed(1)            # for reproducible example
X <- rnorm(100)
Y <- rnorm(100)
plot(X,Y)
XY <- cbind(X,Y)
hull <- ahull(XY,alpha=1)
plot(hull)

# extract the row numbers of the boundary points, in convex order.
indx=hull$arcs[,"end1"]  
points <- XY[indx,]                  # extract the boundary points from XY
points <- rbind(points,points[1,])   # add the closing point
# create the SpatialPolygonsDataFrame
SpP = SpatialPolygons(list(Polygons(list(Polygon(points)),ID="s1"))) 
plot(SpP)
points(XY)

EDIT Response to OP's providing their dataset.

ahull(...) seems to fail, without warning, with your dataset - it does not produce any convex hulls. After a bit if experimentation, it looks like the problem has to do with the magnitude of the x,y values. If I divide everything by 1000, it works. No idea what's going one with that (perhaps someone else will provide an insight??). Anyway, here's the code and the result:

library(sp)
library(alphahull)
df <- read.csv("ahull problem.csv")
hull <- ahull(df[2:3]/1000,alpha=2)
plot(hull)
# extract the row numbers of the boundary points, in convex order.
indx=hull$arcs[,"end1"]  
points <- df[indx,2:3]               # extract the boundary points from df
points <- rbind(points,points[1,])   # add the closing point
# create the SpatialPolygonsDataFrame
SpP = SpatialPolygons(list(Polygons(list(Polygon(points)),ID="s1"))) 
plot(SpP)
points(df[2:3])

Note also that alpha=2. Setting alpha=1 with this dataset actually generates 2 hulls, one with 1 point and one with all the other points. Setting alpha=2 creates 1 hull.

like image 111
jlhoward Avatar answered Oct 20 '22 03:10

jlhoward