rggplot2voronoispatstatdirichlet

Package for Divide Chain of tesselations in R, spatstat package?


I am trying to create pretty figures of clustered points. Is there a package which will create the divide chain between tessellations of points? Ideally it would be fit for plotting in ggplot.

Here is some example code:

#DivideLineExample
library(spatstat)

W=owin(c(0,1),c(0,1))          # Set up the Window
p<-runifpoint(42, win=W)       # Get random points

ll=cbind(p$x,p$y)              # get lat/long for each point
zclust=kmeans(ll,centers=4)    # Cluster the points spatially into 4 clusters

K<-pp<-D<-list()
plot(W,main="Clustered Points")
for (i in 1:4){                   # this breaks up the points into separate ppp objects for each cluster
  K[[i]]=ll[zclust$cluster==i,]   
  pp[[i]]=as.ppp(K[[i]],W)
  plot(pp[[i]],col=i,add=TRUE,cex=1.5,pch=16)
  D[[i]]=dirichlet(pp[[i]])       # This performs the Dirichlet Tessellation and plots
  plot(D[[i]],col=i,add=TRUE)
}

This outputs as such: http://imgur.com/CCXeOEB

Clusters of points without divisions

What I'm looking for is this: https://i.sstatic.net/m7G5i.jpg

Clusters of points with divide chains

I know an algorithm exists.

Any ideas/alternatives?


Solution

  • I have written a function that I think will do what you want:

    divchain <- function (X) {
        stopifnot(is.ppp(X))
        if(!is.multitype(X)) {
            whinge <- paste(deparse(substitute(X)),
                            "must be a marked pattern with",
                            "factor valued marks.\n")
            stop(whinge)
        }
        X <- unique(X, rule = "deldir", warn = TRUE)
        w <- Window(X)
        require(deldir)
        dd <- deldir(X,z=marks(X),rw=c(w$xrange,w$yrange))
        if (is.null(dd)) 
            return(NULL)
        ddd <- dd$dirsgs
        sss <- dd$summary
        z   <- sss[["z"]]
        rslt <- list()
        nsgs <- nrow(ddd)
        K <- 0
        for (i in 1:nsgs) {
             i1 <- ddd[i,5]
             i2 <- ddd[i,6]
             c1 <- z[i1]
             c2 <- z[i2]
             if(c1 != c2) {
                 K <- K+1
                 rslt[[K]] <- unlist(ddd[i,1:4])
             }
        }
        class(rslt) <- "divchain"
        attr(rslt,"rw") <- dd$rw
        rslt
    }
    

    I have also written a plot method for class "divchain":

    plot.divchain <- function(x,add=FALSE,...){
        if(!add) {
            rw <- attr(x,"rw")
            plot(0,0,type="n",ann=FALSE,axes=FALSE,xlim=rw[1:2],ylim=rw[3:4])
            bty <- list(...)$bty
            box(bty=bty)
        }
        lapply(x,function(u){segments(u[1],u[2],u[3],u[4],...)})
        invisible()
    
    }
    

    E.g.:

    require(spatstat)
    set.seed(42)
    X <- runifpoint(50)
    z <- factor(kmeans(with(X,cbind(x,y)),centers=4)$cluster)
    marks(X) <- z
    dcX <- divchain(X)
    plot(dirichlet(X),border="brown",main="")
    plot(X,chars=20,cols=1:4,add=TRUE)
    plot(dcX,add=TRUE,lwd=3)
    

    Dividing chain

    Let me know whether this is satisfactory. Sorry I can't help you with ggplot stuff; I don't do ggplot.