フリーランチ食べたい

機械学習、Python、ソフトウェアエンジニアリング、プロダクティビティなど

matsumotoryさん、kwappaさんの公開ポートフォリオレビューを受けました/自分なりにアウトプットすることのメリットを整理してみた

「エンジニアの日々のアウトプット」に関するForkwellさんのイベントに参加してきました。

  • ForkwellのポートフォリオのぼりーさんのクラウドPodcastのmatsumotoryさんゲスト回を聞いて興味を持ち使い始めました。
  • ただ、埋めやすいところを埋めてGithubを連携させただけで、ちゃんと使いこなせていたかは疑問でした
  • そんなときに、このイベントを前職の同期から聞き、開催場所も前職のオフィスだったので久々に遊びに行きたい気持ちもあって応募しました。

forkwell.connpass.com

ポートフォリオを(強制的に)ちゃんと埋めました

イベントページに

ゲストのお二人が参加者のポートフォリオを抽選で3名レビューいたします!

との記載があったので、せっかくなので…ということで応募したら、なんと選ばれました! 85人の前でレビューされる、というプレッシャーもあり、週末の半日くらいを使って必死に埋めました。ときには追い込まれることも必要だと思いました。

発表: matsumotoryさん・kwappaさん

speakerdeck.com

最初はmatsumotoryの発表でした。

  • 自分自身の経験とアウトプットを続けていて得たもの
  • アウトプットの方法について挑戦を続けてきたこと
  • 様々な視点からアウトプットのメリット

などを伺いました。普段からmatsumotoryさんのBlogを拝読することがあるのですが、記事から伺える価値観や考え方をよりリアルに感じ、刺激になりました。

次はkwappaさんの発表でした。

  • 公開/非公開に関わらずアウトプットしていくことが大事
  • 技術メモや作業ログとして残していき、徐々にアウトプットの量や範囲を広げていく
  • 媒体はGithub・Qiita・はてなブログなんでも良いので自分に合っているものを使う

といった内容を伺いました。現実的でとても参考になるお話でした。

自分なりにアウトプットのメリットを整理してみる

お二人のお話を聞きながら、自分の頭の整理用にアウトプットすることのメリットをまとめていました。せっかくなのでここに書いておきたいと思います。それぞれの項目自体はお二人の発表にあったもので、自分は整理してみただけなのであしからず。

アウトプットすることのメリットは(もちろん色々な整理の仕方があると思いますが)「一時視点」「連続視点」で分けられるのではないかと考えました。

一時視点

f:id:mergyi:20181009232342p:plain

まず「一時視点」です。「一時視点」というのは、時間軸に関わらない、その時点時点でのメリットを指しています。この視点で見たときにはアウトプットすることには下記のメリットがあると思います。

1. 頭の中が整理される

「理解している」と思っていたことが、実際に書いてみるとわかっていなかったり、書いていく中で理解できてくることってたくさんありますよね。

2. 自己を客観的に見れる

1週間前に書いた記事や1年前に書いた文章を読み返すと、まるで他人の文章を読んでいるような気持ちになることがあり、そこで「自分はこういう考えの人間だったのか」ということに客観的に気づけた経験があります。

3. 他者を巻き込める

まだ自分はできていない気がしますが、matsumotoryさんの発表にあったので。 アウトプットから他人を巻き込んでアクションを起こして、そこからさらにアウトプットする、という流れができたら最高ですよね。

連続視点

f:id:mergyi:20181009232415p:plain

もう一つが「連続視点」です。「一時視点」と対象的に、時間が経過していくことで生まれるメリットを指していて、下記のような項目が当てはまると思いました。

1. 実績から未来が予想できる

今までの自分が書いた記事や発表を見返すことで、「自分がどのくらいのスピードで成長しているか」や「興味を持っている技術の方向性」が見えてきて、これから自分がどう進んでいくか、を客観的に予想することができます。

2. 自分の「ストーリー」がわかる

1.の発展形だと思いますが、過去から未来の自分の進捗を観察することで、「自分がどのような価値観で生きているのか」「何を目指して生きていくのか」がわかるようになるのではないでしょうか? 正直、自分はここまで達していないですが、アウトプットを続けていつか、それに気づけたら嬉しいです。

3. 実績を会社に評価してもらう

最後は一番現実的で馴染みがあるかもしれませんが、就職活動や転職活動で実績として認めてもらったり、会社内でも自分が挑戦したい案件があったときに、アウトプットをしていることで周りを説得しやすくなる、ということもあると思います。

特に印象に残った点

どちらかというと、自分の認識ではアウトプットすることのメリットは「一時視点」が中心で、特に「連続視点」の「実績から未来が予想できる」「自分のストーリーがわかる」ことは全然頭になかったので、本当に目からウロコが落ちるような気持ちでした。

その後公開レビューを受けました。

休憩を挟んでいよいよ公開レビューを受けました。自分の前が来る前に、matsumotoryさん、kwappaさんのポートフォリオを見る、というイベントがあり、お二人のポートフォリオに比べると自分のポートフォリオがあまりにしょぼいので、レビューが始まる直前は逃げ出したいような気持ちでした。笑

f:id:mergyi:20181009231047p:plain

しかし、実際に公開レビューが始まるとお二人とForkwellのakagawaさんが、かなり丁寧に、優しく、レビューしてくださったので、比較的にリラックスして話すことができました。 また内容も的確でためになることばかりでした。イベント後、matsumotoryさんがTwitterでコメントしてくださっていてとても嬉しかったです。

自分から質問することもできたので「Blogに対して反応が少ないと、Blogを続けていくのが辛くなっていくと思いますがどう継続するのがおすすめですか?」というような質問をしたところ、 matsumotoryから、matzさんの『俺を踏み台にしてくれ』という言葉を引用しながら、「誰か共有/拡散してくれる人を見つけるといいかもしれない。自分は最初のころmrubyのアウトプットをしていて、matzさんに共有/拡散してもらって人脈と可能性が広がっていった」というような回答をいただきました。なるほど〜と唸りっぱなしでした。

ちなみに僕の次に公開レビューを受けていたsachinさんは若干21歳で僕よりも遥かにエンジニア経験があってOSS活動などもバリバリされていて本当にすごいと思いました。心から尊敬します。

Forkwell Portfolio - エンジニア向けのポートフォリオサービス

さいごに

その後の懇親会でも同年代のエンジニアと何人か情報交換やお互いの仕事内容等を共有にとても刺激になりました。また、Forkwellのakagawaさんが大学の先輩であったことを知り(さらに大学時代も話したことがありました)衝撃を受けました。世界って本当に狭いですね。そういったこともあり、イベントの本筋以外でも非常に実りの多いイベント参加になりました。

matsumotoryさん、kwappaさん、Forkwellの皆さん、本当にありがとうございました。

ちなみに

全然関係ない&完全に偶然ですが、今日誕生日です。笑 28歳になりました。今年もチャレンジをやめない一年にしたいなと思ってます。 また今年一年よろしくお願いします。

ほしいものリスト

10月8日(月)につくばPythonもくもく会 No.1を開催します!

Pythonを使っている方&Pythonを勉強したい方一緒にもくもくしませんか?

  • 宣伝です。
  • これから、つくばでPythonもくもく会を定期的に開いていきます!
  • 第1回は10月8日(月)11:00~@筑波大学内体バチで行います。
  • 将来的にはPythonの情報共有ができたり、交流の輪を広げられるコミュニティにしていけたらいいなと妄想しています。
  • 参加費や条件などは無料なのでぜひぜひお気軽にお申し込みください。

詳細はconnpassページをご参照ください。

申し込みもこちらからできます。途中参加退場なども自由です!

tskubapy.connpass.com

動機など少しだけ

  • 一番直接的な動機はPyConJP 2018 に運営として参加したことです。
  • 今まで自分の会社以外でPythonエンジニアを触れ合う機会がほとんどなかったので、とても刺激的でした。
  • コミュニティがあると、インプットアウトプットの流れができやすく、1人で勉強しているのとは全く違う体験ができることを感じました。
  • つくばには、エンジニアのコミュニティは自分が知る限りそこまで多くなく、特にPythonのコミュニティは存在していない(と思う)ので、つくばでもそういう体験を作り出せたら嬉しいです。
  • また、PyCon JPの一部のPython BootCampにTAとして参加して、Pythonに興味を持っている/勉強したいと思っている人の多さに驚きました。そういった方々に勉強の場の提供や効率的な勉強法なども共有できたら良いなと思ってます。

blog.ikedaosushi.com

  • 最近、RのコミュニティであるTsukuba.R が復活したことも開催するきっかけになりました。 www.meetup.com

最後に

今回のもくもく会に参加したいけどできないという方、今後このもくもく会に参加してみたいという方、参加はするかわからないけどちょっと興味がある方、connpassに tsukuba.py というグループを作ったのでぜひぜひお気軽にご参加ください!Join us!

tskubapy.connpass.com

ISOに従っていないデータをpandas.to_datetime()すると500倍以上遅くなる可能性がある話とその対策

TL;DR

  • pandasの to_datetime メソッドはとても便利で、かなり乱暴にデータを突っ込んでもParseしてくれます
  • でもデータによってはparseに通常の30倍以上時間がかかる可能性があるので注意しましょう
  • ISO_8601の規格に従っていない場合はとりあえず format オプションをつけておくのが得策です。
  • コードはすべてGithubにあがってます

github.com

検証するデータ

こちらのKaggleのデータを使いたいと思います。なぜこのデータかというと実際に痛い目にあったからです。笑

Final project: predict future sales | Kaggle

df = pd.read_csv(Path.home()/'.kaggle/competitions/competitive-data-science-final-project/sales_train.csv.gz', compression='gzip')
print(df.shape) => (2935849, 6)

2935849件のデータがあります。

datetimeのparse

date カラムには下のような形式で日付のデータが入っています。

df['date'].head()
0    02.01.2013
1    03.01.2013
2    05.01.2013
3    06.01.2013
4    15.01.2013

このデータをdatetime型にparseしたい時には to_datetime を使います。

print(df['date'].dtype) # => object
print(pd.to_datetime(df['date']).dtype) # => datetime64[ns]

本当に何も考えなくても簡単にParseできましたね。

to_datetime はISO_8601に従っていない場合に急激に遅くなるケースがある

しかし、この便利な to_datetime には落とし穴があります。データがISO_8601に従っていないとき急激に速度が遅くなることです。ISO_8601とは日付と時刻を表す標準規格です。

参照: ISO 8601 - Wikipedia

今回使ったデータはISO_8601に従っていません。ちょっと時間を計測してみましょう。計測には、timeit decoraterを少しいじったものを使っています。

@timeit(repeat=3, number=10)
def convert_to_datetime(s):
    return pd.to_datetime()

ちょっと怖いので10000件だけで計測してみましょう。

df10000 = df.head(10000)
best_time, parsed_datetime = convert_to_datetime(df10000['date'])
# Best of 3 trials with 10 function calls per trial:
# Function `convert_to_datetime` ran in average of 0.883 seconds.

0.9秒弱かかっています。データは293万件以上あるので(もしO(n)だと仮定しても)Parseだけで263秒もかかってしまうことになります。最初これで全然処理が終わらなく「おかしいな…」と首を傾げていました。笑

解決策: format オプション

それではどうすればよいのかというと format オプションを使います。

@timeit(repeat=3, number=10)
def convert_to_datetime_with_format(s, format_='%d.%m.%Y'):
    return pd.to_datetime(s, format=format_)

format付きバージョンで時間を計測してみましょう。

best_time, parsed_datetime = convert_to_datetime_with_format(df10000['date'])
# Best of 3 trials with 10 function calls per trial:
# Function `convert_to_datetime_with_format` ran in average of 0.024 seconds.

0.02秒!圧倒的ではないか我軍は…!もしO(n)だと仮定すれば、7秒程度でParseが終わります。 このようにformatオプションをつけるだけでParse速度が変わることがわかりました。

念の為、計算量の増加具合(オーダー)も調べてみる

先程から「O(n)だと仮定すれば」と書いてきましたが、本当にそうなのか、調べてみましょう。

n_targets = [100, 1000, 10000, 100000, 1000000]
normal_best_times = []
format_best_times = []

for n in n_targets:
    target_df = df.head(n)

    best_time, parsed_datetime = convert_to_datetime_with_format(target_df['date'])
    format_best_times.append(best_time)

    best_time, parsed_datetime = convert_to_datetime(target_df['date'])
    normal_best_times.append(best_time)

plt.figure(figsize=(16, 8))
ax = plt.subplot()
ax.plot(n_targets, normal_best_times, color='C0', linewidth=5, label='normal')
ax.plot(n_targets, format_best_times, color='C1', linewidth=5, label='with_format')
plt.legend(fontsize=20)

f:id:mergyi:20181003030313p:plain きれいに線形(O(n))であることがわかりました。 それにしてもすごい差ですね。

実際に平均で何倍なのか見てみましょう。

ratio = (np.array(normal_best_times) / np.array(format_best_times)).mean()
print(f'with_format parse take only {ratio:.2f} times faster than normal parse time')
# with_format parse take only 34.00 times faster than normal parse time

ということで34倍早くなることがわかりました。

ちなみに

じゃあISO_8601に近い形式の場合はどのくらい計算時間が変わるのだろう?と疑問に思ったので調べてみました。 2018-10-03 という ISOに則った形式と 02.10.2018 というISOに則っていない形式でformatあり/なしでそれぞれ比べてみます。

target_s1 = pd.Series(['2018-10-03', '2018-10-02', '2018-10-01'] * 10000)
target_s2 = pd.Series(['03.10.2018', '', '01.10.2018'] * 10000)

best_time, parsed_datetime = convert_to_datetime(target_s1)
# Function `convert_to_datetime` ran in average of 0.005 seconds.

best_time, parsed_datetime = convert_to_datetime_with_format(target_s1, format_='%Y-%m-%d')
#  Function `convert_to_datetime_with_format` ran in average of 0.005 seconds.

best_time, parsed_datetime = convert_to_datetime(target_s2)
# Function `convert_to_datetime` ran in average of 2.884 seconds.

best_time, parsed_datetime = convert_to_datetime_with_format(target_s2, format_='%d.%m.%Y')
# Function `convert_to_datetime_with_format` ran in average of 0.072 seconds.

順番でいうと、 ISO = ISO with format > not ISO with format >> (越えられない壁) >> not ISO でした。ISOに従っていない形式だと ISOに従った形式の 576.8倍遅かった のでちょっと衝撃でした。formatオプションを使えば14.4倍にまで緩和できますので必ずformatオプションは使いましょう。

f:id:mergyi:20181003032859p:plain

まとめ

  • to_datetime メソッドはParseする日付データの形式によって34倍も遅くなることがわかりました。
  • それは format オプションを正しく指定することで高速化できます。
  • 遅いParseも早いParseも計算量はO(n)でした
  • 迷ったら to_datetime メソッドにはformat オプションをつけるようにしましょう。

pathlibで見るPythonの演算子オーバーロード活用

pathlibって便利ですよね

最近pathlibの便利さが様々なところで語られています。

pathlibの様々な機能は上記の記事やドキュメントを読んでいただければわかるので、今日はその1つに、Pythonオーバーロードを説明するのに良い機能があるので紹介したいと思います。

pathlibはこんな風にパスを書けます。

from pathlib import Path

etc_dir = Path('/etc')
init_dir = 'init.d'
print(etc_dir/init_dir/'reboot')
# => /etc/init.d/reboot

最初に見ると、ちょっとギョッとするのではないでしょうか?
pathlibでは上記のように、1 / 2 #=> 0.5 のように 変数 / 変数 という書き方でパスを書くことができます。慣れてくるとLinuxのパスを普通に書いているような気持ちで書けるので気持ち良いです。本来は除算の演算子である / でどうしてこんな風に書けるのでしょうか?答えはソースを読むとわかります。

https://github.com/python/cpython/blob/3.6/Lib/pathlib.py#L898

    def __truediv__(self, key):
        return self._make_child((key,))

この __truediv__ メソッドが答えです。

Python演算子オーバーロードできる

公式ドキュメントの 「3. Data model」を読むとPython演算子オーバーロードできることがわかります。

3. Data model — Python 3.7.0 documentation

The following methods can be defined to emulate numeric objects. Methods corresponding to operations that are not supported by the particular kind of number implemented (e.g., bitwise operations for non-integral numbers) should be left undefined.

object.__add__(self, other)  
object.__sub__(self, other)  
object.__mul__(self, other)  
object.__matmul__(self, other)  
object.__truediv__(self, other)
...

これで +-/ といった演算子を再定義できます。それではこの機能を使って擬似的にpathlibを実装してみると次のようになります。

class myPath():
    def __init__(self, path):
        self._path = path
    
    def __truediv__(self, add_path):
        self._path = self._path + '/' + add_path
        return self

    def __str__(self):
        return self._path

etc_dir = myPath('/etc')
init_dir = 'init.d'
print(etc_dir / init_dir / 'reboot')
# => /etc/init.d/reboot

無事、pathlibと同様の挙動が実装できました。 pathlibは演算子オーバーロードをとても上手く使った実装だと思います

他にも演算子オーバーロードはこんなところでも

pathlib以外で自分が良く使う演算子オーバーロードはnumpyの @ です。これを使うと機械学習の実装では欠かすことのできないドット積(内積)の計算を行うことができます。元々は dot メソッドを呼び出さなければいけなかったのですが、途中で演算子オーバーロードとして実装されたので簡単に書けるようになりました。

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
b = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

print(a * b)
# [[ 1  4  9]
#  [16 25 36]
#  [49 64 81]]

# 前は dot メソッドを呼び出す必要があった
print(a.dot(b))
# [[ 30  36  42]
#  [ 66  81  96]
#  [102 126 150]]

# 今は@演算子でdot積が計算できる
print(a @ b)
# [[ 30  36  42]
#  [ 66  81  96]
#  [102 126 150]]

演算子オーバーロード、とても便利ですね。

最後に

研究がESSSB 17th@ベルギー・ヘントでポスター発表されました!

ESSSBとは

正式名称はEuropean Symposium on Suicide & Suicidal Behaviourで、自殺行動/自殺予防について学際的な観点から研究された発表やポスターが多い学会です。今年はベルギーのヘントで開催されました。

www.esssb17.org

ポスター

タイトルは「Regional suicide rates and affecting factors in Japan: an ecological study by using Generalized Linear Model and Machine Learning」です。各都道府県ごとの自殺率と、失業率、中学/高校進学率、一人当たり災害復旧費などの関係について、状態空間モデルとRidge回帰を使って探索的/仮説検証的の両面から検証した研究になります。自分はモデルの設計/実装とグラフの作成を主に行いました。

f:id:mergyi:20180923174656p:plain

発表

発表はFirstAuthoのOeくんがしてくれました。この研究のアピールポイントは、この分野では今まで活用が少なかったMachineLearningを活用して仮説探索/検証を行った点だったのですが、現地でも好評で興味を持ってくれた方が多かったようです。

f:id:mergyi:20180923180657j:plain

さいごに

全国の実際の自殺率データを扱って考察を得るのは非常に興味深い体験でした。 今回はマクロな視点からの研究だったので、この研究で得られたことを活かしつつ細かい、具体的な見識を得ていきたいと感じました

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

TL;DR

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

データ

実際にはアクセスログなどを使うことになると思いますが、テスト用に使える手頃なアクセスログを持ち合わせてないので、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

下記のようなSQLfh-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版としては挙動が怪しい場面もちょくちょくありますが、その際はリロードしたら直りましたので参考にしてください!

PyCon JP 2018のトークのスケジューリングを行い、Blog投稿しました

ついに来週PyCon JP 2018です

今週buildersconに参加してきたのもあって、体が温まっている感じがします。笑 今回、PyCon JPには運営側で参加しますが、そもそも PyCon JPに参加するのがはじめてなので(!)そういう意味でもとてもワクワクしています。

トークのスケジューリングをBlogに投稿しました

自分の担当は以前Blogに投稿したように(PyCon JP2018のHPを作成しました。 - フリーランチ食べたい) システムなのですが、「タイムテーブルを作る人でが足りてない」とのことだったので助っ人として手伝わせていただきました。ただ、自分も当然もともとのシステムチームのタスクもあるので、そこまでガッツリ時間を取ることは難しそうで、そこから今回の最適化問題としてスケジューリングをするというアイディアが浮かびました。

pyconjp.blogspot.com

来週遊びに来てください

PyCon JP 2018のチケットは社会人分は完売しているのですが、学生分はまだ少数あるみたいなので、もし興味がある方いれば是非是非遊びに来てくださいね

pycon.jp