ここの記事はだいぶ古いです.日本地図を描く JapanPrefMap を見ましょう.

ひととおり描いてみる

    1. 次のライブラリをインストールし、読み込む。

library(gpclib)
library(maptools)
library(classInt)

  1. (手動)Rのメニュー:「ファイル」→「ディレクトリの変更」で、作業ディレクトリをjapank.zipを展開したフォルダに変更!
  2. ESRIのサイトから日本の市町村地図のシェープファイル: japan_ver62.zipをダウンロードし、作業ディレクトリに解凍する。
  3. shapeファイルの読み込み

    > d30 <-  readShapeSpatial(“./japan_ver62/japan_ver62.shp”)
    > gpclibPermit()

    [1] TRUE

  4. 島部分で日本地図が小さくなるので割愛する(沖縄の再描画をしたいのだが、まだやり方がわからない。

    > #北方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,]

  5. 市町村境界を都道府県境界に統合する

    > d32 <- unionSpatialPolygons(d31, d31$PREF)
    > # sapply(slot(d32, “polygons”), function(i) slot(i, “ID”))

  6. (手動)都道府県の表記を一致させたcsvデータファイルの作成(negiprd2009.csv)
  7. csvファイルの読み込み

    > d02<-read.csv(“negiprod2009.csv”,h=T)

  8. 出荷率を付け加える(ここは地図描画の手順とは直接関係ない)

    > 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)
    >

  9. データの確認

    > names(d02)
    [1] “pref” “area” “ypa” “yield” “shipments” “spy”
    > d02$pref
    [1] 北海道 青森県 岩手県 宮城県 秋田県 山形県 福島県 茨城県
    [9] 栃木県 群馬県 埼玉県 千葉県 東京都 神奈川県 新潟県 富山県
    [17] 石川県 福井県 山梨県 長野県 岐阜県 静岡県 愛知県 三重県
    [25] 滋賀県 京都府 大阪府 兵庫県 奈良県 和歌山県 鳥取県 島根県
    [33] 岡山県 広島県 山口県 徳島県 香川県 愛媛県 高知県 福岡県
    [41] 佐賀県 長崎県 熊本県 大分県 宮崎県 鹿児島県 沖縄県
    47 Levels: 愛知県 愛媛県 茨城県 岡山県 沖縄県 岩手県 岐阜県 … 和歌山県

  10. 都道府県境地図を描く関数 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)
    }
  11. 地図を描いてみる

    > mappref01(d32,d02[,c(1,6)],xlb=”出荷率”)

  12. 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)