R による地理情報のデータビジュアライゼーション (2)

[昨日](http://qiita.com/ynakayama/items/787fb19ff4e1fd00fb70)は地理情報まわりの基本的な知識を紹介しました。

## シェープファイルの入手

さて GIS を利用するにあたりシェープファイルを自作するのは現実的ではありませんから、どこからか入手することが必要です。

### 主なシェープファイルの入手先

国土交通省 GIS データのダウンロード
http://nrb-www.mlit.go.jp/kokjo/inspect/landclassification/download/index.html

ESRI ジャパン 全国市区町村界データ
http://www.esrij.com/products/data/japan-shp/

Global Administrative Areas
http://www.gadm.org/country

いずれにせよ使用許諾などに目を通し、お使いの用途で利用可能であることを確認してください。今回の例では一番下の Global Administrative Areas を利用します。

### シェープファイルを R で表示してみる

シェープファイルの Viewer はそれほどありません。 GIS の利用がなかなか進まないのはマトモな Viewer が無いからだとも言われているほどです。データ分析者は R を日常的に使っていますから、いっそのこと R そのものを Viewer として使っても良いでしょう。まずはシンプルにシェープファイルを閲覧してみます。

```r
library(gpclib) # maptools の前提パッケージ gpclib を R で使うのに必要
library(maptools)
library(RColorBrewer) # 統計グラフで使える便利なカラーパレット

# 経度と緯度を設定する
xlim <- c(128, 146)
ylim <- c(30, 46)

# まずはシンプルに JPN_adm1.shp を表示してみる
png("image.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
jpn <- readShapePoly("JPN_adm1.shp")
plot(jpn, xlim=xlim, ylim=ylim)
dev.off()
```

![image.png](https://qiita-image-store.s3.amazonaws.com/0/38501/b013addb-1783-ebe8-6cbe-1ab82cdf0f8a.png)

これが素のシェープファイルです。
次に、県の境界ごとにカラーパレットを割り当ててみます。

```r
col <- sample(1:8,size=47,replace=TRUE)
png("image2.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plot(jpn, xlim=xlim, ylim=ylim, col=brewer.pal(8,"Accent")[col])
dev.off()
```

![image2.png](https://qiita-image-store.s3.amazonaws.com/0/38501/890372bd-1e6d-bd11-af90-03e93a9eea99.png)

今度はシェープファイルの都心付近を拡大して、都市をポイントしてみます。

```r
jpn <- readShapePoly("JPN_adm2.shp")
png("image3.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plot(jpn[jpn$NAME_1=="Tokyo",])
points(139.7036, 35.69389, lwd=20, col="red")
text(139.7036, 35.69389, "新宿", col="blue", adj = c(-0.3,0.5), cex=2)
dev.off()
```

![image3.png](https://qiita-image-store.s3.amazonaws.com/0/38501/7e392f3e-5a9e-e1b6-31a1-71b6f44643d2.png)

R のインタプリタで表示を確認しながら、描画したいデータを絞り込んでいくと良いでしょう。

## オレゴン気候局のデータを使う

次に[昨日](http://qiita.com/ynakayama/items/787fb19ff4e1fd00fb70)も紹介した、オレゴン気候局のデータを利用してシェープファイル描画の流れを追ってみます。

Maps in R -- Examples
http://geog.uoregon.edu/geogr/topics/maps.htm

```r
# まずは依存するライブラリをロードする
library(gpclib)
library(maptools)
library(RColorBrewer)
library(classInt)

# シェープファイルを読み込む
# オレゴン州のアウトライン
orotl.shp <- readShapeLines('orotl.shp',
    proj4string=CRS("+proj=longlat"))
# オレゴン州の気候局データ
orstations.shp <- readShapePoints('orstations.shp',
    proj4string=CRS("+proj=longlat"))
# オレゴン州郡の国勢調査データ（ポリゴン）
orcounty.shp <- readShapePoly('orcounty.shp',
    proj4string=CRS("+proj=longlat"))
```

これら地理情報に対し、国勢調査データなどをマッピングしていきます。

```r
# 緯度経度を含む CSV ファイル
orstationc <- read.csv("orstationc.csv")
# オレゴン州郡国勢調査データ
orcountyp <- read.csv("orcountyp.csv")
# 国名を持つデータ·セット
cities <- read.csv("cities2.csv")  

# まずはデータの中身を表示してみる
summary(orcounty.shp)
attributes(orcounty.shp)
attributes(orcounty.shp@data)
attr(orcounty.shp,"polygons")
```

### 元のシェープファイルの中身を覗いてみる

```r
png("image0-1.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plot(orotl.shp, xlim=c(-124.5, -115), ylim=c(42,47))
dev.off()
```

![image0-1.png](https://qiita-image-store.s3.amazonaws.com/0/38501/eb30d372-efb2-1266-0bba-b64b41dc7749.png)

```r
png("image0-2.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
dev.off()
```

![image0-2.png](https://qiita-image-store.s3.amazonaws.com/0/38501/29be0d3e-a6cb-f87e-1a78-da3970b291c4.png)

```r
png("image0-3.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plot(orstations.shp, xlim=c(-124.5, -115), ylim=c(42,47))
dev.off()
```
![image0-3.png](https://qiita-image-store.s3.amazonaws.com/0/38501/8e1f3cc5-69d7-1c0d-0497-afc96cd75845.png)

### シンプルなマッピング

ここからマッピングしていきます。

#### オレゴン州郡の国勢調査データ - orcounty.shpシェイプファイル内の属性データ

```r
png("image1-1.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plotvar <- orcounty.shp@data$POP1990
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu") 
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
plot(orcounty.shp, col=colcode, add=T)
title(main="Population 1990",
    sub="Quantile (Equal-Frequency) Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image1-1.png](https://qiita-image-store.s3.amazonaws.com/0/38501/cbcbe107-4909-a034-071a-ea4279bb8b62.png)

#### オレゴン州の気候局データ - orstationc.csvファイルのデータ、orotl.shpにおけるベースマップ

```r
png("image1-2.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plotvar <- orstationc$tann
nclr <- 8
plotclr <- brewer.pal(nclr,"PuOr") 
plotclr <- plotclr[nclr:1] # reorder colors
class <- classIntervals(plotvar, nclr, style="equal")
colcode <- findColours(class, plotclr)
plot(orotl.shp, xlim=c(-124.5, -115), ylim=c(42,47))
points(orstationc$lon, orstationc$lat, pch=16, col=colcode, cex=2)
points(orstationc$lon, orstationc$lat, cex=2)
title("Oregon Climate Station Data -- Annual Temperature",
    sub="Equal-Interval Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image1-2.png](https://qiita-image-store.s3.amazonaws.com/0/38501/92241f83-5bb0-fd06-1458-dcea1613e851.png)

#### オレゴン州の気候局データ - 形状ファイル内の場所とデータ

```r
png("image1-3.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plotvar <- orstations.shp@data$pann
nclr <- 5
plotclr <- brewer.pal(nclr,"BuPu") 
class <- classIntervals(plotvar, nclr, style="fixed",
fixedBreaks=c(0,200,500,1000,2000,5000))
colcode <- findColours(class, plotclr)
orstations.pts <- orstations.shp@coords # get point data
plot(orotl.shp, xlim=c(-124.5, -115), ylim=c(42,47))
points(orstations.pts, pch=16, col=colcode, cex=2)
points(orstations.pts, cex=2)
title("Oregon Climate Station Data -- Annual Precipitation",
    sub="Fixed-Interval Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image1-3.png](https://qiita-image-store.s3.amazonaws.com/0/38501/7118731c-0e86-5bcf-8c9d-696db51697cb.png)

### カラースケールと表現のバリエーション

#### オレゴン州郡の国勢調査データ - 等しい周波数のクラス間隔

```r
png("image2-1.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plotvar <- orcounty.shp@data$AREA
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu")
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
plot(orcounty.shp, col=colcode, add=T)
title(main="Area",
sub="Quantile (Equal-Frequency) Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image2-1.png](https://qiita-image-store.s3.amazonaws.com/0/38501/0e76b3ce-1925-a234-6c99-0ef6d257e10f.png)

#### オレゴン州郡の国勢調査データ - 等幅クラス毎

```r
png("image2-2.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
plotvar <- orcounty.shp@data$AREA
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu")
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
class <- classIntervals(plotvar, nclr, style="equal")
colcode <- findColours(class, plotclr)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
plot(orcounty.shp, col=colcode, add=T)
title(main="Area",
    sub=" Equal-Width Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
#equal-width class intervals of 1990 population
plotvar <- orcounty.shp@data$POP1990
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu") 
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
class <- classIntervals(plotvar, nclr, style="equal")
colcode <- findColours(class, plotclr)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
plot(orcounty.shp, col=colcode, add=T)
title(main="Population 1990",
    sub=" Equal-Width Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image2-2.png](https://qiita-image-store.s3.amazonaws.com/0/38501/38eacdc2-887c-7313-3626-a4132fe81fcd.png)

#### オレゴン州郡の国勢調査データ - バブルプロット

```r
png("image2-3.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
# bubble plot equal-frequency class intervals
plotvar <- orcounty.shp@data$AREA  
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu") 
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
max.symbol.size=12
min.symbol.size=1
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
symbol.size <- ((plotvar-min(plotvar))/
    (max(plotvar)-min(plotvar))*(max.symbol.size-min.symbol.size)
    +min.symbol.size)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
orcounty.cntr <- coordinates(orcounty.shp)
points(orcounty.cntr, pch=16, col=colcode, cex=symbol.size)
points(orcounty.cntr, cex=symbol.size)
    text(-120, 46.5, "Area: Equal-Frequency Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
# bubble plot equal-frequency class intervals
plotvar <- orcounty.shp@data$POP1990
nclr <- 8
plotclr <- brewer.pal(nclr,"BuPu") 
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
max.symbol.size=12
min.symbol.size=1
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
symbol.size <- ((plotvar-min(plotvar))/
(max(plotvar)-min(plotvar))*(max.symbol.size-min.symbol.size)
+min.symbol.size)
plot(orcounty.shp, xlim=c(-124.5, -115), ylim=c(42,47))
orcounty.cntr <- coordinates(orcounty.shp)
points(orcounty.cntr, pch=16, col=colcode, cex=symbol.size)
points(orcounty.cntr, cex=symbol.size)
text(-120, 46.5, "Population: Equal-Frequency Class Intervals")
legend(-117, 44, legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")
dev.off()
```

![image2-3.png](https://qiita-image-store.s3.amazonaws.com/0/38501/9b85ea48-d720-4a64-2cad-6fabcd186355.png)

#### オレゴン州郡の国勢調査データ - （擬似）ドット密度マップ

```r
png("image2-4.png", width = 480, height = 480, pointsize = 12, bg = "white", res = NA)
# maptools dot-density maps
# warning: this can take a little while
plottitle <- "Population 1990, each dot=1000"
orpolys <- SpatialPolygonsDataFrame(orcounty.shp, data=as(orcounty.shp, "data.frame"))
plotvar <- orpolys@data$POP1990/1000.0
dots.rand <- dotsInPolys(orpolys, as.integer(plotvar), f="random")
plot(orpolys, xlim=c(-124.5, -115), ylim=c(42,47))
plot(dots.rand, add=T, pch=19, cex=0.5, col="magenta")
plot(orpolys, add=T)
title(plottitle)
dots.reg <- dotsInPolys(orpolys, as.integer(plotvar), f="regular")
plot(orpolys, xlim=c(-124.5, -115), ylim=c(42,47))
plot(dots.reg, add=T, pch=19, cex=0.5, col="purple")
plot(orpolys, add=T)
title(plottitle)
dev.off()
```

![image2-4.png](https://qiita-image-store.s3.amazonaws.com/0/38501/7f06eaae-6997-aa40-cf70-4c3bf89e6b67.png)

## まとめ

シェープファイルで地理情報を描画し、その上にデータを散りばめることができることがわかりました。

## 参考

Rの基本グラフィックス機能またはggplot2を使って地図を描くには
http://sudillap.hatenablog.com/entry/2013/03/26/210202

Maps in R -- Examples
http://geog.uoregon.edu/geogr/topics/maps.htm