ある地理上での感染伝播的なものをやろうとか思ってみた。

MikuHatsune2011-02-11

感染症学を学んでいる。
本日の講義は、感染症の国際情勢と国内対応についての、お話し的な講義だった。
日本国内で、ある感染症が発生したときの患者数を地図上にプロットしたいな、とか思ったので。

 
Rで地図の項があるけど、Rmapというものはもう使っていないらしい。こちらのページは使いにくい。
こちらでRを使った地図描出、およびこちらの派生リンクで地図を描いている。今回はこれらを参考にさせていただいた。
というわけで、描く手順の覚書き。
 
spsurveyというパッケージをRにインストールする。
 
spsurveyを使って地図を描くには、shp(シェイプ)というファイルが必要らしい。これを取ってくるのだが、無償で提供して下さっているサイトがあるらしい。ひとつがESRIジャパン株式会社で、ここから取ったデータは図の左側である。
都道府県から、市町村、人口、世帯数などが同梱されている模様。
今回は都道府県単位で描きたいので、Global Administrative Areaからちょっと違うshpファイルをもらう。これで図の右側が描ける。
ファイルのダウンロードだが、上のDownloadから
Country:Wanted Area
File format:Shape file
を選んでOKすれば、任意の国を(たぶん州とか都道府県単位で)描ける。
ダウンロードしたらおそらくzipファイルなので、解凍して、Rのワーキングディレクトリにおいておくのが簡単だと思う。
ちなみにESRIジャパン株式会社のほうは、解凍したらさらにもうひとつディレクトリがあるという面倒くささで、余計な外箱からファイルを出してRのワーキングディレクトリに置く。じゃないと以下のスクリプトが走りません。
 
これで準備完了。さっそく描いてみよう。
 

library(spsurvey)

これをしないと始まりません。
地図座標データっぽいものを読み込みます。

jpn    <- read.shape("japan_ver62/japan_ver62.shp")
jpn.GA <- read.shape("JPN.adm/JPN.adm.shp")

上が左図、下が右図を描くデータ。
今のRワーキングディレクトリは

getwd()
"C:/Users/MikuHatsune/Documents/R"

だと思うので、read.shapeの読み込みのときにディレクトリ指定を付加しておく。
ちなみにけっこう時間かかる。
あとは

plot(jpn)
plot(jpn.GA)

で日本地図は描ける。
 
いろいろいじくろうと思ったら
plotの範囲は経度(xlim)・緯度(ylim)で指定できる。

jpn.GA[[5]]

都道府県が収まっているから、これの順番を維持して色づけしたりできる。
 
色付けにはカラーパレットを用いた。
NK的には、拡散と組み合わせても楽しいんじゃないか、との意見。
 
スクリプト

library(spsurvey)
jpn    <- read.shape("japan_ver62/japan_ver62.shp")	# 15分かかった。
jpn.GA <- read.shape("JPN.adm/JPN.adm.shp")             #  5分かかった。
xlim<- c(128, 146)                                      # 経度
ylim<- c(30, 45)                                        # 緯度
par(mfrow=c(1,2))
flag<- sample(0:1,size=47,replace=TRUE)                 # 色づけされる都道府県を適当に選ぶ。
plot(jpn,xlim=xlim,ylim=ylim,col=flag*5)
points(135.75,35,lwd=5,col="red")
# もっとうまい選び方ないかな?
Pandemic<- c(which(jpn.GA[[5]]=="Kyoto"),
             which(jpn.GA[[5]]=="Osaka"),
             which(jpn.GA[[5]]=="Shiga"),
             which(jpn.GA[[5]]=="Nara"),
             which(jpn.GA[[5]]=="Mie"),
             which(jpn.GA[[5]]=="Wakayama"),
             which(jpn.GA[[5]]=="Hyogo"),
             which(jpn.GA[[5]]=="Okayama"),
             which(jpn.GA[[5]]=="Tottori"),
             which(jpn.GA[[5]]=="Shimane"),
             which(jpn.GA[[5]]=="Hiroshima"),
             which(jpn.GA[[5]]=="Fukui"),
             which(jpn.GA[[5]]=="Gifu"),
             which(jpn.GA[[5]]=="Aichi"),
             which(jpn.GA[[5]]=="Nagano"),
             which(jpn.GA[[5]]=="Niigata"),
             which(jpn.GA[[5]]=="Ishikawa"),
             which(jpn.GA[[5]]=="Shizuoka"),
             which(jpn.GA[[5]]=="Toyama"))
flag[Pandemic]<- c(rep(3,5),rep(2,6),rep(1,8))
# 赤いところほどヤバい、くらいの意味合い。
col.pallet<- rev(heat.colors(max(flag)+1))
col.pref<- rep(0,47)
  for(i in 0:max(flag)){
    col.pref[which(flag==i)]<- col.pallet[i+1]
  }
plot(jpn.GA,xlim=xlim,ylim=ylim,col=col.pref)
points(135.75,35,lwd=5,col=3)
title("Infected patients")