ここの記事はだいぶ古いです.日本地図を描く JapanPrefMap を見ましょう.
ひととおり描いてみる
-
- 次のライブラリをインストールし、読み込む。
library(gpclib)
library(maptools)
library(classInt)
- (手動)Rのメニュー:「ファイル」→「ディレクトリの変更」で、作業ディレクトリをjapank.zipを展開したフォルダに変更!
- ESRIのサイトから日本の市町村地図のシェープファイル: japan_ver62.zipをダウンロードし、作業ディレクトリに解凍する。
- shapeファイルの読み込み
> d30 <- readShapeSpatial(“./japan_ver62/japan_ver62.shp”)
> gpclibPermit()
[1] TRUE - 島部分で日本地図が小さくなるので割愛する(沖縄の再描画をしたいのだが、まだやり方がわからない。
> #北方4島を除く
> d31<-d30[d30\$JCODE>1000,]
> d31<-d31[d31\$JCODE<1700|d31\$JCODE>2000,]
> #八丈島とかを除く
> d31<-d31[d31\$JCODE<13310|d31\$JCODE>14000,]
> #沖縄を除く
> d31<-d31[floor(d31\$JCODE/1000)!=47,]
> #鹿児島の島部分を削除
> d31<-d31[d31$JCODE<46500,]
> d31<-d31[d31\$JCODE!=46222,]
> d31<-d31[d31$JCODE!=46304,] - 市町村境界を都道府県境界に統合する
> d32 <- unionSpatialPolygons(d31, d31$PREF)
> # sapply(slot(d32, “polygons”), function(i) slot(i, “ID”)) - (手動)都道府県の表記を一致させたcsvデータファイルの作成(negiprd2009.csv)
- csvファイルの読み込み
> d02<-read.csv(“negiprod2009.csv”,h=T)
- 出荷率を付け加える(ここは地図描画の手順とは直接関係ない)
> d02<-cbind(d02,spy=d02[,5]/d02[,4])
> d03<-d02[,c(2,6)]
> plot(d03,type=”n”,main=”秋冬ネギの作付面積と出荷率の相関”,ylim=c(0,1))
> text(d03,labels=d02[1:46,1],cex=0.7)
> - データの確認
> names(d02)
[1] “pref” “area” “ypa” “yield” “shipments” “spy”
> d02$pref
[1] 北海道 青森県 岩手県 宮城県 秋田県 山形県 福島県 茨城県
[9] 栃木県 群馬県 埼玉県 千葉県 東京都 神奈川県 新潟県 富山県
[17] 石川県 福井県 山梨県 長野県 岐阜県 静岡県 愛知県 三重県
[25] 滋賀県 京都府 大阪府 兵庫県 奈良県 和歌山県 鳥取県 島根県
[33] 岡山県 広島県 山口県 徳島県 香川県 愛媛県 高知県 福岡県
[41] 佐賀県 長崎県 熊本県 大分県 宮崎県 鹿児島県 沖縄県
47 Levels: 愛知県 愛媛県 茨城県 岡山県 沖縄県 岩手県 岐阜県 … 和歌山県 - 都道府県境地図を描く関数 mappref01 の定義
mappref01<-function(mdata,ddata,drwNo=c(1:47),nbr=6,xlb=colnames(ddata)[2],brcol=c("white", "red2")){ #mdata:japank.shpをreadShapePolyで読み込んだもの #ddata:描画データ。1列目はrow.names(mdata)と同じ表記の都道府県名、2列めは描画データ #drwNo:描画する都道府県名。row.names(mdata)の番号。 #nbr:階級区分数-1 #xlb:グラフラベル #brcol: 描画色#都道府県名を一致させる d32n<-match(row.names(mdata),ddata$pref) d03<-ddata[d32n,] row.names(d03)<-row.names(mdata) d04<-SpatialPolygonsDataFrame(mdata,d03) #出荷量の分布地図をつくる #ここでは、n=2,3,4,5,...の等間隔に分けるが、 #区切りをキリの良いところにして、 #しかも、0付近だけは別に小さく区切る(他の1階級の1/5幅)という #n+1階級を求める関数をつくってみる。 #ただし、x>0の場合のみ fbr01<-function(x,n){ mx01<-max(x,na.rm=T) vc01<-ceiling((mx01/10^floor(log10(mx01))*10)/n)*10^(floor(log10(mx01))-1) br01<-c(0,vc01/5,vc01) for(i in 2:n) br01<-c(br01,vc01*i) br01 }#区分点を作成する brp01<<-fbr01(d03[,2],nbr)#作図に含まない都道府県を除外 #この後の処理で、NAが問題にされるようなレコーはここで除外。 d06<-d03[drwNo,] d06<-na.omit(d06) d05<-d04[row.names(d06),] #階級分け nbr1<-nbr+1 d23<-classIntervals(d06[,2], n=nbr1,style="fixed", fixedBreaks=brp01)#色の作成 d24 <- findColours(d23,brcol) #作図 plot(d05,col=d24) #タイトルと凡例 title(xlb, cex=0.5) legend("topleft", fill=attr(d24,"palette"), cex=0.6,legend=names(attr(d24,"table")), bty=nbr1) }
- 地図を描いてみる
> mappref01(d32,d02[,c(1,6)],xlb=”出荷率”)
- 3つの地図を並べる
> par(mfrow=c(2,2))
> par(mai=c(0,0,0.6,0))
> mappref01(d32,d02[,c(1,2)],xlb=”作付け面積”)
> mappref01(d32,d02[,c(1,5)],xlb=”出荷量”)
> mappref01(d32,d02[,c(1,6)],xlb=”出荷率”)
> par(oma=c(3,0,0,0))
> mtext(“ネギ H21年(作物統計)”,outer=T,side=1)
ソースコード
#ライブラリの読み込み
library(gpclib)
library(maptools)
library(classInt)
##################
#全国市町村shapeファイル
#
#ESRIのHPから日本市町村地図にshapeファイルを入手する
#http://www.esrij.com/products/gis_data/japanshp/japanshp.html
# japan_ver62.zipをDL。
#作業ディレクトリに解凍する。
#maptool用のデータに変換
d30 <- readShapeSpatial("./japan_ver62/japan_ver62.shp")
gpclibPermit()
##################
# 都道府県地図を描く
#
####
#島部分で日本地図が小さくなるので割愛する
#沖縄の再描画をしたいが、まだできない。
#
#北方4島を除く
d31<-d30[d30$JCODE>1000,]
d31<-d31[d31$JCODE<1700|d31$JCODE>2000,]
#八丈島とかを除く
d31<-d31[d31$JCODE<13310|d31$JCODE>14000,]
#沖縄を除く
d31<-d31[floor(d31$JCODE/1000)!=47,]
#鹿児島の島部分を削除
d31<-d31[d31$JCODE<46500,]
d31<-d31[d31$JCODE!=46222,]
d31<-d31[d31$JCODE!=46304,]
#沖縄を別枠で表示したいのだが、やり方がわからない・・・
#????
#都道府県に統合する
d32 <- unionSpatialPolygons(d31, d31$PREF)
# sapply(slot(d32, "polygons"), function(i) slot(i, "ID"))
#(手動)都道府県の表記を一致させたcsvデータファイルの作成(negiprd2009.csv)
#csvファイルの読み込み
d02<-read.csv("negiprod2009.csv",h=T)
#出荷率を付け加える(ここは地図描画の手順とは直接関係ない)
d02<-cbind(d02,spy=d02[,5]/d02[,4])
d03<-d02[,c(2,6)]
plot(d03,type="n",main="秋冬ネギの作付面積と出荷率の相関",ylim=c(0,1))
text(d03,labels=d02[1:46,1],cex=0.7)
#データの確認
names(d02)
d02$pref
###########################################
#都道府県境地図を描く関数 mappref01 の定義
#ESRIの全国市町村境界シェープファイルを読み
#込んだもの(を加工したもの)を使って描く。
###########################################
mappref01<-function(mdata,ddata,drwNo=c(1:47),nbr=6,xlb=colnames(ddata)[2],brcol=c("white", "red2")){
#mdata:japank.shpをreadShapePolyで読み込んだもの
#ddata:描画データ。1列目はrow.names(mdata)と同じ表記の都道府県名、2列めは描画データ
#drwNo:描画する都道府県名。row.names(mdata)の番号。
#nbr:階級区分数-1
#xlb:グラフラベル
#brcol: 描画色
#都道府県名を一致させる
d32n<-match(row.names(mdata),ddata$pref)
d03<-ddata[d32n,]
row.names(d03)<-row.names(mdata)
d04<-SpatialPolygonsDataFrame(mdata,d03) #出荷量の分布地図をつくる #ここでは、n=2,3,4,5,...の等間隔に分けるが、 #区切りをキリの良いところにして、 #しかも、0付近だけは別に小さく区切る(他の1階級の1/5幅)という #n+1階級を求める関数をつくってみる。 #ただし、x>0の場合のみ
fbr01<-function(x,n){
mx01<-max(x,na.rm=T)
vc01<-ceiling((mx01/10^floor(log10(mx01))*10)/n)*10^(floor(log10(mx01))-1)
br01<-c(0,vc01/5,vc01)
for(i in 2:n) br01<-c(br01,vc01*i)
br01
}
#区分点を作成する
brp01<<-fbr01(d03[,2],nbr)
#作図に含まない都道府県を除外
#この後の処理で、NAが問題にされるようなレコーはここで除外。
d06<-d03[drwNo,]
d06<-na.omit(d06)
d05<-d04[row.names(d06),]
#階級分け
nbr1<-nbr+1
d23<-classIntervals(d06[,2], n=nbr1,style="fixed", fixedBreaks=brp01)
#色の作成
d24 <- findColours(d23,brcol)
#作図
plot(d05,col=d24)
#タイトルと凡例
title(xlb, cex=0.5)
legend("topleft",
fill=attr(d24,"palette"),
cex=0.6,legend=names(attr(d24,"table")),
bty=nbr1)
}
#地図を描いてみる
mappref01(d32,d02[,c(1,6)],xlb="出荷率")
#3つの地図を描く
par(mfrow=c(2,2))
par(mai=c(0,0,0.6,0))
mappref01(d32,d02[,c(1,2)],xlb="作付け面積")
mappref01(d32,d02[,c(1,5)],xlb="出荷量")
mappref01(d32,d02[,c(1,6)],xlb="出荷率")
par(oma=c(3,0,0,0))
mtext("ネギ H21年(作物統計)",outer=T,side=1)