[R-sig-Geo] alpha hull (ahull) to polygon shapefile

[Andrew Bevan]

Dear Kevin

As far as I remember, the ahull class object stores its information in arcs that are not a straight conversion to standard GIS-style lines or polygons. I played with the following function some time ago and it is very rough (apologies in advance), but it works for my purposes. Its starting point was a script for conversion to lines only by Dylan Beaudette that the latter may be sufficient for you (http://casoilresource.lawr.ucdavis.edu/drupal/node/919). In contrast, the version below models the alphahull arcs with a user-defined line segment increment and tries to borrow a polygon-building method from the ASDAR book to create a SpatialPolygonsDataFrame object. Both the function and a toy example below.

Andy

########################

ah2sp <- function(x, increment=360, rnd=10, proj4string=CRS(as.character(NA))){
  require(alphahull)
  require(maptools)
  if (class(x) != "ahull"){
    stop("x needs to be an ahull class object")
  }
  # Extract the edges from the ahull object as a dataframe
  xdf <- as.data.frame(x$arcs)
  # Remove all cases where the coordinates are all the same      
  xdf <- subset(xdf,xdf$r > 0)
  res <- NULL
  if (nrow(xdf) > 0){
    # Convert each arc to a line segment
    linesj <- list()
    prevx<-NULL
    prevy<-NULL
    j<-1
    for(i in 1:nrow(xdf)){
      rowi <- xdf[i,]
      v <- c(rowi$v.x, rowi$v.y)
      theta <- rowi$theta
      r <- rowi$r
      cc <- c(rowi$c1, rowi$c2)
      # Arcs need to be redefined as strings of points. Work out the number of points to allocate in this arc segment.
      ipoints <- 2 + round(increment * (rowi$theta / 2),0)
      # Calculate coordinates from arc() description for ipoints along the arc.
      angles <- anglesArc(v, theta)
      seqang <- seq(angles[1], angles[2], length = ipoints)
      x <- round(cc[1] + r * cos(seqang),rnd)
      y <- round(cc[2] + r * sin(seqang),rnd)
      # Check for line segments that should be joined up and combine their coordinates
      if (is.null(prevx)){
        prevx<-x
        prevy<-y
      } else if (x[1] == round(prevx[length(prevx)],rnd) && y[1] == round(prevy[length(prevy)],rnd)){
          if (i == nrow(xdf)){
          #We have got to the end of the dataset
            prevx<-append(prevx,x[2:ipoints])
            prevy<-append(prevy,y[2:ipoints])
            prevx[length(prevx)]<-prevx[1]
            prevy[length(prevy)]<-prevy[1]
            coordsj<-cbind(prevx,prevy)
            colnames(coordsj)<-NULL
            # Build as Line and then Lines class
            linej <- Line(coordsj)
            linesj[[j]] <- Lines(linej, ID = as.character(j))
          } else {
            prevx<-append(prevx,x[2:ipoints])
            prevy<-append(prevy,y[2:ipoints])
          }
        } else {
      # We have got to the end of a set of lines, and there are several such sets, so convert the whole of this one to a line segment and reset.
          prevx[length(prevx)]<-prevx[1]
          prevy[length(prevy)]<-prevy[1]
          coordsj<-cbind(prevx,prevy)
          colnames(coordsj)<-NULL
      # Build as Line and then Lines class
          linej <- Line(coordsj)
          linesj[[j]] <- Lines(linej, ID = as.character(j))
          j<-j+1
          prevx<-NULL
          prevy<-NULL
        }
    }
    # Promote to SpatialLines
    lspl <- SpatialLines(linesj)
    # Convert lines to polygons
    # Pull out Lines slot and check which lines have start and end points that are the same
    lns <- slot(lspl, "lines")
    polys <- sapply(lns, function(x) { 
      crds <- slot(slot(x, "Lines")[[1]], "coords")
      identical(crds[1, ], crds[nrow(crds), ])
    }) 
    # Select those that do and convert to SpatialPolygons
    polyssl <- lspl[polys]
    list_of_Lines <- slot(polyssl, "lines")
    sppolys <- SpatialPolygons(list(Polygons(lapply(list_of_Lines, function(x) { Polygon(slot(slot(x, "Lines")[[1]], "coords")) }), ID = "1")), proj4string=proj4string)
    # Create a set of ids in a dataframe, then promote to SpatialPolygonsDataFrame
    hid <- sapply(slot(sppolys, "polygons"), function(x) slot(x, "ID"))
    areas <- sapply(slot(sppolys, "polygons"), function(x) slot(x, "area"))
    df <- data.frame(hid,areas)
    names(df) <- c("HID","Area")
    rownames(df) <- df$HID
    res <- SpatialPolygonsDataFrame(sppolys, data=df)
    res <- res[which(res at data$Area > 0),]
  }  
  return(res)
}


########################

#So for example:
library(alphahull)
library(maptools)
n <- 300
theta <- runif(n,0,2*pi)
r <- sqrt(runif(n,0.25^2,0.5^2))
x <- cbind(0.5+r*cos(theta),0.5+r*sin(theta))
a <- ahull(x, alpha = 0.1)
plot(a)

#Convert
b <- ah2sp(a)
plot(b, col="lightgrey", pbg="white")

#Write to shapefile
writeSpatialShape(b, "b")



On 4 Mar 2012, at 18:52, Kevin Goulding wrote:

> I would like to convert a point shapefile into a polygon shapefile
> containing an alpha hull of the points.  I've imported the point shapefile
> and calculated the desired alpha hull:
> 
> #######################################
> require(spdep)
> require(alphahull)
> 
> ccPoints=readShapeSpatial('myshapefile')
> a=5000
> 
> x.coords <- coordinates(ccPoints)
> x.ah <- ahull(x.coords[,1], x.coords[,2], alpha=a)
> 
> # plot the concave alpha hull:
> # plot(x.ah,col='purple',add=T,asp=1,cex=0, pch=4)
> # looks good...
> 
> # now, would like to export alpha hull as a shapefile with a single polygon
> # something like this:
> # write.shp(x.ah, file = 'ahull-shapefile.shp')
> #######################################
> 
> How do I output a polygon shapefile containing the alpha hull?
> 
> Many thanks, Kevin
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>