【R】いろいろな地図をプロット・作成する方法まとめ【地理空間データ分析】

この記事ではRで様々な地図をプロットする方法を紹介します。Rではドットプロットやチャートマップ・コロプレス図といったものをプロットすることができます。

今回は「maptools」というパッケージを使用して地図をプロットしていきます。

maptoolsをインストールしていない場合は下記のコマンドでインストールを済ませておいてください。

install.packages("maptools")

サンプルデータの読み込み

今回はmaptoolsのパッケージに含まれる米国ノースカロライナ州の乳幼児突然死症候群のデータを使用します。

library(maptools)
pj <- CRS("+proj=longlat +ellps=clrk66")
f <- system.file("shapes/sids.shp", package = "maptools")[1]
nc <- readShapePoly(f, IDvar = "FIPSNO", proj4string = pj)

CRS関数は、オープンソースの座標操作系ライブラリPROJ.4で定められている書式で、CRSクラスを作成するコマンドになります。詳しくはこちらを参照してください。

readShapePolyは、side.shpというESRIシェープファイル形式のファイルを読み込んで、ncというSpetialPolygonsDataFrameオブジェクトを作成しています。(こちらの詳細については別記事で紹介したいと思います。)

 

これでデータの読み込みが完了したので、ここからは地図のプロットをしてみましょう。

 

ドットマップのプロット

先ほど読み込んだデータ、ノースカロライナ州の幾何中心点をドットマップとしてプロットする。

plot(coordinates(nc), asp = 1, xlab = "Longitude", ylab = "Latitude")

次に、地点データがある属性を持つことを考えてみます。1974年から1978年の間に生まれた新生児数を上記の幾何中心に付与してプロットしてみます。

# spplotを使うため、spクラスに変換する
nc.b74 <- SpatialPointsDataFrame(SpatialPoints(coordinates(nc)), data = data.frame(nc@data$BIR74))
spplot(nc.b74)

ドットマップをプロットすることはできましたが、このドットマップの問題点として、点の重なりをうまく表現できない問題があります。例えば、多くの点が同じ地点に重なり合う場合、その重なりが2つ分なのかそれとも10個分なのかはドットマップからは判断ができないです。

 

この解決方法としてサーフェイス化という物があり、こちらは別の記事で紹介したいと思います。

 

ドットマップのプロット

ドット密度マップは領域内で数えられた数値を表現するとき、領域内にランダムに点をうつことで、その地理的分布を表現することができます。

# この関数に与える変数は数え上げた数値を前提としているので整数型に変換
d <- dotsInPolys(nc, as.integer(nc@data$SID74))
plot(nc, axes = TRUE)
title(main = "SIDS 1974-78")
plot(d, pch = 19, cex = 0.5, col = "red", add = TRUE)

チャートマップのプロット

チャートマップは、統計グラフを地図上に配置することで、グラフのパターンと地理的分布パターンを同時に把握することができます

統計グラフに適した種類のデータであれば、どのようなデータでもこの方法を使うことができるようです。

 

今回は、ncオブジェクトにあるdataスロットに、ノースカロライナ州における1974-1978年の群別出生数(BIR74)とSIDS死亡数(SID74)、1979-1984年の群別出生数(BIR79)とSIDS死亡数(SID79)が含まれているので、こちらのデータを使って死亡率を計算します。

d74 <- nc@data$SID74 / nc@data$BIR74 * 1e+05
summary(d74)
 Min. 1st Qu. Median Mean 
0.0 108.4 185.5 204.6 
3rd Qu. Max. 
260.4 955.4 
d79 <- nc@data$SID79 / nc@data$BIR79 * 1e+05
summary(d79)
 Min. 1st Qu. Median Mean 
0.0 124.9 207.5 203.9 
3rd Qu. Max. 
253.9 611.4 

この2つの期間の死亡率を棒グラフとして地図上にプロットしてみます。

library(TeachingDemos)
d <- cbind(d74, d79)
d.m <- max(d)
cols <- c(gray(0.3), gray(0.8))
plot(nc)
nc.c <- coordinates(nc)
for (i in 1:nrow(nc.c)) {
  subplot(
    barplot(
    d[i, ], 
    ylim = c(0, d.m), 
    col = cols,
    names = c("", ""), 
    yaxt = "n"
  ), 
  x = nc.c[i, 1],
  y = nc.c[i, 2], 
  size = c(0.1, 0.5), 
  vadj = 0
  )
}
legend(-83.5, 35, c("1974-1978", "1979-1984"), fill = cols, bty = "n", cex = 0.8)

コロプレス図のプロット

コロプレス図は、段彩図、階級区分地図、塗り分け地図とも呼ばれる基本的な地図表現です。領域データかつ集計データを扱う場合に使用します。

spplot(nc, zcol = c("BIR74", "NWBIR74", "BIR79", "NWBIR79"),
names.attr = c("Births, 1974-78", "Non-white births, 1974-78",
"Births, 1979-84", "Non-white births, 1979-84"))

spplotを使うことでSpatialPolygonsDataFrameオブジェクトから簡単にコロプレス図を作成することができます。

自分で作成した変数に対してコロプレス図を作成したい場合は下記のようにすることで実現できます。(死亡率のコロプレス図)

# デフォルトの色がダサいのでカラーパレットを使う
library(RColorBrewer)
d74 <- nc@data$SID74 / nc@data$BIR74 * 1e+05
nc@data <- cbind(nc@data, d74)
palette <- brewer.pal(n = 7, name = "OrRd")
spplot(nc, zcol = "d74", col.regions = palette)