放射能の測定結果の表し方にはお約束があって、検出下限値未満だったときは、検出下限値の数値の左側に不等号(小なり記号)をつけて表現します。
つまり
8.6
と書いてあれば検出値が8.6だったということですが、
<8.6
と書いてあれば、検出下限値が8.6という環境(測定器機の精度、試料の量、計測時間 etc.)で測定したが、値が小さいせいで検出されなかった、ということになります。
これを知らない人がやらかしてしまったのが、↓こちらの誤報だったりします。
新潟県:「女性セブン」誌の放射能関連記事に誤った内容が掲載されたため、訂正と謝罪を求める抗議文を送付しました。
「出典は厚生労働省のホームページであり、不等号(<)が付されている値が検出下限値であることを理解せず転記した」 by小学館
とのこと。記者一人が勘違いしたとしても、チェック機構はなかったんでしょうかね。不等号を見て、「これ何だろう?」って思わなかったんでしょうか。
で、話を戻しまして・・・
このような表記法となっておりますので、データはこんな感じになっていたりします。
> data
食品 結果
1 A <7.8
2 B 8.9
3 C <6.7
4 D 9.0
下限値未満をどう扱うかは分析の方針次第ですが、仮に「下限値未満はゼロとして扱う」としたかったとしましょう。その場合、不等号付きの結果を加工してやる必要があります。
また、data$結果はcharacter型かfactor型(因子型)になっているので、このままでは平均を出したり、ヒストグラムを書いたりはできません。
「結果」をnumeric型に変換した「結果2」という列を追加してみましょう。
for ( i in 1:nrow(data) ) {
if ( substr(data$結果[i], 1, 1) == "<" ) {
# 不等号付きのときは、ゼロにする
data$結果2[i] <- 0
} else {
# 不等号付きでないときは、そのまま数値にする
data$結果2[i] <- as.numeric(data$結果[i])
}
}
>data
食品 結果 結果2
1 A <7.8 0.0
2 B 8.9 8.9
3 C <6.7 0.0
4 D 9.0 9.0
結果2はnumeric型になっていますので、mean(data$結果2) とか、hist(data$結果2)とかすることも可能です。
が、このようにfor文とif文で処理する人はあまりいないと思います。Rにはifelseという便利な記法があり、しかもifelseは「ベクトル化された演算」(詳しくは「アート・オブ・Rプログラミング」参照)であるため、for文で記述するよりも高速に実行できるというメリットもあります。
先ほどのfor文によるループをifelseで書き直すと以下のようになります。
data$結果2 <- ifelse( substr(data$結果, 1, 1) == "<", 0, as.numeric(data$結果) )
ifelseは第1引数を評価して、真ならば第2引数を、偽ならば第3引数を返すという動きをします。引数としてベクトルを渡せば、すべての要素について処理を繰り返してくれます。
実行結果として生成される、データフレームの列は、先ほどのfor文の場合と同じになります。
が、しかし、上記の書き方だと、↓こんなwarningが出てしまうんですよね。
警告メッセージ:
In ifelse(substr(data$結果, 1, 1) == "<", 0, as.numeric(data$結果)) :
強制変換により NA が生成されました
結果の中にNAなどないし、強制変換しなきゃいけないような処理もしていないはずだし・・・、なぜでしょうか?
どうやらこれは、ifelseが「第1引数の真偽にかかわらず、第2引数・第3引数を評価する」という動きをしているのことに起因しているようです。分岐はせずにすべて評価しておいて、不要な結果は捨てるという動きをしているのでしょう。
つまり「<7.8」に対しても、as.numericしてしまい結果がNAとなる、でも第1引数の条件式は偽になるので、その結果は使われない。なので、最終結果としては問題なしです。
でも、警告が出るのが気持ち悪い。なんとかならないものか・・・
で、結論としては、下記のように書けばOKということが分かりました。(前置きが長くなりました・・・)
data$結果2 <- as.numeric( ifelse( substr(data$結果, 1, 1) == "<", 0, data$結果 ) )
最後にas.numericすればいいだけですね。0に対してもas.numericすることになってしまいますが、これについては警告は出ないので、こちらの方がベターかなと思います。
【重要な補足】
「結果」の列がcharacter型ではなく、factor型になっている場合があります。いやむしろデフォルトでCSVファイルを読み込むとそうなるので、そちらのケースが多いかもしれません。
その場合は・・・
data$結果2 <- as.numeric( ifelse( substr(data$結果, 1, 1) == "<", 0, as.character(data$結果) ) )
としてやる必要があります。
ファクタをそのまま as.numericすると文字列の表す数字ではなく、水準の値が返ってきてしまうんですよね。
水準というのは、"男"、"女"、"その他" というファクタがあったときに、内部的にはそれぞれに 1、2、3が割り当たっていたりするんですが、その数字のことです。
ファクタの水準値を使われてしまうと、数値が明らかにおかしくなるので、変換後のデータフレームの内容をちゃんとチェックしてから、次の処理を行いましょう。
2014年7月29日火曜日
2014年7月23日水曜日
Rでshapeを出力するときはPDFファイルにすると便利
表題は正確に言うと、shapeファイルを読み込んで、加工やいろいろな可視化を施した後、画像(的なもの)として出力するときは、PDFにすると便利ですよ、という感じの趣旨です。
shapeファイルを読み込んで、plotする方法は↓以前に書きました。
e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ
ただ、shapeファイルってものによっては、ものすごく大量のポリゴンを含んでいますよね。地名などを重畳した場合、文字サイズが大きいと互いに重なってしまう、文字サイズが小さいと読めなくなってしまう、なんて問題が起こります。
前述のエントリでは、高解像度のpngに出力するという対処法を書いたんですけど、高解像度とはいえラスタである以上、拡大すれば文字はギザギザになるし、解像度を大きくするのにも限度がある。
使っている環境(メモリとか? CPUのビット数とか?)によるんでしょうけど、あまり大きくすると、
> png("output.png", width=10000, height=10000)
以下にエラー png("output.png", width = 10000, height = 10000) :
デバイス png() を開始できませんでした
追加情報: 警告メッセージ:
1: In png("output.png", width = 10000, height = 10000) :
ビットマップを割り当てられません
2: In png("output.png", width = 10000, height = 10000) :
opening device failed
みたいな感じにエラーになってしまいます。
で、結論としてはPDFに出力するといい感じじゃん、ってことになりました。
コードは↓こんな感じ。
# GIS関係のパッケージ
library(maptools)
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")
# シェープファイルの読み込み
x <- readShapePoly("h22ka14137.shp", proj4string=pj)
# PDFファイルをオープン
pdf("output.pdf", family="Japan1GothicBBB")
# プロット(出力先はPDFファイルとなる)
plot(x, lwd=0.1)
text(coordinates(x), as.character(x@data$MOJI), cex=0.1)
# ファイルを閉じる
dev.off()
ちなみにPDFファイルに出力するときは、フォントを指定(family="Japan1GothicBBB"の箇所)しないと、文字化けが起きるというワナがあります。
出力されたPDFファイルは↓こんな感じ。
川崎市麻生区の町丁目境界と町丁目名称(PDFファイル)
↓全体として見ると、「字ちっちゃっ!」って感じですが。
こまごまとしたところでも文字が重ならないようにするにはこれくらいのサイズにする必要がありました。text関数のcexオプションで調整できます。
線もデフォルトだと拡大したときに太く感じられたので、細めにしました。plot関数のlwdオプションで調整できます。
境界も地名も画像ではなく、ベクタの情報として埋め込まれているので、拡大しても粗くなりません↓
「上麻生〇丁目」に混ざって、「上麻生」というのが点在するのは、たぶんこれらの地域に対して住居表示が実施されていないからですね。こんなにちっちゃく残すくらいなら、全部住居表示にしちゃえばよかったのに、って気もします。何か事情があったのかもしれませんが。■
[2014年9月10日追記]
PDFファイルのサイズが大きくなりすぎた場合は、↓こんな方法もお試しください。
SkyPDFを使ってPDFファイルのサイズを圧縮する: 主張
■
shapeファイルを読み込んで、plotする方法は↓以前に書きました。
e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ
ただ、shapeファイルってものによっては、ものすごく大量のポリゴンを含んでいますよね。地名などを重畳した場合、文字サイズが大きいと互いに重なってしまう、文字サイズが小さいと読めなくなってしまう、なんて問題が起こります。
前述のエントリでは、高解像度のpngに出力するという対処法を書いたんですけど、高解像度とはいえラスタである以上、拡大すれば文字はギザギザになるし、解像度を大きくするのにも限度がある。
使っている環境(メモリとか? CPUのビット数とか?)によるんでしょうけど、あまり大きくすると、
> png("output.png", width=10000, height=10000)
以下にエラー png("output.png", width = 10000, height = 10000) :
デバイス png() を開始できませんでした
追加情報: 警告メッセージ:
1: In png("output.png", width = 10000, height = 10000) :
ビットマップを割り当てられません
2: In png("output.png", width = 10000, height = 10000) :
opening device failed
みたいな感じにエラーになってしまいます。
で、結論としてはPDFに出力するといい感じじゃん、ってことになりました。
コードは↓こんな感じ。
# GIS関係のパッケージ
library(maptools)
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")
# シェープファイルの読み込み
x <- readShapePoly("h22ka14137.shp", proj4string=pj)
# PDFファイルをオープン
pdf("output.pdf", family="Japan1GothicBBB")
# プロット(出力先はPDFファイルとなる)
plot(x, lwd=0.1)
text(coordinates(x), as.character(x@data$MOJI), cex=0.1)
# ファイルを閉じる
dev.off()
ちなみにPDFファイルに出力するときは、フォントを指定(family="Japan1GothicBBB"の箇所)しないと、文字化けが起きるというワナがあります。
出力されたPDFファイルは↓こんな感じ。
川崎市麻生区の町丁目境界と町丁目名称(PDFファイル)
↓全体として見ると、「字ちっちゃっ!」って感じですが。
麻生区全体のPDFファイル |
こまごまとしたところでも文字が重ならないようにするにはこれくらいのサイズにする必要がありました。text関数のcexオプションで調整できます。
線もデフォルトだと拡大したときに太く感じられたので、細めにしました。plot関数のlwdオプションで調整できます。
境界も地名も画像ではなく、ベクタの情報として埋め込まれているので、拡大しても粗くなりません↓
PDFファイルなので、拡大しても線や文字は滑らかなまま |
「上麻生〇丁目」に混ざって、「上麻生」というのが点在するのは、たぶんこれらの地域に対して住居表示が実施されていないからですね。こんなにちっちゃく残すくらいなら、全部住居表示にしちゃえばよかったのに、って気もします。何か事情があったのかもしれませんが。■
[2014年9月10日追記]
PDFファイルのサイズが大きくなりすぎた場合は、↓こんな方法もお試しください。
SkyPDFを使ってPDFファイルのサイズを圧縮する: 主張
■
2014年7月22日火曜日
エクセルの日付(シリアル値)を、Rで使えるように変換する
例えば、セルの1つに「2014年7月22日」と入力してあったとしましょう。
これをXLConnectパッケージを使ってRで読み込むと、
> library(XLConnect)
> wb <- loadWorkbook("input.xlsx")
> ws <- readWorksheet(wb, sheet=1, header=F)
> ws
Col1
1 41842
のように、「2014年7月22日」と入力してあったセルの値が、「41842」となっています。
これはシリアル値といって、1900年1月1日を「1」、1月2日を「2」、・・・として日付を数値化したもので、エクセルは内部的にはこのように日付データを保持しています。というのを知っている人は多いと思います。
じゃあ、この値を日付としてRから使いたいときにどうやるのか? というのが今回のエントリの目的です。
Rには日付のシリアル表現のための関数があって、整数値と起源日付を指定することによって、シリアル値を日付(Date)型に変換することができます。
ということは、↓これでいけるはずだが・・・
> as.Date(41842, origin="1900-01-01") # 2014年7月22日のつもり
[1] "2014-07-24"
あれ? 2日ほどずれている。
このずれの原因の1つ目は以下のコードを実行することで、すぐに理解できます。
> as.Date(0, origin="1900-01-01")
[1] "1900-01-01"
> as.Date(1, origin="1900-01-01")
[1] "1900-01-02"
> as.Date(2, origin="1900-01-01")
[1] "1900-01-03"
エクセルでは1900年1月1日を「1」としているのに対し、Rのas.Date関数はoriginオプションで指定された日付(つまり1900年1月1日)を「0」としているわけですね。0オリジン vs 1オリジン のゴタゴタがこんなところにまで・・・
そして、ずれの原因の2つ目は、エクセルでシリアル値を実際に表示してみると判明します。
エクセル上でシリアル値を確認したい場合は、日付が入っているセルの書式設定で、表示形式タブを選択、ユーザー定義で「G/標準」を指定すればOKです。
絶対的な真理として、1900年はうるう年ではありません。つまり、1900年2月29日は存在せず、2月28日が59なら、3月1日を60にするべきなのです。なのに、そうなっていないせいで、1900年3月1日以降(あなたが扱うであろうほとんどの日付)に関しては、このずれが発生してしまうことになります。
このバグとも思えるような変な仕様は、(当時有力だった)「Lotus1-2-3」との互換性を維持するための意図的なものであったとのこと。
なんて、いきさつのことは知らなくても、やり方だけ覚えておけばOK。
【A】origin="1900-01-01"と指定するなら、シリアル値から2を引く または
というのでもいいけど、↓こっちの方がシンプルかも
【B】origin="1899-12-30"と指定する
【B】の方は、マニュアルにも書いてあって(?as.Dateと打てば見られます)、
## So for dates (post-1901) from Windows Excel
as.Date(35981, origin = "1899-12-30") # 1998-07-05
というサンプルが載っています。
マッキントッシュは起源日付が違うので、マック版エクセルでは、以下のようにすればいいらしいです。
## and Mac Excel
as.Date(34519, origin = "1904-01-01") # 1998-07-05
ウィンドウズ版: origin = "1899-12-30"
マック版 : origin = "1904-01-01"
というオプション指定だけを覚えてしまえば十分ですね。
これをXLConnectパッケージを使ってRで読み込むと、
> library(XLConnect)
> wb <- loadWorkbook("input.xlsx")
> ws <- readWorksheet(wb, sheet=1, header=F)
> ws
Col1
1 41842
のように、「2014年7月22日」と入力してあったセルの値が、「41842」となっています。
これはシリアル値といって、1900年1月1日を「1」、1月2日を「2」、・・・として日付を数値化したもので、エクセルは内部的にはこのように日付データを保持しています。というのを知っている人は多いと思います。
じゃあ、この値を日付としてRから使いたいときにどうやるのか? というのが今回のエントリの目的です。
Rには日付のシリアル表現のための関数があって、整数値と起源日付を指定することによって、シリアル値を日付(Date)型に変換することができます。
ということは、↓これでいけるはずだが・・・
> as.Date(41842, origin="1900-01-01") # 2014年7月22日のつもり
[1] "2014-07-24"
あれ? 2日ほどずれている。
このずれの原因の1つ目は以下のコードを実行することで、すぐに理解できます。
> as.Date(0, origin="1900-01-01")
[1] "1900-01-01"
> as.Date(1, origin="1900-01-01")
[1] "1900-01-02"
> as.Date(2, origin="1900-01-01")
[1] "1900-01-03"
エクセルでは1900年1月1日を「1」としているのに対し、Rのas.Date関数はoriginオプションで指定された日付(つまり1900年1月1日)を「0」としているわけですね。0オリジン vs 1オリジン のゴタゴタがこんなところにまで・・・
そして、ずれの原因の2つ目は、エクセルでシリアル値を実際に表示してみると判明します。
エクセル上でシリアル値を確認したい場合は、日付が入っているセルの書式設定で、表示形式タブを選択、ユーザー定義で「G/標準」を指定すればOKです。
絶対的な真理として、1900年はうるう年ではありません。つまり、1900年2月29日は存在せず、2月28日が59なら、3月1日を60にするべきなのです。なのに、そうなっていないせいで、1900年3月1日以降(あなたが扱うであろうほとんどの日付)に関しては、このずれが発生してしまうことになります。
このバグとも思えるような変な仕様は、(当時有力だった)「Lotus1-2-3」との互換性を維持するための意図的なものであったとのこと。
なんて、いきさつのことは知らなくても、やり方だけ覚えておけばOK。
【A】origin="1900-01-01"と指定するなら、シリアル値から2を引く または
というのでもいいけど、↓こっちの方がシンプルかも
【B】origin="1899-12-30"と指定する
【B】の方は、マニュアルにも書いてあって(?as.Dateと打てば見られます)、
## So for dates (post-1901) from Windows Excel
as.Date(35981, origin = "1899-12-30") # 1998-07-05
というサンプルが載っています。
マッキントッシュは起源日付が違うので、マック版エクセルでは、以下のようにすればいいらしいです。
## and Mac Excel
as.Date(34519, origin = "1904-01-01") # 1998-07-05
ウィンドウズ版: origin = "1899-12-30"
マック版 : origin = "1904-01-01"
というオプション指定だけを覚えてしまえば十分ですね。
2014年7月18日金曜日
Rの自己組織化マップ(SOM)で文字が重ならないようにする方法
Rでsom関数を使って自己組織化マップを作った後、各ユニットの座標に文字列をplotした場合、同じユニットに配置された文字列の重なりが気になることってありますよね。
下記は、放送大学の「データからの知識発見」という授業のテキストからの引用です。
入門者向けの講義なので、サンプルコードのシンプルさを優先したのか、文字の重なりについての対処はありません。用いられていた入力のデータ数(動物の数)が14、同じ場所に配置されてしまったデータが最大で2個なので、なんとか読み取れます。ライオンとトラが重なると「ライトオン」っぽくなるというのはトリビアですね。
ただ、データ数が多くなって、3つも4つも重なってくると、読み取れなくなって実害が出てくると思います。
「データマイニング入門」(わたせせいぞうの表紙のシリーズのやつ)では、そこに一工夫されていて、$visualで得られる座標をそのまま使うのではなく、乱数を加えることによって、文字が全く同じ位置で重なることを避けています↓
でも、そこはやはり乱数なので、うまく分離されているところもあれば、これじゃ読み取れないでしょってレベルのところもあります。
↓以下はデータマイニング入門に載っていたサンプルをそのまま動かしたものです。
カナリヤとスズメのところは書籍の図よりも見やすいようですが、キツネとタヌキのところは完全に重なってしまっていますね。乱数だからこんなこともあります。
で、私としては不確実なものに頼るよりは、座標が重なった(=同じユニットに配置された)なら、定量だけ座標をずらしてやればいいじゃんと思ったわけです。
上記の考えにしたがって作ったサンプルコードと実行例の図が↓以下です。
grpには、x座標とy座標を組にした因子でグループ分けされたリストが入ります。表示させると↓こんな感じです。
> grp
$`7:0`
[1] 5 13
$`3:1`
[1] 18
$`6:1`
[1] 15
・・・(略)・・・
(7, 0)に配置されたのは5番目と13番目のデータの2つ、(3, 1)に配置されたのは18番目のデータの1つ、・・・という具合です。
for文の中で、それぞれのグループに対し、1つ目のものはそのまま、2つ目のものには0.2を足す、3つ目のものには0.4を足す、という具合にy軸の表示位置を下にずらしていき、重なりを避けています。
0.2という数値は試行錯誤した結果で、cex=0.6と合わせて微調整してください。データ数が多いと、隣のユニットの領域にまではみ出すことになるので、その場合はpngやpdfに描き出すようにして、描画領域を大きく設定し、cexを小さくするようなやり方で対処できると思います。
このやり方で、理屈の上では、一切重なりのない自己組織化マップがRで描けるということになります。お試しあれ。
データからの知識発見 (放送大学教材)
データマイニング入門
下記は、放送大学の「データからの知識発見」という授業のテキストからの引用です。
重なり対処を全く行っていない表示例 |
入門者向けの講義なので、サンプルコードのシンプルさを優先したのか、文字の重なりについての対処はありません。用いられていた入力のデータ数(動物の数)が14、同じ場所に配置されてしまったデータが最大で2個なので、なんとか読み取れます。ライオンとトラが重なると「ライトオン」っぽくなるというのはトリビアですね。
ただ、データ数が多くなって、3つも4つも重なってくると、読み取れなくなって実害が出てくると思います。
「データマイニング入門」(わたせせいぞうの表紙のシリーズのやつ)では、そこに一工夫されていて、$visualで得られる座標をそのまま使うのではなく、乱数を加えることによって、文字が全く同じ位置で重なることを避けています↓
乱数による重なり対処を行った表示例 |
でも、そこはやはり乱数なので、うまく分離されているところもあれば、これじゃ読み取れないでしょってレベルのところもあります。
↓以下はデータマイニング入門に載っていたサンプルをそのまま動かしたものです。
乱数による重なり対処を行った表示例(私の環境で実行させたもの) |
library(som) 動物データ1 <- read.csv("animal1.csv", header=T) 標準化動物データ <- normalize(動物データ1[,2:14], byrow=F) 動物SOM <- som(標準化動物データ, xdim=10, ydim=10, topol="rect") 乱数 <- cbind(rnorm(nrow(動物データ1), 0, 0.15), rnorm(nrow(動物データ1), 0, 0.15)) 動物マップ <- 動物SOM$visual[, 1:2] + 乱数 + 0.5 plot(動物マップ, xlim=c(0,10), ylim=c(0,10)) text(動物マップ[, 1], 動物マップ[, 2], 動物データ1$動物名)
カナリヤとスズメのところは書籍の図よりも見やすいようですが、キツネとタヌキのところは完全に重なってしまっていますね。乱数だからこんなこともあります。
で、私としては不確実なものに頼るよりは、座標が重なった(=同じユニットに配置された)なら、定量だけ座標をずらしてやればいいじゃんと思ったわけです。
上記の考えにしたがって作ったサンプルコードと実行例の図が↓以下です。
定量だけずらすことによる重なり対処を行った表示例 |
# 同じユニットに配置されたものがあれば、定量だけずらすことにより # 重なりを避ける vis <- 動物SOM$visual # x座標とy座標の組を因子(ファクタ)として、グループに分ける grp <- split(1:nrow(vis), vis[,1:2], sep=":", drop=T ) for (i in 1:length(grp)){ # 同じ場所に配置されたものの添え字のベクトル idx <- grp[[i]] # 位置ずらし用の数 0.0, +0.2, +0.4, ... を足す vis[idx, ]$y <- vis[idx, ]$y + seq(0.0, 10, 0.2)[1:length(idx)] } plot(vis$x, vis$y, type="n", xlim=c(0,10), ylim=c(0,10)) text(vis$x, vis$y, 動物データ1$動物名, cex=0.6)
grpには、x座標とy座標を組にした因子でグループ分けされたリストが入ります。表示させると↓こんな感じです。
> grp
$`7:0`
[1] 5 13
$`3:1`
[1] 18
$`6:1`
[1] 15
・・・(略)・・・
(7, 0)に配置されたのは5番目と13番目のデータの2つ、(3, 1)に配置されたのは18番目のデータの1つ、・・・という具合です。
for文の中で、それぞれのグループに対し、1つ目のものはそのまま、2つ目のものには0.2を足す、3つ目のものには0.4を足す、という具合にy軸の表示位置を下にずらしていき、重なりを避けています。
0.2という数値は試行錯誤した結果で、cex=0.6と合わせて微調整してください。データ数が多いと、隣のユニットの領域にまではみ出すことになるので、その場合はpngやpdfに描き出すようにして、描画領域を大きく設定し、cexを小さくするようなやり方で対処できると思います。
このやり方で、理屈の上では、一切重なりのない自己組織化マップがRで描けるということになります。お試しあれ。
データからの知識発見 (放送大学教材)
データマイニング入門
2014年7月16日水曜日
RからYahoo!のジオコーディングを利用する方法
ジオコーディングというのは、
入力: 住所の文字列
出力: 緯度経度の数値
という処理のことです。YahooやGoogleがAPIを公開しています。
YahooのAPIの仕様は↓下記のページにあります。
YOLP(地図):Yahoo!ジオコーダAPI - Yahoo!デベロッパーネットワーク
このページの説明にしたがって、住所文字列をエンコードしてURLに埋め込んで、ブラウザのアドレス欄に入れれば結果が返ってきます。
でも、こういうのはプログラムから使いたいですよね。で、私が多少なりとも使いなれているRでのやり方を書きとめておきます。
普通に使うときは、一連の処理を関数の形にして、それを呼び出すような書き方をする方がすっきりするのですが、今回は使い方の紹介がメインなので、あえてベタ書きしています。
あと、アプリケーションIDは各人で取得する必要があります。YahooのIDを持っていれば(持っていない場合は作れば)、無料で取得できます。
RESTについての前提知識があって、前述のリンク先の仕様説明が理解できて、RCurlの使い方を知っていれば、普通にできることなんでしょうけど、私はそうではなかったので多少苦労しました。
私がつまづいたのは、文字コード指定とプロキシ指定です。
私はWindowsを使っていて、YahooのAPIはデフォルトではUTF-8を入力することになっているので、iconvしてやる必要がありました。もしくは、リクエストの「ei」パラメタで文字コード(WindowsならSJIS)を指定する方法でもいいです。
プロキシを使っていなければ、そのまま getURL(url.u)のようにリクエストを投げればOKなのですが、私の環境ではプロキシを使っているので指定する必要がありました。サンプルでlocalhostになっているのは、私の環境ではプロキシ中継ツールを使っているからです。この箇所にご自分の環境に合わせてproxyのURLを指定してやればOKです。
プロキシサーバが認証を必要とする場合は、
XXXXX: ID
YYYYY: パスワード
ZZZZZ: プロキシURL
のようにプロキシURLの前にIDとパスワードをつければOKです。
fromJSON()関数を使うと、JSON形式のオブジェクトを、リスト形式に変換してくれるようです。正直、JSON形式のことはよく理解していないのですが、リストになってしまえば、私にもアクセス方法が分かりました。
x(JSONからリストに変換されたもの)を表示させてみると、リストがどういう内容を持っているか分かるので、あとは
x$Feature[[1]]$Geometry$Coordinates
のように、各要素にアクセスすればOKですね。
サンプルに書いている
addr <- "東京都墨田区押上1-1-13"
は、東京スカイツリーの住所なのですが、これで試してみると、
> (lat <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[2])
[1] "35.71000629"
> (lng <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[1])
[1] "139.81070961"
> (location_type <- x$Feature[[1]]$Property$AddressType)
[1] "地番・戸番"
> (formatted_address <- x$Feature[[1]]$Property$Address)
[1] "東京都墨田区押上1丁目1-13"
のような感じで結果が返ってきます。
AddressTypeを見ると情報の粒度が分かります。
「東京都墨田区押上1-1-13」 → 「地番・戸番」
「東京都墨田区押上」 → 「町・大字」
「東京都墨田区」 → 「特別区」
のように値が設定されます。
また、Addressを見ると、どういう住所だと認識して結果を返してきたのかを確認することができます。
入力: 住所の文字列
出力: 緯度経度の数値
という処理のことです。YahooやGoogleがAPIを公開しています。
YahooのAPIの仕様は↓下記のページにあります。
YOLP(地図):Yahoo!ジオコーダAPI - Yahoo!デベロッパーネットワーク
このページの説明にしたがって、住所文字列をエンコードしてURLに埋め込んで、ブラウザのアドレス欄に入れれば結果が返ってきます。
でも、こういうのはプログラムから使いたいですよね。で、私が多少なりとも使いなれているRでのやり方を書きとめておきます。
普通に使うときは、一連の処理を関数の形にして、それを呼び出すような書き方をする方がすっきりするのですが、今回は使い方の紹介がメインなので、あえてベタ書きしています。
あと、アプリケーションIDは各人で取得する必要があります。YahooのIDを持っていれば(持っていない場合は作れば)、無料で取得できます。
library(RCurl) library(RJSONIO) # appidのところは自分で取得したものを入れてください root <- "http://geo.search.olp.yahooapis.jp/OpenLocalPlatform/V1/geoCoder?appid=xxxxxxxxxx&output=json" addr <- "東京都墨田区押上1-1-13" url <- paste(root, "&query=", addr, sep="") url.u <- URLencode(iconv(url, "", "UTF-8")) # UTF-8に変換 # プロキシへの対処(使っていない場合はgetURL(url.u)とそのまま呼んでOK) curl <- getCurlHandle() curlSetOpt(.opts=list(proxy="ZZZZZ:8080"), curl=curl) #プロキシURLは環境に合わせて doc <- getURL(url.u, curl=curl) # curlの指定を追加 x <- fromJSON(doc, simplify = FALSE) lat <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[2] lng <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[1] location_type <- x$Feature[[1]]$Property$AddressType formatted_address <- x$Feature[[1]]$Property$Address
RESTについての前提知識があって、前述のリンク先の仕様説明が理解できて、RCurlの使い方を知っていれば、普通にできることなんでしょうけど、私はそうではなかったので多少苦労しました。
私がつまづいたのは、文字コード指定とプロキシ指定です。
私はWindowsを使っていて、YahooのAPIはデフォルトではUTF-8を入力することになっているので、iconvしてやる必要がありました。もしくは、リクエストの「ei」パラメタで文字コード(WindowsならSJIS)を指定する方法でもいいです。
プロキシを使っていなければ、そのまま getURL(url.u)のようにリクエストを投げればOKなのですが、私の環境ではプロキシを使っているので指定する必要がありました。サンプルでlocalhostになっているのは、私の環境ではプロキシ中継ツールを使っているからです。この箇所にご自分の環境に合わせてproxyのURLを指定してやればOKです。
プロキシサーバが認証を必要とする場合は、
XXXXX: ID
YYYYY: パスワード
ZZZZZ: プロキシURL
のようにプロキシURLの前にIDとパスワードをつければOKです。
fromJSON()関数を使うと、JSON形式のオブジェクトを、リスト形式に変換してくれるようです。正直、JSON形式のことはよく理解していないのですが、リストになってしまえば、私にもアクセス方法が分かりました。
x(JSONからリストに変換されたもの)を表示させてみると、リストがどういう内容を持っているか分かるので、あとは
x$Feature[[1]]$Geometry$Coordinates
のように、各要素にアクセスすればOKですね。
サンプルに書いている
addr <- "東京都墨田区押上1-1-13"
は、東京スカイツリーの住所なのですが、これで試してみると、
> (lat <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[2])
[1] "35.71000629"
> (lng <- unlist(strsplit(x$Feature[[1]]$Geometry$Coordinates, ","))[1])
[1] "139.81070961"
> (location_type <- x$Feature[[1]]$Property$AddressType)
[1] "地番・戸番"
> (formatted_address <- x$Feature[[1]]$Property$Address)
[1] "東京都墨田区押上1丁目1-13"
のような感じで結果が返ってきます。
AddressTypeを見ると情報の粒度が分かります。
「東京都墨田区押上1-1-13」 → 「地番・戸番」
「東京都墨田区押上」 → 「町・大字」
「東京都墨田区」 → 「特別区」
のように値が設定されます。
また、Addressを見ると、どういう住所だと認識して結果を返してきたのかを確認することができます。
2014年7月1日火曜日
e-Stat(統計局)で公開されているShapeファイルを、Rで表示する
まずはShapeファイルを入手するところから。
総務省統計局がe-Statというサイトを運営していて、そこに全国各地域のShapeファイルが公開されています。
地図で見る統計(統計GIS)
(1) 「平成22年国勢調査(小地域) 2010/10/01」を選択
(2) 「男女別人口総数及び世帯総数」にチェック
(3) 「統計表各種データダウンロードへ」ボタンをクリック
(1) 「都道府県」を選択
(2) 「市区町村」を選択
(3) 「検索」ボタンをクリック
(4) 境界データの欄に表示されている
「世界測地系緯度経度・Shape形式」のデータをクリック
ダウンロードしたzipファイルを展開すると、
xxx.shp
xxx.dbf
xxx.shx
xxx.prj
という、拡張子だけが違うファイルが4つできるんですが、これがshapeデータの1セットになります。
それぞれのファイルの中身は、シェープファイルとは? | 用語集とGISの使い方 | 株式会社パスコによると、
shp: 図形の座標が保存
dbf: 属性の情報が保存
shx: shpの図形とdbfの属性の対応関係が保存
だそうです。あと、「.prj」は測地系の情報が入っています。(テキストとして開くとそれっぽい文言が入っている)
では単純にshapeファイルをRで表示させてみましょう。
あら簡単。実質4行で表示できました。
じゃあ、いろいろやってみましょう。
coordinates()という関数は、ポリゴンの重心の位置を返してくれます。
> coordinates(x)
[,1] [,2]
0 139.4639 35.62169
1 139.5178 35.62177
2 139.5138 35.62332
3 139.5057 35.62163
4 139.4617 35.61308
... ...
そのままプロットすると↓こんな感じですね。
plot(x)
points(coordinates(x))
このshapeを読み込んだオブジェクト(SpatialPolygonsDataFrame)には、境界のポリゴンだけでなく、属性情報も入っています。
どうやら、この「MOJI」というカラムにそのポリゴンに対応する地名が入っているようです。これを地図上にプロットしてみましょう。プロット位置はさっきのやり方で求めた重心の位置にしてみます。あと、MOJIは因子型で入っているので、as.character()で文字列に変換してやる必要がありました。
plot(x)
text(coordinates(x), as.character(x@data$MOJI))
うわ、ぐちゃぐちゃ。
字を小さくしてみよう↓
plot(x)
text(coordinates(x), as.character(x@data$MOJI), cex=0.5)
今度は、字が小さくて見づらいですね。
画面に出すだけのときは、描画用のウィンドウをググッと広げると、ポリゴンは大きく、文字はそのままで表示される、見やすくなります。
直接、画像ファイルとして出力し、その際に解像度を高めに指定するというやり方もあります。
解像度やら文字の大きさやらは、お好みで調整してみてください。
他にもいろいろと遊べるのですが、長くなったので、今回はここまで。■
[2016年3月30日追記]
広いエリアの行政界シェープファイルのようなものを画像ファイルに出力するのは、けっこう無理があります。そのようなときはPDFファイルに出力すると、いい感じになります↓
Rでshapeを出力するときはPDFファイルにすると便利 - Rプログラミングの小ネタ
■
総務省統計局がe-Statというサイトを運営していて、そこに全国各地域のShapeファイルが公開されています。
地図で見る統計(統計GIS)
(1) 「平成22年国勢調査(小地域) 2010/10/01」を選択
(2) 「男女別人口総数及び世帯総数」にチェック
(3) 「統計表各種データダウンロードへ」ボタンをクリック
(1) 「都道府県」を選択
(2) 「市区町村」を選択
(3) 「検索」ボタンをクリック
(4) 境界データの欄に表示されている
「世界測地系緯度経度・Shape形式」のデータをクリック
ダウンロードしたzipファイルを展開すると、
xxx.shp
xxx.dbf
xxx.shx
xxx.prj
という、拡張子だけが違うファイルが4つできるんですが、これがshapeデータの1セットになります。
それぞれのファイルの中身は、シェープファイルとは? | 用語集とGISの使い方 | 株式会社パスコによると、
shp: 図形の座標が保存
dbf: 属性の情報が保存
shx: shpの図形とdbfの属性の対応関係が保存
だそうです。あと、「.prj」は測地系の情報が入っています。(テキストとして開くとそれっぽい文言が入っている)
では単純にshapeファイルをRで表示させてみましょう。
境界情報だけをプロット |
あら簡単。実質4行で表示できました。
じゃあ、いろいろやってみましょう。
coordinates()という関数は、ポリゴンの重心の位置を返してくれます。
> coordinates(x)
[,1] [,2]
0 139.4639 35.62169
1 139.5178 35.62177
2 139.5138 35.62332
3 139.5057 35.62163
4 139.4617 35.61308
... ...
そのままプロットすると↓こんな感じですね。
重心の位置をプロット |
points(coordinates(x))
このshapeを読み込んだオブジェクト(SpatialPolygonsDataFrame)には、境界のポリゴンだけでなく、属性情報も入っています。
> head(x@data) AREA PERIMETER ... KEN_NAME ... MOJI ... 0 206286.40 2421.068 ... 神奈川県 ... はるひ野5丁目 ... 1 67760.92 1765.906 ... 神奈川県 ... 細山6丁目 ... 2 119563.80 1457.849 ... 神奈川県 ... 細山7丁目 ... 3 225729.80 2633.794 ... 神奈川県 ... 細山 ... 4 2017504.00 11135.110 ... 神奈川県 ... 黒川 ... 5 95170.74 1569.407 ... 神奈川県 ... 細山8丁目 ... ... ... ... ... ... ... ...
どうやら、この「MOJI」というカラムにそのポリゴンに対応する地名が入っているようです。これを地図上にプロットしてみましょう。プロット位置はさっきのやり方で求めた重心の位置にしてみます。あと、MOJIは因子型で入っているので、as.character()で文字列に変換してやる必要がありました。
地名を表示させたが、文字サイズが大きすぎた |
text(coordinates(x), as.character(x@data$MOJI))
うわ、ぐちゃぐちゃ。
字を小さくしてみよう↓
文字サイズが小さいと見にくい・・・ |
text(coordinates(x), as.character(x@data$MOJI), cex=0.5)
今度は、字が小さくて見づらいですね。
画面に出すだけのときは、描画用のウィンドウをググッと広げると、ポリゴンは大きく、文字はそのままで表示される、見やすくなります。
直接、画像ファイルとして出力し、その際に解像度を高めに指定するというやり方もあります。
高解像度のpngに出力(拡大してみればいい感じです) |
解像度やら文字の大きさやらは、お好みで調整してみてください。
他にもいろいろと遊べるのですが、長くなったので、今回はここまで。■
[2016年3月30日追記]
広いエリアの行政界シェープファイルのようなものを画像ファイルに出力するのは、けっこう無理があります。そのようなときはPDFファイルに出力すると、いい感じになります↓
Rでshapeを出力するときはPDFファイルにすると便利 - Rプログラミングの小ネタ
■
登録:
投稿 (Atom)