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.
