イントロダクション

このチュートリアルでは、とあるデータサイエンティストの典型的な1日の過ごし方をご案内しましょう。まず地理空間のデータセットを入手し、不要なものを整理し、補強し、可視化します。使用するツールはPython、BeautifulSoup、pandasとNominatimライブラリ、そして地理情報システムの組織で広く使われているオープンソースの地図ソフトウェア、 QGIS です。

データセットになるのは、全米UFO情報センター（NUFORC）の このページ に掲載されているアメリカ全土のUFO目撃情報です。目標は、過去12カ月間に目撃されたUFOの地図を可視化することです。可視化によりデータセットをはっきりと示し、調査して、目撃されたとされるUFOの行動をよりよく理解することができます。可視化は地図作成プログラム内で行われます。QGISは地理空間データの手軽な試験的分析に特に向いているのです。また、可視化したものはビデオやアニメーションに出力し、同プログラムを使う他のユーザに共有することもできます。

最初のタスク： Webからデータを抽出する

UFO目撃情報のwebサイトを訪問すると、データはかなり構造化されたフォーマットだということが分かります。月毎のページがあり、そこにはUFO目撃情報の行と、時間、都市名、州、形状と説明のカラムがあります。ですから過去12カ月分のページをダウンロードして、ページのHTMLからデータを抽出すればいいのです。BeautifulSoupはDOMをパースして簡単にクエリできる適切なライブラリで、これを使ってメインページから過去12カ月のデータへのリンクを入手し、その後情報を抽出します。

コード

from datetime import datetime
from itertools import chain
from urllib import urlopen

from bs4 import BeautifulSoup
import pandas as pd

base_url = "http://www.nuforc.org/webreports/"
index_url = "http://www.nuforc.org/webreports/ndxevent.html"

# get the index page
source = BeautifulSoup(urlopen(index_url), "html5lib")
# get all the links in the index page
monthly_urls = map(lambda x: (x.text, base_url + x['href']),source('a'))
# get the last 12 links that have a text like 06/2015
last_year_urls = filter(lambda x: can_cast_as_dt(x[0], "%m/%Y"), monthly_urls)[:12]
# extract the data from each monthly page and flatten the lists of tuples
last_year_ufos = list(chain(*map(lambda x: get_data_from_url(x[1]), monthly_urls)))
# initialize a pandas DataFrame with the list of tuples
ufos_df = pd.DataFrame(last_year_ufos, columns=["start","city","state","shape","duration_description"])

def can_cast_as_dt(dateStr, fmt):
try:
datetime.strptime(dateStr, fmt)
return True
except ValueError:
return False

上記のコードで注目してほしいのは、htmllibを使っている点です。どうもNUFORCのHTMLはフォームが完全ではなく、HTMLの終了タグが1個ではなく2個あったりします。BeautifulSoupのデフォルトのパーサはこういうケースに対処できないので、引数”htmllib”を挿入しています。下記に、月毎のページからデータを抽出する関数を記載しました。ページの構造をうまく活用しBeautifulSoupを使っています。それから、日付文字列をPythonのdatetimeに変換します。このページの日付のフォーマットは「01/21/15」という形式になっています。これは実際のところあいまいな表記ではあるのですが、このチュートリアルの範囲で言えば今世紀の出来事と考えて間違いないため、2015年が正しい年だと予測してかまいません。

def get_data_from_url(url):
print "Processing {}".format(url)
data = []
source = BeautifulSoup(urlopen(url), "html5lib")
for row in source('tr'):
if not row('td'):
continue # header row
row_data = row('td')
# parse the datetime from the string
datetime = parse_dt(row_data[0].text)
city = row_data[1].text
state = row_data[2].text
shape = row_data[3].text
duration = row_data[4].text
data.append((datetime, city, state, shape, duration))
return data

def parse_dt(dateStr):
# the data in the website comes in two different formats, try both
for fmt in ["%m/%d/%y %H:%M", "%m/%d/%y"]:
try:
return datetime.strptime(dateStr, fmt)
except ValueError:
continue

2つ目のタスク： データを構造化して補強する

さて、webから抽出したUFO目撃情報を基にしたDataFrameオブジェクト(ufos_df)が手に入りました。説明の情報については無視することにします。この情報はあまり構造化されていないため、私たちの用途には使えないからです。ここで2つ問題があります。都市と州の情報から、その都市が世界のどこにあるのかを推測するには十分なのですが、位置の座標がありません。加えて、UFOが目撃された継続時間の情報はあまり正確なルールに従っていません。例えば”25 minutes “や、”25+ Min”、そして”約25mins”などとフォーマットされています。ですから、このデータを使えるようにするため、ここで余分なものを整理し、標準化する必要があります。

UFO目撃の継続時間と終了時刻を抽出する

import re

# extract the duration in seconds
ufos_df["duration_secs"] = ufos_df["duration_description"].apply(infer_duration_in_seconds)
# now we can infer the end time of the UFO sighting as well
# which will be useful for the animation later
ufos_df["end"] = ufos_df.apply(lambda x:x["start"] + timedelta(seconds=x["duration_secs"]),axis=1)

# function that infers the duration from the text
def infer_duration_in_seconds(text):
# try different regexps to extract the total seconds
metric_text = ["second","s","Second","segundo","minute","m","min","Minute","hour","h","Hour"]
metric_seconds = [1,1,1,1,60,60,60,3600,3600,3600]
for metric,mult in zip(metric_text, metric_seconds):
regex = "\s*(\d+)\+?\s*{}s?".format(metric)
res = re.findall(regex,text)
if len(res)>0:
return int(float(res[0]) * mult)
return None

上記のコードでは推測された継続時間の単位を”秒”で返す関数を各行に適用しました。通常、継続時間のレポートは大体次のようになっています。”（約、多分、大体、などの修飾子）＋数値＋任意でプラスのマーク＋（複数形だったりそうでなかったり、標準化されていないフォーマットでの時間の単位）”。

これをパースするには、正規表現が使えます。様々な単位を試し、時間の単位にマッチする秒数を乗じた数値を正規表現でキャプチャして返します。興味深いのは、正規表現は前のステップにおいてHTMLをパースするのにいいツールではなかったのですが、文字列から継続時間の情報を抽出するには適切なツールだといえることです。任意の文字の存在と、特殊な文字クラス（例えば桁など）が、正規表現に向いているからです。データの中のおよそ85～90％の継続時間の記述が上記の簡単な関数によってパースできる、つまりNoneは返さないということです。

位置をジオコーディング（地理座標を付加）する

都市と州の情報をより有益な地理的情報に変換するために、ジオコーディングサービスを利用する必要があります。APIが分かりやすいので、Nominatimを選びました。

from nominatim import Nominatim
geolocator = Nominatim()
geolocator.query("Houston, TX")

クエリオブジェクトは可能性のある結果のリストを返します。各結果は”lat（緯度）”と”lon（経度）”のカラムを含んだディクショナリです。もしクエリによって位置を特定できなかった場合、最初の結果の座標を採用するかNoneを返すことにします。シンプルなPythonの辞書を使って結果をキャッシュしておくと、”Houston, TX”（ヒューストン, テキサス）などというクエリを何度も送らずに済みます。NominatimはREST APIを毎回クエリする必要があり、少々時間がかかるからです。チュートリアルがごちゃごちゃしないよう、実際に使われているコードはここではお見せしていませんが、pandasの”apply”機能を使っています。都市・州のカラムを入力値として、座標を表すカラムを新しく作成するために用いています。

最終的には、DataFrameのすべての行に”lat”と”lon”列が入ります。もちろん、これらの座標はUFO目撃がレポートされた都市を反映していますが、目撃者の正確な位置ではありません。そもそも、そんなに正確な情報は記録されていないのですから。きっと宇宙人は気にしないでしょう。

データの統合

更に処理を進めるために、全てのデータをCSVフォーマットに出力します。

# Note: dropna will drop any columns with None values, which is desirable
ufos_df[["start","end","lon","lat","shape"]].dropna().to_csv("ufo_data.csv",index=False, encoding="utf-8")

3つ目のタスク： QGISを使ってデータを可視化する

QGISの設定

手元には、整理されたUFO目撃データのCSVファイルがあります。QGISを使ってこのデータを可視化するには、地図の背景を表示させるための OpenLayers プラグインと、指定した時間内でデータを処理するためのバージョン2.0以上の TimeManager プライグインを追加でインストールする必要があります。新しいプラグインをインストールするには、ツールバーから”プラグイン” > “プラグインの管理とインストール”を選択し、”OpenLavers”と”TimeManager”を検索します。もし可能であれば、バージョン2.9以上のQGISを使ってください。TimeManagerは、元々 Anita Graser によって開発されたプラグインですが、数年をかけて、私、 carolinux も、リファクタリングや維持管理、新しい機能の開発など、大いに貢献してきました。

プラグインがインストールできたら、次は、ツールバーから”Web” > “OpenLayers Plugin”を選択し、低階層から地図を選んでください。このチュートリアルでは、OpenStreetMapsを使用します。もし、プログラムの一番下にTimeManagerのスライダが表示されていなければ、ツールバーの”プラグイン” > “TimeManger”を選択し、”Toggle visibility”をクリックしてください。

データの読み込み

ツールバーから”レイヤ” > “レイヤの追加” > “デリミティッドテキストレイヤの追加”を選択し、CSVファイルを選択した上で、下記のように設定します。

そして、TimeManagerにある”Settings”ボタンをクリックし、”Add Layer”をクリックします。下記のように設定を行ったら、”OK”をクリックします。

これで、UFO目撃データの時間をスライドさせたり、タイムフレームの範囲を変更したりできるようになりました。例えば、タイムフレームの範囲を2時間に設定し、現時刻から2時間前までの間に報告されたUFOの目撃情報を地図上に表示させることができます。スライダを使えば、データセットを検証することができ、同時刻または別の場所での群発など、特定の状況下で報告された複数の目撃情報を見ることもできます。どうやら、2015年3月6日の夜は、未知なる地球外生命体でにぎわっていたようですね。

ここで使用したQGISのクールな機能が、可視化するためのカスタムのSVGmaker（エイリアンの顔）です。QGIS上でカスタムのSVGmakerを設定する方法については、 こちら の説明書をご覧ください。

スタイルをデータに対応させる

ここまで順調に進んできましたが、QGISのとても便利な機能をまだ使っていません。それは、データからスタイルを定義する機能です（多くの可視化ライブラリで使われています）。この機能を活用して、UFOの形状情報を可視化していきます。QGISで”ルール”を定義することは可能で、もし、ある点の属性が条件を満たすのであれば、レンダリングの際には、特定のマーカースタイルが使われます。

これらのルールを設定するには、ufo-dataレイヤ上で右クリックし、プロパティを選んだら、Styleの項目で”Rule-based”を選択します。プラス記号をクリックすると新しいルールを追加することができます。条件を満たす属性は、Filterのテキストボックスで定義されます。円形物体として報告された全てのUFO目撃情報に対するスタイルの定義は、円形メーカーとしてフィルタ、選択されるように”shape=’circle'”と書きます。以下は、私が定義したルールです。

出力と共有

TimeManagerプラグインの”Export Video”をクリックすることで、データのフレームを出力することができます。それぞれのフレームは時間間隔ごとのデータのスナップショットになっています。これらのフレームは、ImageMagickやffmpegを使った動画への変換といったコマンドラインユーティリティでアニメーション化することができます。TimeManagerプラグインの設定にある”Do not export empty frames”にチェックを入れることで、UFOが現れていないフレームの出力を避けることができ、”Display frame start time on map”にチェックを入れると、現在時刻を表示するテキストラベルを得ることができます。これは、カスタマイズが可能です。以下の画像は、3月のとある2週間のUFO目撃情報をアニメーション化したものです。GIFが機能しない場合は、こちらの 直リンク から読み込んでみてください。

QGISプロジェクトのファイルそのものを共有することもできます。必要な情報全てがファイルに含まれているので、XMLフォーマットでこのプロジェクトを再現できます。更に、データソースに関連するパスもあるので、CSVやqgsファイルが含まれたフォルダを、チーム内やソースコントロールで共有することができるのです。

まとめと改善点

今回行ったのは、およそ構造化されたデータをWebサイトからいくつか抽出したものをデータ変換、整理し、QGISに読み込んで時間ごとに区切り、そしてそれをアニメーションまたはビデオとして出力するということでした。地理空間のデータを可視化するための機能を持ったプログラムを使い、時空間の側面を操作すれば、時間を節約することができます。地図上にポイントを置いたり、時間を操作したりするようなロジックがすでに実装されているので便利です。代替の提案があれば、是非コメントを寄せてください。2週間のスナップショットの整理されたデータセットをいじってみたいという方は、 こちら からダウンロードしてください。

このワークフローで明らかに改善しうる点は、パイプラインでのデータ変換の完全自動化です。これはデータサイエンティストが心に留めておくべきことです。この内容については以前に 記事 にしたので、読んでみてください。最初のバージョンである現状においては、QGISを起動したり、クリックしたり、出力したフレームをアニメーションやビデオに変換したり……といった手動操作が必要になっています。もしQGISプロジェクト生成をデータパイプラインの一環と捉えるのであれば（つまり、何度も行わなければならないルーチンの一部と捉える場合です。私の環境においてはそうでした）、QGISをスクリプトしたり、事前に作成したQGISプロジェクトのテンプレートからほとんどの設定を読み込んだりしなければなりません。幸運なことに、これはPyQGISバインディングを使ってPythonで対応することがきますし、全てプログラムで処理することができます。次のバージョンでは、 TimeManagerプラグインからワンクリックでビデオやアニメーションを出力できないかというリクエストを受けましたので、次回リリースするバージョンでこの機能が追加されることを期待していてください。

謝辞

TimeManagerプラグインの共同制作者であるAnita Graser、特定のデータセットにおいて可視化のアイデアを提供してくれたJulijana Busic、そして校正を行ってくれたYannis Chatzimichosに心からお礼を言いたいと思います。