自分のための「Rで地理空間3D樹状図」

このコードは、シェープファイルを読み込んで、
平面上の距離から階層的クラスタリングを行い、
さらに、その結果を地理空間上に投影するというもの。

このソースコードは、私の「思い付き」に無理を言って後輩を付き合わせ、
私は、ソースコードをチェックして書き直しただけ。
同志社大学大学院文化情報学研究科の谷岡さん、長々とありがとう。

本人に許可をもらったので、このままUP。
もちろん、私も谷岡さんも一切の責任を取りません。
距離法のところは注意。Ward Method のときは、ユークリッド「平方」距離。

上記注意事項の意味が解らない人は、使用を避けた方が良いかもしれない。
他にも、細々とした所に問題があるのは理解している。
比較的、可読性は高いと思うので、見れば解るだろう。

なお、二人で相談した結果、論文クオリティに足りないということで、
このソースコードは「使い道」を知っている人に使って頂きたい。
使う時には、このブログを参考文献に...。

ソースコードを自分の用途に書き換えれば、多少は使えるかもしれない。

とりあえず、今回のサンプルデータは以下の場所に。
http://docs.google.com/open?id=0B9OtIs5BjRVhdWJmWU43V0Z3eVE

実行コードは以下の通り。

---

# Read shape file.
hoge <- readShapePoints("Krugans.shp")

# Number cluster. k <- 3
clus.method <- "average"

# Build 3D dendrogram. dendro.3d(hoge, k, clus.method)

ソースコードは以下の通り。

---

dendro.3d <- function(shp, k, c.method){
    # Read packages for this function.
    require(rgl)
    require(maptools)

    # Get coordinates of points data.
    coords <- shp@coords
    oldnames <- rownames(coords)
    rownames(coords) <- 1:nrow(coords)
    
    # Add nodes corresponding clusters.
    rowcount <- nrow(coords)

    # Conduct cluster analysis
    coords.dist <- dist(coords)
    coords.clus <- hclust(coords.dist, method=c.method)

    # Get a list of name of merge.
    coords.merge <- coords.clus$merge
    coords.merge2 <- coords.clus$merge
    coords.height <- coords.clus$height

    # Add a new column to store Z values.
    coords <- cbind(coords, rep(0, nrow(coords)))
    colnames(coords) <- c("x","y","z")

    # Change merged nodes name
    for(i in 1:nrow(coords.merge)){
            if(as.numeric(coords.merge[i,1]) > 0){
                rownum <- as.numeric(coords.merge[i,1])
                node_org <- coords.merge[rownum,1]
                node_dst <- coords.merge[rownum,2]
                node_new <- paste(node_org, node_dst, sep="")
                coords.merge[i,1] <- node_new
            }
            if(as.numeric(coords.merge[i,2] > 0)){
                rownum <- as.numeric(coords.merge[i,2])
                print(rownum)
                print(i)
                node_org <- coords.merge[rownum,1]
                node_dst <- coords.merge[rownum,2]
                node_new <- paste(node_org, node_dst, sep="")
                coords.merge[i,2] <- node_new
            }
    }

    # Combine a list of merged nodes and a list of height.
    coords.merge_and_height <- cbind(coords.merge, coords.height)

    # Rename rownames.
    rownames(coords) <- -1*as.numeric(rownames(coords))
    
    for(i in 1:(rowcount-1)){
       # print(i)
        # Get z value from height of cluster.
        z <- coords.merge_and_height[i,3]

        # Create a new node name.
        node_org <- coords.merge_and_height[i, 1]
        node_dst <- coords.merge_and_height[i, 2]
        newname <- paste(node_org, node_dst, sep="")

        # Add new column with new node name.
        coords <- rbind(coords, as.numeric(c(0,0,z)))
        rownames(coords)[i+rowcount] <- newname

        x_org <- as.numeric(coords[rownames(coords)==node_org,1])
        x_dst <- as.numeric(coords[rownames(coords)==node_dst,1])
        x_new <- (x_org + x_dst)/2

        y_org <- as.numeric(coords[rownames(coords)==node_org,2])
        y_dst <- as.numeric(coords[rownames(coords)==node_dst,2])
        y_new <- (y_org + y_dst)/2

        coords[i+rowcount, 1] <- x_new
        coords[i+rowcount, 2] <- y_new
        coords[i+rowcount, 3] <- as.numeric(z)
    }
    # Define colours for nodes.
    nodes.col <- as.numeric(cutree(coords.clus,k))+1
    nodes.marged.col <- c(nodes.col, rep(1,length(coords.height)))
    
    # Plot nodes.
    plot3d(coords,col=nodes.marged.col)
    
    # Plot branches.
    for(i in 1:(rowcount-1)){
        hoge1 <- coords[(i+nrow(shp)),]
        segments3d(rbind(coords[rownames(coords)==coords.merge[i,1],],c(coords[rownames(coords)==coords.merge[i,1],1:2],hoge1[3])),col="orange")
        segments3d(rbind(coords[rownames(coords)==coords.merge[i,2],],c(coords[rownames(coords)==coords.merge[i,2],1:2],hoge1[3])),col="orange")
        segments3d(rbind(c(coords[rownames(coords)==coords.merge[i,1],1:2],hoge1[3]),c(coords[rownames(coords)==coords.merge[i,2],1:2],hoge1[3])),col="orange")
    }
}