ホテル街を見つける（２）RとGoogleAPIでジオコーディング
前回に続いて、ホテル街を見つける技術編第二弾。前回は住所一覧のベクトルをcsvファイルにしまうところまでやりました。
今回は住所一覧を元に緯度経度をゲットする ジオコーディング のやり方について。方法としてはRからGoogle APIを利用します。
流れ
今回の最終目標は以下の？を埋めたcsvファイルを吐き出すところまで。
| 住所 | 緯度 | 経度 |
|---|---|---|
| 北海道なんたら１ | ？ | ？ |
| 北海道なんたら２ | ？ | ？ |
| 北海道なんたら３ | ？ | ？ |
| … | ？ | ？ |
参照
ジオコーディングをやってみる
ジオコーディングとは、住所を緯度・経度の座標値に変換すること。住所のままだとプロットの都合上塩梅が悪いので、緯度経度のセットに変換したい。例えば、東京駅の住所がわかっていたとして、わからない緯度と経度をひっぱってきたい。
| 住所 | 緯度 | 経度 |
|---|---|---|
| 東京都千代田区丸の内一丁目 | ？ | ？ |
何はなくともジオコーディングしてみよう（ここは先の参照サイト)のコードをほぼそのままお借りしてます）。
# ライブラリ読み込み library(rjson) library(RCurl) # 住所 hoge <- "東京都千代田区丸の内一丁目" # 欲しい住所を指定したURLを作成 url <- paste0("http://maps.googleapis.com/maps/api/geocode/json?address=", hoge, "&sensor=false") # URLを適切にエンコード url_utf8 <- URLencode(iconv(url, "", "UTF-8")) # URLからJSONデータを取得し、リスト形式に変換 getJSON <- fromJSON(getURL(url_utf8)) # リストの中から緯度と経度を引っ張りだす getJSON$results[[1]]$geometry$location
結果として以下の様に緯度と経度を獲得することができる。
| 住所 | 緯度 | 経度 |
|---|---|---|
| 東京都千代田区丸の内一丁目 | 35.68332 | 139.7674 |
paste0()は引数同士の間に何も挟まず、文字列を結合する（参照）。incov()は第1引数に与えた文字列を、第2引数のコードから第3引数のコードへ変換する（参照）。ここではUTF-8へと変換している。さらにURLencode()でURLエンコードへ変更する。よくわからないが、URLに入っていてはまずい文字を変換するっぽい（参照)。
ホテルの緯度・経度を獲得する
前回の記事で以下のように住所一覧を得て address.csv にしまった。今回はこの緯度と経度の？をすべて埋めたcsvファイルを作ることが目標である。
| 住所 | 緯度 | 経度 |
|---|---|---|
| 北海道なんたら１ | ？ | ？ |
| 北海道なんたら２ | ？ | ？ |
| 北海道なんたら３ | ？ | ？ |
| … | ？ | ？ |
やることは至って単純で、先ほどの「東京都千代田区丸の内一丁目」(= hoge)のところへ、順番に住所を投げ込んでいけば良い。乗っかかりっきりですが、Rで(逆)ジオコーディングさんのコードのマイナーチェンジで挑みます（感謝）。
# データの読み込み text <- read.csv("address.csv",header=T,stringsAsFactors = FALSE) text <- text$x # ライブラリ読み込み library(RCurl) library(rjson) library(RgoogleMaps) # 関数作成 GEO_new <- function(address, para = list(sensor = "false", language = "ja", region = "JP"), proxy = NULL) { curl_set <- getCurlHandle() if(!is.null(proxy)) curlSetOpt(.opts = list(proxy = proxy), curl = curl_set) parameters <- paste0("&", names(para), "=", unname(unlist(para)), collapse = "") GeoRequest <- iconv(paste0("http://maps.googleapis.com/maps/api/geocode/json?address=", address, parameters), "", "UTF-8") geodata <- vector("list", length(address)) for(i in seq_along(address)) { getJSON <- fromJSON(getURL(URLencode(GeoRequest[i]), curl = curl_set)) if (length(getJSON$results)!=0){ geodata[[i]] <- as.data.frame(getJSON$results[[1]]$geometry$location) Sys.sleep(0.5) print(paste(i,"=================="))} else { geodata[[i]] <- NA Sys.sleep(0.5) print(paste(i,"失敗した…"))}} cbind(address, do.call("rbind", geodata)) } # ジオコーディング：3回に分けて3つのセットにしまう(ちなみにうまくいかないと思うので下記参照) location0001_2000 <- GEO_new(text[1:2000]) location2001_4000 <- GEO_new(text[2001:4000]) location4001_6088 <- GEO_new(text[4001:6088]) # 縦につないで一つのデータセットに location <- rbind(location0001_2000,location2001_4000,location4001_6088) # 書き出し write.csv(location, "location.csv", quote=FALSE, row.names=FALSE)
大きな流れは、データを読み込んで text ベクトルに住所を入れ込む。関数 GEO_new()を定義する。2000ずつジオコーディングした住所をしまいこんで、縦にくっつける。csvで書き出す。
変更点は以下の3点。第1にエラーの処理した。第2に現在の進行状況の確認と住所の獲得に失敗したかどうかの表示した。第3にAPIのリクエスト間に時間を設けた。重要な注意点もあるので順に説明する。
① エラーを吐いても止めない
当初のコードだと途中で エラーを吐いた場合止まる 。該当住所がないのか何なのか、多くはないが住所の獲得に失敗しエラーを吐いてしまった。そこでfor文の中にif文をはさみ、エラーを吐いた場合は該当箇所にNAを入れこむ仕様にした。これで止まらない。中身としては先に述べたように、getJSONにリストが入る仕様になっているのだが、失敗するとgetJSON$resultsにlist()が入ってしまって[[1]]にアクセスできずエラーが起きる。これを防いでいる。
② 何ケースやったか/どこで失敗したか
print()を挟み込んで 何ケース進んだかを表示 するようにした。また、 失敗した時は「失敗した…」と表示 するようにした。そのため、どのケースがNAになっているのかとりあえずわかる。
③ APIリクエスト制限
リクエストを過剰に送るとgoogle側からstopが入りやはりエラーする 。具体的には、
IPあたり24時間に2500リクエスト、秒間5リクエスト
参照 意外に多いGoogle Maps関連機能制限をまとめた、Google Geocoding APIの制約に関するよくある誤解
そこで、まず2500件制限に対策すべく上記のように 2000件ずつのアクセス にしている。上記のコードを一度に回すとおそらくGEO_new(text[2001:4000])が500辺り以降ですべてエラーになるためNAが入ってしまう。3日間に分けて獲得するか、他のIPアドレスを用いて獲得しよう。yahooのサービスは制限がないためそちらを使っても良い、みたいな話もある。
次に秒間リクエストだが、Sys.sleep(0.5)で アクセスごとに0.5秒止める ことにした。理論上は0.2でもいいはずなのだが、時々失敗するので試行錯誤の末0.5とした。この辺りはちょっといじったほうがいいかもしれない。
まとめ
以上で以って以下の様な location.csv ファイルがディレクトリにできたはず。やったね！
| address | lat | lng |
|---|---|---|
| 北海道なんたら１ | xx.xxxxx | xx.xxxxx |
| 北海道なんたら２ | xx.xxxxx | xx.xxxxx |
| 北海道なんたら３ | NA | NA |
| … | xx.xxxxx | xx.xxxxx |
次回はいよいよこのデータを使ってQGISで地図上へのプロット、および「GISデータを用いた分析」っぽい分析を少し行います。QGISのダウンロードはあちこちに良い解説があるのでそちらを参照してください。ではでは。
人口・駅データをプロット
目的はQGISを使って地図上に人口および駅をプロットすること。いろいろな分析の基礎になりますからやっておいて損はなさそうですね。QGISはインストール済みを想定しています。
QGISはレイヤを重ねていくことで表示するので順番に3つを重ねていきます。
- 日本地図
- 人口
- 駅・路線図
日本地図
国土数値情報の行政区域からshapeファイルをゲット。
レイヤ > レイヤの追加 > ベクタレイヤの追加 > ブラウズ
で該当する.shpを開けばこんな感じになります。ひとまず日本地図のプロットに成功！
人口：国勢調査データ
つづいて人口について。ここ参照です。国勢調査のデータを用います。流れは以下。
国勢調査データのダウンロード
e-Statから、国勢調査500mメッシュを取得。自分の欲しい該当コードに当たるデータをダウンロード。今回は東京を落としてみます。番号は5339。
shape形式をダウンロード。すると、「H002005112009DDSWH05339」にはshapeデータとidが、「tblT000609H53390」の方には.txt形式で国勢調査データが入っている。このとき.txtファイルの拡張子を.csvに書き換えておく。
読み込みとひも付け
読み込み
shapeファイル、csvをそれぞれベクターデータで読み込む。これで各地域を500m単位で四角く区切ったメッシュデータ(shapeファイル)と、各四角区域に何人がいるのか(csvファイル)がそれぞれ読み込まれた。
ひも付け
つづいて、各四角(shape)と人口(csv)を紐付ける必要がある。
メッシュの方(MESH5339)を右クリック。
プロパティ > 結合 > +
結合するレイヤに「tblT000609H53390」、結合するフィールド・ターゲットフィールドに「KEY_CODE」を選択してひも付け完了。MESH5339を右クリック、属性テーブルを開いてやると紐付いているはず。
人口の数値化
とここで注意なのが、人口データがstring形式になっている…。stringは文字なので塗り分けのときに数値として認識してくれない！ちょっとうまいやり方がわからないので、新しい数値の変数を作成し、そこに結合した人口を代入してやることにする。
MESH5339を右クリック、属性テーブルを開く。左上の鉛筆のアイコン「編集モード切替」を押す。続いて「新規カラムを作る」を押す。名前に「people」とかつけて、 タイプに「integer」を選ぶ （このために変数を作ってるのでココ重要）。
これで変数ができるはず。上にあるボックスで左から「people」を選び、右側に「tblT000609H53390_T000609001」(要するにひも付けた最初の列)を選んで右にある「全部更新」のボタン。暫し待つ（反映までちょっと時間かかる）。最後に左上の鉛筆のアイコン「編集モード切替」を押して編集モードを修了する。
人口で塗り分け
いよいよMESH5339を右クリック、オプションでスタイル。「段階にわけられた」を選んでカラム「people」を選ぶ。色階調を適当に選んで、分類数とか適当にいじり分類ボタン。OKをおすとこうなります（やったね！）。
駅・路線図：国土数値情報
駅と路線図をプロット
最後に駅と路線図を重ねる。国土数値情報で、鉄道を押して流れに任せてダウンロード。駅と路線を先ほどと同じ要領でベクターで読み込み。
駅のデータと路線のラインがそれぞれ別のデータで得られるのでレイヤーでチェックを入れて重ねる。
駅名を表示
Station 右クリック > オプション > ラベル
ラベルにチェックを入れて「N02_005」を選んでOKを押すと駅名も表示されます。結果が以下！
あとは人口の人数区切りを変えたり、地域を変えたり、路線等々の色を変えたりは個人の好みで。QGISはプロットがきれいやね。
過去50ツイート時間をリプで返す
リプライで「おい」と飛んできたら、飛ばしてきたアカウントの過去50ツイートのツイート時間を返すbotの作成。下記の機能をtwitterに実装したものです。
こんな感じを目指します。
#!/usr/local/bin/ruby -Ku require 'kconv' require 'twitter' client = Twitter::REST::Client.new do |config| config.consumer_key = "なんとか" config.consumer_secret = "かんとか" config.access_token = "うんとか" config.access_token_secret = "すんとか" end recentid = ツイート取得id初期値 #client.mentions_timeline[0].id とかで得られます loop{ ripu = client.mentions_timeline(:since_id => recentid) #リプライを取得 if ripu.size != 0 for i in 0..(ripu.size-1) if ripu[i].text.index(/おい/).class==Fixnum #「おい」を含むか判定 user_id = ripu[i].user.id #カウント配列に0をセット tweet_counts = [] (0..23).each { |j| tweet_counts[j] = 0 } #順番に該当時をカウント（本当は50より大きくしたいけどこれ以上にすると140字を超える） client.user_timeline(user_id, { count: 50 } ).each do |tweet| tweet_counts[tweet.created_at.hour] += 1 end #リプ内容を作る（#のヒストグラム） hoge = "" (0..23).each do |j| hoge = hoge + "%02d"% j+"#" * tweet_counts[j] + "\n" end post = "@#{ripu[i].user.screen_name} \n#{hoge}" twi_id=ripu[i].id client.update(post, {:in_reply_to_status_id => twi_id }) end end # recentidの更新 recentid = ripu[0].id end # 60秒ごとに稼働させる sleep(60) }
bot を heroku で動かす
Rubyでbot を作ったはいいが crontab は PC が起きてないと動かない。そこでサーバー上で bot を動かすことを試みる。
上記を参照して同じフォルダに5つのファイルを準備し、ターミナルにて以下のプログラムを走らせる。ブラウザでherokuへログインしてbotのスイッチをオンすれば完了。
- Procfile
- Gemfile
- app.rb
- config.ru
- botfile.rb(自分のプログラム)
git init heroku create botfile bundle install git add . git commit -m'commit' git push heroku master heroku open
heroku上のアプリを削除
失敗したら以下のコマンドで消してやり直そう。
heroku apps:destroy --app アプリ名 --confirm アプリ名
注意すべきポイント
何度も heroku へのデプロイを試みている際、フォルダ内の「Gemfile.lock」ファイルは毎度消してから bundle を行うこと。
Stataでシークエンス分析
Sequence Analysis のやり方を日本語で紹介したものがみあたらないので、Stataのパッケージマニュアルを参考に、紹介してみる。
Sequence Analysisとは？
Seaquence Analysisとは、並び順のパターンを発見する分析である。つまり、
11111
11311
11122
77777
のように4つの並び順があるとすると、上二つはかなり似ていて、四つ目はかなり違う。三番目はそこそこ似ている。みたいなのを発見する分析。手順としては以下の二段階を踏む。
- 配列同士の距離を計算する
- 距離からクラスター分析でいくつかのパターンに分類する
今回はStataを用いてSAを行う。
データの構築とパッケージの導入
データの構築
まずはパッケージの導入からデータの加工まで。とりあえず、定型データ（全てのサンプルが同じ長さ）形式にする^1。つまり、
11111
11311
221
のように３つ目だけ短い、みたいなのを避けて全部同じ長さのデータを作る。ひとまず、３つ目の空白に9とか代入して以下の形へ。
11111
11311
22199
パッケージの導入
なんはなくともパッケージの導入。
*** パッケージ導入 ssc install sq
データはロング型で与える。つまり、以下みたいな形でズラーッと長く^2。
| ID | order | みたい変数 |
|---|---|---|
| 1 | 1 | 1 |
| 1 | 2 | 1 |
| 1 | 3 | 1 |
| 1 | 4 | 1 |
| 1 | 5 | 1 |
| 2 | 1 | 1 |
| 2 | 2 | 1 |
| 2 | 3 | 3 |
| … | … | … |
つづいて、配列が入ってる変数と、各配列のid、配列の順序orderを指定して、sqsetの形を指定する（xtsetやstsetみたいにidや順序を指定したデータセット)。
*** データ形式指定
sqset みたい変数 id order, ltrim
分析
基礎分析
まずは基礎的な分析を行う。いわゆる記述統計的な分析。
| 関数 | 意味 | 例 |
|---|---|---|
| sqdes | 全要素をばらしたときの頻度を表しているはず^3 | |
| sqtab | 出現配列を出力。オプションは、rankは上位のみを出力。soは例えば11112と12222を同じものとしてカウント。seは12121と11112を同じものとしてカウント。 | |
| sqlength | 配列の長さ | 12121であれば5である。elementオプションは |
| sqelemcount | 配列内で異なる要素の数 | 11111であれば1、12121であれば2。 |
| sqepicount | 異なるエピソードの数。つまり数値の変化が何回起こったか-1 | 11111であれば1、12121であれば5。 |
*** 要素の頻度（長さ1の要素の頻度） sqdes *** 出現配列一覧 sqtab * top10のみ表示 sqtab,ranks(1/10) * 長さを考慮せずどういう順列か（ same order） sqtab,ranks(1/10) so * パターン構成要素が同じ配列の頻度（順不同でいかなる要素によって構成されるか same elements） sqtab,ranks(1/10) se * シークエンスの長さ egen length = sqlength() egen length1 = sqlength(),element(2) * 異なる要素の数 egen elemnum = sqelemcount() * エピソード単位のカウント egen epinum = sqepicount() * リスト表示で変数一覧 sqstatlist * 記述統計 sqstatsum * 1変数 sqstattab1 * クロス表 sqstattab2
図示
シークエンス分析の強みは何と言っても図示。
*** Sequence index plots
sqindexplot
シークエンス分析
以降では、いよいよ配列間の距離を計算しクラスタリングする。ここで重要なのは以下の点。
- 配列間の距離をいかに決定するか
- いかにクラスタリングするか
距離を求める
まず、配列間の距離は Levenshtein距離 を用いる。Levenshtein距離は、 indel-cost と substitution cost の2つの合計である。indel-cost(insert delete)とは、配列のある要素を挿入するないし、消すことである。例えば1231→131(消去)や、1231→12231(挿入)を指す。一方 substitution cost はある要素を他の要素に変更することである。例えば1231→1241（3を4に変更）。これらの操作それぞれにコストを決め（例えばindel-costを1、substitution costを2）、「コスト × 行われた手順」が Levenshtein距離となる。
ただし、これらの手順が一通りとは限らない。そこで、最小距離を求めるアルゴリズムが Needleman-Wunsch アルゴリズムである。
まずsqomで距離を求める。全ての距離を計算するオプションがfull。結構時間がかかる。名前を指定しなければSQDistというマトリックスに距離行列がしまわれる。
よくわからないが、データの形式上 sqclusterdat をして特殊な画面？に移動する。そこでクラスター分析を行い、sqclusterdat,returnで元の画面に戻って来る。クラスター分析の際、何クラスで割り当てるかはgroupの中で指定する。indexplotをクラスターごとに出力したところで今日はここまで。
*** 距離を計算（OM） sqom, full *** クラスター分析 sqclusterdat * 階層的クラスター分析（今回はウォード法で行う） * クラスターを計算 clustermat wardslinkage SQdist, name(wards) add * デンドログラムの出力 cluster tree wards, cutnumber(20) * カテゴリー値の割り当て cluster generate cl3=group(3), name(wards) cluster generate cl4=group(4), name(wards) cluster generate cl5=group(5), name(wards) cluster generate cl6=group(6), name(wards) sqclusterdat, return sqindexplot, by(cl3)
