Rでネパールの地図の塗り分けをする

Dr. Fukuyo
Professor, Yamaguchi University
Last updated: April 6, 2018

以前，Rのパッケージ"maptools"とDIVA-GISのデータを使って，カンボジアの地図を塗り分けてみたわけです（繰り返しになりますが塗り分け地図のことは英語でコロプレス・マップ(choropleth map)と言います）。

今回，およそ6年半ぶりにコロプレス・マップに関する記事を書きます。今回はネパールです。

準備作業（以前と同じ）

地図の塗り分け作業をする前に，余分なデータを消去しておく。Rを立ち上げて，以下のコマンドを打ち込む。

ls()
rm(list=ls(all=TRUE))

これで今までのデータセットは全てクリアされる。

あと，以前の記事を読んで，DIVA-GISからDistrict形状データ"NPL_adm3.shp"他一式を入手し，適当なフォルダに保存しておくこと。描画のために最低限必要なファイルは次の通りである：


ちなみに，上記のファイルのうち，"NPL_adm3.csv"と"NPL_adm3.dbf"に関しては著者が世帯数，人口，面積のデータを書き加えたものを用意しているので，DIVA-GISのデータそのものではなく，上述のリンク先から世帯数，人口，面積のデータを加筆した後のファイルをダウンロードしていただきたい。

つぎに，以下のコマンドを打ち込んでmaptoolsを呼びだす。

library(maptools)

人口密度によるネパールの塗り分け地図の作成

シェープファイルを読み込み，district毎に人口密度に応じて色付けをする。色付けはheat color（レモン色から段階的に黄色，オレンジ，赤，茶色と変化するカラーパレット）で行う。

スクリプトは以下の通り。

なお，下線の部分は州境形状データや水域データを保存したフォルダに遭わせて書きなおすこと。


x <- readShapePoly("c:\\rtemp\\NPL_adm\\NPL_adm3.shp")
				#シェープファイルの読み込み

ix <- length(x)			# データのサイズ（行＝districtの数）
zdata <- x$POP_2011/x$AREA	# 人口密度の計算；POP_2011はdistrict別人口，AREAはdistrict面積
colnam <- rep("white", length=ix)
				# 塗分けのための変数
t <- c(50, 100, 500, 1000)	# 塗分けのための閾値
color.set <- c("lemonchiffon", "yellow", "orange", "red", "darkred")
				# カラーパレット(heat color)
unit <- "Persons/km^2" # zdataの単位

# zdata（人口密度）に応じて，カラーパレットから色を選んで
# 塗分けのための変数colnamに代入する
for (ii in 1:ix) {
	if (zdata[ii] > t[4]) {
		colnam[ii] <- color.set[5]
	} else if (zdata[ii] > t[3]) {
		colnam[ii] <- color.set[4]
	} else if (zdata[ii] > t[2]) {
		colnam[ii] <- color.set[3]
	} else if (zdata[ii] > t[1]) {
		colnam[ii] <- color.set[2]
	} else {
		colnam[ii] <- color.set[1]
	}
}


plot(x, col=colnam) # 人口密度で色分け



# 以下は凡例の描画用
oo <- par("usr")	# プロット範囲（経度・緯度）をooに保存
ox <- oo[1]			# 経度の値の小さい方を原点x座標に
oy <- (oo[3] + oo[4])/2	# 緯度の値の小さい方と大きいほうの中間を原点y座標に
oxl <- oo[2] - oo[1]		#
oyl <- oo[4] - oo[3]		#
par(new = T)			# 現在の図に新たな図を重ね書きする
for (ii in 1:5) {
	aa <- c(ox + oxl * 15 / 20,
		   ox + oxl * 15 / 20,
		   ox + oxl * 16 / 20,
		   ox + oxl * 16 / 20)
	bb <- c(oy + oyl / 20 * ii,
		   oy + oyl / 20 * (ii + 1),
		   oy + oyl / 20 * (ii + 1),
		   oy + oyl / 20 * ii)
	polygon(aa, bb, col=color.set[ii])
	if (ii < 5) text(ox + oxl * 16.5 / 20, oy + oyl / 20 * (ii + 1), labels = t[ii], adj = 0)
	else text(ox + oxl * 15 / 20, oy + oyl / 20 * (ii + 1.5), labels = unit, adj = 0)
}

実際に使ってみる


次のような画像が現れればOK。
Population Density
参考：
<back>