せっかくだから俺は現在地から小選挙区を検索できるようにするぜ

ネット選挙が解禁されて早一年、公式の政党関係者のネット発信が、支持者のSAN値をカツオ節のごとく削り続けるエクストリームゲームと化した今日この頃、皆様いかがおすごしでしょうか。
さて、今回は現在地から小選挙区の情報を取得する方法を考えてみます。

##結果
http://needtec.sakura.ne.jp/analyze_election/page/ElectionArea/shuin_47

このページでは、県を選択するか、現在地を取得することで、小選挙区の候補が表示されます。
さらに小選挙区を選択することで、選挙区のだいたいの位置と候補者の一覧が表示されます。

##ソースコード
https://github.com/mima3/analyze_election

依存しているライブラリ
lxml-3.4.0-py2.7-freebsd-9.1-RELEASE-p15-amd64.egg
rdp-0.5-py2.7.egg
numpy-1.9.1-py2.7-freebsd-9.1-RELEASE-p15-amd64.egg
sympy-0.7.5-py2.7.egg
Beaker-1.6.4-py2.7.egg

sympyは0.7.5でないと動きません

##使用データ
国土数値情報　行政区域データ
http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03.html

衆議院小選挙区選出議員の選挙区（都道府県別）
http://www.soumu.go.jp/senkyo/senkyo_s/news/senkyo/shu_kuwari/

上記をCSVに手入力で変換したデータ
https://github.com/mima3/analyze_election/blob/master/election_area.csv

朝日新聞デジタル>2014衆院選>候補者
http://www.asahi.com/senkyo/sousenkyo47/kouho/

上記をCSVに変換したデータ
https://github.com/mima3/analyze_election/blob/master/script/candidate_shuin_47.csv

##データ作成の手順

# データベースの作成
python create_db.py election.sqlite

# 国土数値情報　行政区域データのインポート　数時間で終わるよ！
python import_administrative_boundary.py election.sqlite area\N03-14_140401.xml

# 行政区域データをsympyのPolygonに変換する 24時間ぐらいで終わるよ！
python convert_poly.py election.sqlite

# 小選挙区の情報を登録
python import_election_area.py election.sqlite election_area.csv

# 小選挙区の候補者情報を登録
python import_candidate.py election.sqlite shuin_47 script\candidate_shuin_47.csv
 

##解説
###行政区域データ
国土数値情報には行政区域データがXMLで提供されています。
これを用いれば、行政区割りをGoogleMap上に表示することができます。
ただし、このデータは莫大な大きさなので、取扱いにはいくつか注意する点があります。

####大きなサイズXMLを解析する
大きなサイズのXMLを解析する場合、XMLファイルを文字に変換してパースするというやりかたをすると、メモリ使用量が激増し、処理しきれなくなります。

そのため、lxml.etree.iterparseを用いて、順次処理するようにします。
実際の処理をみてみましょう。

election_db.py

    def ImportAdministrativeBoundary(self, xml):
        f = None
        contents = None
        namespaces = {
            'ksj': 'http://nlftp.mlit.go.jp/ksj/schemas/ksj-app',
            'gml': 'http://www.opengis.net/gml/3.2',
            'xlink': 'http://www.w3.org/1999/xlink',
            'xsi': 'http://www.w3.org/2001/XMLSchema-instance'
        }
        self._conn.execute('begin')

        print ('admins....')
        context = etree.iterparse(xml, events=('end',), tag='{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}AdministrativeBoundary')
        for event, admin in context:
            adminId = admin.get('{http://www.opengis.net/gml/3.2}id')
            print (adminId)
            bounds = admin.find('ksj:bounds', namespaces=namespaces).get('{http://www.w3.org/1999/xlink}href')[1:]
            prefectureName = admin.find('ksj:prefectureName', namespaces=namespaces).text
            subPrefectureName = admin.find('ksj:subPrefectureName', namespaces=namespaces).text
            countyName = admin.find('ksj:countyName', namespaces=namespaces).text
            cityName = admin.find('ksj:cityName', namespaces=namespaces).text
            areaCode = admin.find('ksj:administrativeAreaCode', namespaces=namespaces).text
            sql = '''INSERT INTO administrative_boundary
                     (gml_id, bounds, prefecture_name, sub_prefecture_name, county_name, city_name, area_code)
                     VALUES(?, ?, ?, ?, ?, ?, ?);'''
            self._conn.execute(sql, [adminId, bounds, prefectureName, subPrefectureName, countyName, cityName, areaCode ])

            admin.clear()
            # Also eliminate now-empty references from the root node to <Title> 
            while admin.getprevious() is not None:
                del admin.getparent()[0]
        del context

上記のコードは国土数値情報を解析してDBに格納している処理の一部です。
この例では、etree.iterparseを用いて、「AdministrativeBoundary」タグが検知されるたびに、処理を行っています。

詳細は、以下のページを参考にすると良いです。
lxml を使用して Python での XML 構文解析をハイパフォーマンスにする
http://www.ibm.com/developerworks/jp/xml/library/x-hiperfparse/

####座標情報を間引く
国土数値データは正確な座標情報が格納されているため、サイズが大きくなります。
今回は、ざっくりした座標がわかればいいので、その情報を間引きます。

線を間引くアルゴリズムとしては，Ramer-Douglas-Peuckerというアルゴリズムがあります。
このアルゴリズムの説明は下記のページを見ると良いでしょう。

[Mathematica] 折れ線を間引く
http://www.330k.info/essay/oresenwomabiku

Pythonの Ramer-Douglas-Peucker アルゴリズムのモジュールとしては下記が利用できます。
https://pypi.python.org/pypi/rdp

###小選挙区の情報登録
総務省は例によって例のごとく、小選挙区の情報をPDFでしか公開しておりません。
なので、頑張って手入力するか、別の方法をためす必要があります。

また、小選挙区の情報が意外とファジーな感じになっています。
たとえば、岩手県の第一区と第二区を見てみましょう。
http://www.soumu.go.jp/main_content/000240041.pdf

第一区　盛岡市（本庁管内、盛岡市役所青山支所管内、盛岡市役所簗川支所管内、盛岡市役所太田支所管内、盛岡市役所支所管内、盛岡市都南総合支所管内）

第二区　盛岡市（第一区に属しない区域）

とあります。
当然の権利のように、本庁管内がどこを表すかの情報は、ここで提供されていません。
そのため、今回は、盛岡市だった場合は第一区と第二区が小選挙区の候補としてあがるようにしました。

このようにデータ解析をさせないという強い意志を持った総務省に対して呪詛を吐きながら作成したCSVが以下の通りになります。
https://github.com/mima3/analyze_election/blob/master/election_area.csv

また、このあたりで心がおれたので、候補者の情報については朝日新聞から引っ張ってきています。
下記のスクリプトで、それっぽいデータが抜けるので、必要に応じて手で直し使用しています。

###現在地から小選挙区を取得する。

ブラウザで現在地を取得する。

ブラウザで現在地を使用するには、navigator.geolocationを使用します。
以下のようなコードで経度、緯度が取得できるでしょう。

if (!navigator.geolocation) {
  //Geolocation APIを利用できない環境向けの処理
  alert('GeolocateAPIが使用できません');
}
navigator.geolocation.getCurrentPosition(function(position) {
  console.log(position)
}

なお、IEの場合、現在位置が、ひどくずれます。
これは、Microsoftが使用している位置用法のデータベースがChromeやFirefoxが利用している位置情報のデータベースより精度が悪いためだと考えられます。（おれは悪くねー！おれは悪くねー！）
https://social.technet.microsoft.com/Forums/ie/en-US/aea4db4e-0720-44fe-a9b8-09917e345080/geolocation-wrong-in-ie9-but-not-other-browsers

特定の座標が行政区域に含まれているか調べる。

特定の座標が行政区域に含まれているか調べるには、多角形中に点が含まれているか調べるのと同じ意味になります。
Pythonの場合、SymPy で 点の内外判定を行えば 基本的 には良いです。

[Python]SymPy で 点の内外判定
http://d.hatena.ne.jp/factal/20121013

基本的にはと言ったのは、この処理がくっそ遅いからです。
数万ある行政区画に対して一々やっていたのでは時間がかかりすぎてしょうがないです。

このために、ここでは２つの対策を施しました。

１つ目、点の内外判定を行う前に、対象の多角形の数を絞る。
これは、多角形の頂点と指定の点の距離が、ある範囲内にある場合のみ検索対象としています。

具体的には下記のようなコードになります。

election_db.py

    def GetPos(self, lat, long):
        """
        現在の経度緯度に紐づくcurveのデータを取得
        """
        m =0.005
        while 1:
            rows = self._getCurveId(lat, long, m).fetchall()
            if rows:
              break
            m = m * 2
            if m > 0.1:
              return None

        dict = {}
        pt = Point(lat, long)

        for r in rows:
            key = r[0]
            ret = self._isCrossCurveId(pt, key)
            if ret:
                return ret
        return None

２つ目として、あらかじめPolygonのオブジェクトを作成し、データベースに登録しておくことです。
具体的には以下のような実装になります。

    def ConvertPoly(self):
        """
        curveテーブルからpolygonの作成
        """
        #gc.enable()
        #gc.set_debug(gc.DEBUG_LEAK)
        sql = '''DELETE FROM polygon'''
        self._conn.execute(sql)

        sql = '''select curve_id,lat,lng from curve order by curve_id'''
        rows = self._conn.execute(sql)
        dict = {}
        for r in rows:
            key = r[0]
            if dict.get(key) is None:
                dict[key] = []
            dict[key].append((r[1], r[2]))
        i = 0
        self.Commit()

        self._conn.execute('begin')
        for key in dict:
            print (key + ":" + str(i))
            #b = len(gc.get_objects())
            self._createPoly(key, dict[key])
            i = i + 1
            #gc.collect()
            #print str(b) + ":" + str(len(gc.get_objects()))
            if i % 100 == 0:
                clear_cache()
                self.Commit()

    def _createPoly(self, key, list):
        poly = Polygon(*list)
        obj = pickle.dumps(poly)
        sql = '''INSERT INTO polygon
                       (curve_id, object)
                     VALUES(?, ?);'''
        self._conn.execute(sql, [key, obj ])
        del poly
        del obj

オブジェクトのシリアライズはpickle.dumpsを用いて行います。
また、sympyはキャッシュ機構を使用しているので、大量にオブジェクトを作成するとメモリがすぐに数GByte消費されます。
それを避けるため、100個作成するため、clear_chacheを用いてキャッシュを削除しています。

また、最新のsympyのPolygonで"Polygon has intersecting sides."を例外を出力する条件が0.7.5と0.7.6で異なっています。
0.7.6では条件が厳しくなり、このデータでは正常に作成できないです。
そのため、sympyは0.7.5を使用してください。