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>