ggplot2+シェープファイル+穴

ggplot2でポリゴンのドテッ腹に風穴はあくか？ http://blogs.yahoo.co.jp/igproj_fusion/16633531.html

↑良い記事です。

ggplot2でポリゴンに穴を開けます。

gridでのgrid.pathで穴あきポリゴンを書けます。 fillルールはR Journal最新号“It's Not What You Draw,It's What You Don't Draw” by Paulで。 http://journal.r-project.org/current.html

ggplot2のgeom_polygonはgrid.polygonで実装されているので、使えません。 代わりにここではgeom_holygonを作って使います（と言ってもgrid.polygonをgrid.pathに変えてrule引数を作っただけです）。

geom_holygonを定義

library(ggplot2)
library(proto)

geom_holygon <- function(mapping = NULL, data = NULL, stat = "identity", position = "identity", 
    rule = "winding", ...) {
    GeomHolygon$new(mapping = mapping, data = data, stat = stat, position = position, 
        rule = rule, ...)
}

GeomHolygon <- proto(ggplot2:::GeomPolygon, {
    objname <- "holygon"

    draw <- function(., data, scales, coordinates, rule, ...) {
        n <- nrow(data)
        if (n == 1) 
            return()

        munched <- coord_munch(coordinates, data, scales)
        munched <- munched[order(munched$group), ]

        first_idx <- !duplicated(munched$group)
        first_rows <- munched[first_idx, ]

        ggname(.$my_name(), gTree(children = gList(pathGrob(munched$x, munched$y, 
            default.units = "native", id = munched$group, rule = rule, gp = gpar(col = first_rows$colour, 
                fill = alpha(first_rows$fill, first_rows$alpha), lwd = first_rows$size * 
                  .pt, lty = first_rows$linetype)))))
    }
})

描画

library(maptools)

## Loading required package: foreign

## Loading required package: sp

## Loading required package: grid

## Loading required package: lattice

## Checking rgeos availability: TRUE

p524072 <- readShapePoly("p524072.shp")
p524072p <- subset(p524072, p524072$HANREI_C == 570400)
p524072p_df <- fortify(p524072p)

## Regions defined for each Polygons

ggplot(p524072p_df, aes(long, lat, group = group)) + geom_holygon(fill = "skyblue") + 
    coord_map()

polygon版

library(maptools)
p524072 <- readShapePoly("p524072.shp")
p524072p <- subset(p524072, p524072$HANREI_C == 570400)
p524072p_df <- fortify(p524072p)

## Regions defined for each Polygons

ggplot(p524072p_df, aes(long, lat, group = group)) + geom_polygon(fill = "skyblue") + 
    coord_map()

holes are emphasized.

library(maptools)
p524072 <- readShapePoly("p524072.shp")
p524072p <- subset(p524072, p524072$HANREI_C == 570400)
p524072p_df <- fortify(p524072p)

## Regions defined for each Polygons

ggplot(p524072p_df, aes(long, lat, group = group, fill = hole)) + geom_polygon(alpha = 0.5) + 
    coord_map()

overlap on google map

library(ggmap)
gm <- get_googlemap(c(mean(range(p524072p_df$long)), mean(range(p524072p_df$lat))), 
    maptype = "hybrid", zoom = 12)
ggmap(gm) + geom_holygon(data = p524072p_df, mapping = aes(long, lat, group = group), 
    fill = "skyblue") + coord_map(xlim = range(p524072p_df$long), ylim = range(p524072p_df$lat))