2012/09/23

自分のための「RでShape file の表示」

Rでは、plot()関数やspplot()関数など、空間データの表示方法がいくつかある。
ただ、中々、痒い所に手が届かず、色々と不満。
そんな訳で、自分用の描画デバイスを作ってみようかと...。

今回のサンプルデータは以下の場所に。
https://docs.google.com/folder/d/0B9OtIs5BjRVhWm8zV1hCUFh6YnM/edit

とりあえず、表示は上手く行く。あとは、凝った表現。
どうやら、図の重ね合わせ等もできそう。
ためしに、どのようなものか挑戦してみた。やはり、若干面倒か?

以下が実行コード。

# Load a required library.
library(maptools)

# Read a shape file.
shp <- readShapePoly('/home/yufujimoto/Desktop/japan.shp')

# Draw polygons.
drawShape.poly(shp)

そして、以下がソースコード

drawShape.poly <- function(shp){
    require(grid)
 
    # Get range from bounding box.
    range.x <- c(shp@bbox[1,1],shp@bbox[1,2])
    range.y <- c(shp@bbox[2,1],shp@bbox[2,2])

    # Get widths and heights from data.
    width <- diff(range.x)
    height <- diff(range.y)

    pushViewport(
        viewport(name="layer",
                layout=grid.layout(1,1,
                       widths=width,
                       heights=height,
                       respect=TRUE)
                 )
    )
    pushViewport(
        viewport(name="viewport",
                layout.pos.row=1,
                layout.pos.col=1,
                xscale=range.x,
                yscale=range.y,
                clip=TRUE)
    )
    index <- 1

    gp <- gpar(col <- 3)

    # Draw lines for shape files.
    for(i in shp@polygons) {
        for(j in i@Polygons){
            name <- i@ID
            x <- j@coords[,1]
            y <- j@coords[,2]
            grid.polygon(x, y, default.unit="native",name=name, gp=gpar(fill="red"))
            index <- index + 1
        }
    }

    # Draw bounding box.
    name <- "bbox"
    x <- c(shp@bbox[1,1], shp@bbox[1,1], shp@bbox[1,2], shp@bbox[1,2], shp@bbox[1,1])
    y <- c(shp@bbox[2,1], shp@bbox[2,2], shp@bbox[2,2], shp@bbox[2,1], shp@bbox[2,1])
    grid.polyline(x, y, default.unit="native",name=name)
    upViewport(2)
}

0 件のコメント:

コメントを投稿