TL;DR

BigQueryの fh-bigquery:geocode.geolite_city_bq_b2b Tableを使い IPアドレスから緯度経度を取得できます
国土地理院の「地球地図日本」から地図データを取得することができます
GeoPandasを使うことでお手軽に地理情報を可視化できます
BigQuery GISを使うと制限はありますが、もっと簡単に可視化することができます
↓BigQuery GISを使った可視化です
コードは全てGithubに置いてあります。

アクセスログ解析でのIPアドレス

アクセスログには基本的にIPアドレスが格納されていますが、アクセスログ分析する際に使っていないことも多いのではないでしょうか？
実はIPアドレスとBigQueryを上手く活用することで簡単に地理情報解析が行えるのです。

データ

実際にはアクセスログなどを使うことになると思いますが、テスト用に使える手頃なアクセスログを持ち合わせてないので、BigQuery内にあるsampledataのwikipediaテーブルを使います。wikipediaテーブルには編集者のIPなどが記録されているcontributor_ipというカラムがあります。

https://bigquery.cloud.google.com/table/publicdata:samples.wikipedia

contributor_ip: Typically, either ip or (id and username) will be set. IP information is unavailable for edits from registered accounts. A (very) small fraction of edits have neither ip or (id and username). They show up on Wikipedia as "(Username or IP removed)".

今回は、これを使って「どの地域からのWikipediaの編集が多いか」を可視化してみたいと思います。

Geoip geolocation with Google BigQuery

まずIPアドレスと位置情報の紐づけを行います。紐づけを行えるサービスというものが幾つかあって、最初はipstackなどを使っていました。ipstackは無料で月1000アクセス紐づけを行えるサービスで、手早く試せたので良かったのですが…

ipstack.com

途中で、BigQueryにIPアドレスと位置情報の紐づけ用のTableがあることを知り、こちらを使うことにしました。無料ですし紐づけも高速なので、もしデータがBigQuery上にあるなら使わない手はないと思います。

cloudplatform.googleblog.com

下記のようなSQLで fh-bigquery.geocode.geolite_city_bq_b2b にJOINし、緯度経度と都市名、それぞれの場所ごとに何件データがあったかを取得します。IPアドレスと位置情報は都市ごとに紐づけが行われます。

WITH safe AS (
  SELECT
    NET.SAFE_IP_FROM_STRING(contributor_ip) AS ip,  
    LENGTH(NET.SAFE_IP_FROM_STRING(contributor_ip)) AS len 
  FROM
      `publicdata.samples.wikipedia` 
  WHERE contributor_ip IS NOT NULL
),
wiki_ips AS (
  SELECT 
    NET.IPV4_TO_INT64(ip) AS clientIpNum,
    TRUNC(NET.IPV4_TO_INT64(ip)/(256*256)) AS classB_
  FROM safe 
  WHERE ip IS NOT NULL 
    AND len = 4
)
SELECT MAX(b.longitude) AS longitude, MAX(b.latitude) AS latitude, b.city 
FROM wiki_ips  
LEFT JOIN `fh-bigquery.geocode.geolite_city_bq_b2b`AS b ON wiki_ips.classB_ = b.classB
WHERE wiki_ips.clientIpNum BETWEEN b.startIpNum AND b.endIpNum
GROUP BY b.city 
HAVING MAX(b.country) = 'JP'
;

都市ごとにカウントすると下記のようになります。当然といえば当然ですが、東京が圧倒的ですね。5位に高知が来ているのは少し気になります。とてもアクティブな編集者の方がいらっしゃるのでしょうか？

f:id:mergyi:20180916021248p:plain

ちなみに世界版だとこんな感じでした。ロンドンが1位なのはちょっと予想外でした。

f:id:mergyi:20180916030709p:plain

国土地理院地図データ「地球地図日本」

それではこのデータを地理情報として可視化してみたいと思います。Pythonで地図をPlotするためには地図データが必要です。地図データは国土地理院が「地球地図日本」として公開してくれているのでこちらを使います。

地球地図日本｜国土地理院

下記のダウンロードスクリプトを使うとローカルに日本地図だけでなく路線地図などもダウンロードしてこれます。

import os
import zipfile
from pathlib import Path
import requests
from tqdm import tqdm
import pandas as pd 

data_dir = Path(os.getcwd())/'data'

# ZipFileがあるURLエンドポイント
gm_jpn_url = 'http://www1.gsi.go.jp/geowww/globalmap-gsi/download/data/gm-japan/gm-jpn-all_u_2_2.zip'

r = requests.get(gm_jpn_url, stream=True)
filename = data_dir/'gm-jpn.zip'

# 容量が大きいのでchunk_sizeごとにダウンロード
with open(filename, 'wb') as f:
    for chunk in tqdm(r.iter_content(chunk_size=1024)):
        if chunk:
            f.write(chunk)
            f.flush()

# Zipファイルを解答
zfile = zipfile.ZipFile(filename)
zfile.extractall(data_dir)

# ダウンロードファイルをチェック
os.listdir(data_dir/'gm-jpn-all_u_2_2')
# => ['trans_jpn.met',
# 'builtupp_jpn.prj',
# 'coastl_jpn.shx',
# 'polbnda_jpn.dbf',
# 'miscl_jpn.dbf',
# 'rstatp_jpn.dbf',
# 'raill_jpn.shp',
# 'coastl_jpn.shp',
# 'inwatera_jpn.prj',
# 'miscp_jpn.dbf',
# 'raill_jpn.shx', ...

GeoPandas

GeoPandas 0.4.0 — GeoPandas 0.4.0 documentation

GeoPandasというPythonのライブラリを使います。GeoPandasはPandasのデータ形式を活用しつつ、地図がPlotできるとても便利なライブラリです。

GeoPandas 0.4.0 — GeoPandas 0.4.0 documentation

※ GeoPandasの基本的な説明は僕が勝手に大尊敬しているsinhrksさんのBlogのエントリがオススメです

sinhrks.hatenablog.com

GeoPandasを使ってどの都市から編集されているかをPlotしてみます。

# 先程のクエリの結果をDataFrameに読み込み
df = pd.read_gbq(query, project_id=os.environ.get('PROJECT_ID'), dialect='standard')

# 日本地図の読み込み
gdf_jpn = gpd.read_file(str(data_dir/'gm-jpn/gm-jpn-all_u_2_2/polbnda_jpn.shp'))

# GeoDataFrame形式に変換
df['coordinates'] = list(zip(df['longitude'].astype(float), df['latitude'].astype(float)))
df['coordinates'] = df['coordinates'].apply(Point)
gdf = gpd.GeoDataFrame(df, geometry='coordinates')

# 日本地図に重ね合わねて表示
fig, axes = plt.subplots(figsize=(24, 8))
gdf_jpn.plot(color='C7', axes=axes)
gdf.plot(axes=axes, alpha=.3, color='C0', markersize=100)
plt.xlim(129, 145)
plt.ylim(30, 45)
plt.grid(False)

f:id:mergyi:20180916021718p:plain

人工が多い都市圏に固まっているのは予想通りですが、予想以上に全国的に満遍なく分布していることがわかりました。また北陸地方が少し多めなことも気になりました。もしかしたらJAISTなどがあることが影響しているのかもしれません。

BigQuery GIS

さて、実は本当に最近なのですが、BigQuery GISというサービス/機能群がGoogleから発表されて、今GeoPandasでやったようなことを BigQueryを書くだけでできる ようになりました。簡単に使い方を説明しながら地図上にPlotしてみたいと思います。

www.youtube.com

Geo Vizにアクセス

可視化には、BigQuery GISの中のGeo Vizを使います。まず下のURLにアクセスします。

BigQuery Geo Viz

Queryの入力

下記画像のSelect Dataに下記のようなQueryを入力します。

f:id:mergyi:20180916024041p:plain

WITH safe AS (
  SELECT
    NET.SAFE_IP_FROM_STRING(contributor_ip) AS ip,  
    LENGTH(NET.SAFE_IP_FROM_STRING(contributor_ip)) AS len 
  FROM
      `publicdata.samples.wikipedia` 
  WHERE contributor_ip IS NOT NULL
),
wiki_ips AS (
  SELECT 
    NET.IPV4_TO_INT64(ip) AS clientIpNum,
    TRUNC(NET.IPV4_TO_INT64(ip)/(256*256)) AS classB_
  FROM safe 
  WHERE ip IS NOT NULL 
    AND len = 4
), 
address AS (
  SELECT MAX(b.longitude) AS longitude, MAX(b.latitude) AS latitude, b.city, COUNT(*) AS count
  FROM wiki_ips  
  LEFT JOIN `fh-bigquery.geocode.geolite_city_bq_b2b`AS b ON wiki_ips.classB_ = b.classB
  WHERE wiki_ips.clientIpNum BETWEEN b.startIpNum AND b.endIpNum
  GROUP BY b.city 
  HAVING MAX(b.country) = 'JP'
)
SELECT 
    ST_GeogPoint(SAFE_CAST(longitude AS FLOAT64), SAFE_CAST(latitude AS FLOAT64))  AS wkt, count
FROM address 
ORDER BY count DESC
;

Queryは先程のものに少しだけクエリを追加しています。ポイントは ST_GeogPoint(SAFE_CAST(longitude AS FLOAT64), SAFE_CAST(latitude AS FLOAT64)) です。ここでPOINT型に変換しています。BigQueryのWebConsoleで見ると変換されていることがわかります。

f:id:mergyi:20180916024228p:plain

他の型や変換方法についてはDocumentに記載されています。 Geography Functions in Standard SQL | BigQuery | Google Cloud

カラムを選択

POINT型が入っているカラムを選びます

f:id:mergyi:20180916024314p:plain

スタイルを選択します

マーカーの色や大きさ、透過度などを設定できます

f:id:mergyi:20180916024537p:plain

完成

BigQueryを書いて少し設定するだけで地理情報を可視化することができました。 f:id:mergyi:20180916030135g:plain

ちなみに世界版の同様に簡単に作れます。こちらは次の項目で書きますが2000件の上限に引っかかってしまい全ての点をPlotできていない点にご注意ください。アメリカの西側がやけに少ないのはそれが原因だと思います。

f:id:mergyi:20180916031246p:plain

制限

Geo Vizには幾つかの制限があるのでご注意ください。例えば実際に業務で使う際には下記のような制限が影響してくる可能性が高いです。

Geo Vizには上限2000個の結果しかPlotできません
Geo Vizはシェアや保存はできません

こういった場合は、GeoPandasなど他の選択肢を検討するのが良さそうです。

まとめ＆感想

BigQueryに入っているIPデータからGeoPandasやBigQuery GISを使ってお手軽に地理空間分析を行えることがわかりました
Geo VizにはJupyterNotebookのExtensionもあって、Jupyterからも表示できるみたいなのでこちらも触ってみたいと思いました
Geo Vizは注意点としてはまだBeta版としては挙動が怪しい場面もちょくちょくありますが、その際はリロードしたら直りましたので参考にしてください！

フリーランチ食べたい

No Free Lunch in ML and Life. Pythonや機械学習のことを書きます。

BigQuery GIS/GeoPandasを使ってお手軽にIPアドレスで地理空間分析を行う「Wikipediaはどこから編集されている？」