R/concave.hull.r

In jae0/ecmgis: Package that interacts with stmv.

  concave.hull = function( xy, k=5, ub=NULL, plot=FALSE, random.start=FALSE ) {

    if (k < 3) return (NULL)
    p = as.matrix( xy )
    dimnames(p) =NULL
    p = p[ ! duplicated( p) , ]
    if (is.null(ub)) {
      dx = diff(range(xy[,1]))
      dy = diff(range(xy[,2]))
      dd = max( dx, dy )
      ub = dd / 8
    }

    nsets = nrow (p)
    if ( nsets < 3 ) return ( NULL )
    if ( nsets == 3) return( p )

    # initiate list with the first few points
    convex.hull <- spatstat::convexhull.xy( p )

    p0i = which.min( convex.hull$bdry[[1]]$y )
    if (random.start) {
      p0i = floor(runif(1) * length(convex.hull$bdry[[1]]$y) ) + 1
    }
    p.start = as.matrix( cbind( convex.hull$bdry[[1]]$x[ p0i ], convex.hull$bdry[[1]]$y[ p0i ] ) )
    hull = p0 = p.start
    to.remove = which( p0[,1] == p[,1] & p0[,2] == p[,2] )
    if (length( to.remove) > 0 ) p = p[ - to.remove ,]

    finished = F
    while ( ! finished ) {
      if (nrow(hull)==3) p = rbind( p.start, p)  # return p.start after a few steps
      p.new=NULL
      p.new = nearest.sorted.angles( p, hull, k, eps=pi/10, ub=ub ) # eps determines the acute angle solutions to ignore
      if (is.null(p.new) ) {
        kadd = k
        while (is.null(p.new) ) {
          kadd = kadd + 1
          p.new = nearest.sorted.angles( p, hull, k=kadd, eps=pi/10, ub=ub )
          if ((kadd-k) > 10) stop()
        }
      }
      if (p.new[,1] == p.start[,1] & p.new[,2] == p.start[,2]) finished = T
      p0 = p.new
      hull = rbind( hull, p0 )
      to.remove = which( p0[,1] == p[,1] & p0[,2] == p[,2] )
      if (length( to.remove) > 0 ) p = p[ - to.remove ,]
#      lines(hull[,1], hull[,2])
    }

    sa = signif(splancs::areapl(hull),3)
    densit = signif( nsets/splancs::areapl(hull),3)

    print( paste("Surface area:", sa))
    print( paste("Mean density:", densit))
    print( paste("Neighbourhood radius:", ub))
    print( "Check the plot and modify ub to alter the complexity of the outline." )

    if (plot) {
      plot( xy, pch=20, col="gray",
        sub = paste("Criteria:", k, "neighbours and", ub, "km radius", "\n No. sets =", nsets, "; SA = ", sa, "km^2 \n Density = ", densit, "stations / km^2" ),
        ylab="", xlab="" )
      lines(hull, lwd=2)
    }

    return (hull)
  }