オープンデータでゴニョゴニョする準備として国土数値情報を使ってpythonのfoliumでコロプレス図書いて遊んでみました。

作業環境

Mac OS X 10.10 Yosemite

参考情報

やったこと

gdal, pyshpのインストール
国土数値情報のDLページからデータを行政区域のShapefileを取得
ShapefileをGeoJsonに変換
ポリゴン融合(各地物をまとめる)
コロプレス図色付け用データの取得
GeoJsonとweb上で取得したデータのヒモ付準備
MAPの表示

gdalのインストール

macの場合、以下のコマンドでインストールできます。

brew install gdal
pip install GDAL={インストールしたgdalのversion}
pip install pyshp

ここで、gdalのバージョンを合わせないとpythonパッケージのインストールでエラーが発生するので注意して下さい。
pyshpはShapefileを扱うために入れています。

行政区域のshapefileを取得

行政区域から都道府県の各市区町村の境界データがshapefileとして取得出来ます。
shapefileだと扱いづらいのでGeoJsonに変換しています。
とりあえず今回は島根県のデータを使ってみています。

ShapefileをGeoJsonに変換

以下のようなpython scriptで変換します。

import os
import osgeo.ogr
import shapefile

# Flashメッセージ用の関数
def msg(s): print (s)
def dashes(): msg(40*'-')
def msgt(s): dashes(); msg(s); dashes()
def msgx(s): dashes(); msg('ERROR'); msg(s); dashes(); sys.exit(0)

# shapefile->GeoJson関数
def convert_shp_to_geojson(shape_fname):
    if not os.path.isfile(shape_fname):
        msgx('File not found: %s' % shape_fname)

    # shapefileの読み込み
    try:
        reader = shapefile.Reader(shape_fname)
    except:
        msgx('Failed to read shapefile: %s' % shape_fname)

    fields = reader.fields[1:]
    field_names = [field[0].decode('shift-jis').encode('utf-8') for field in fields]

    geom_dict = {}
    uniq_key = 'N03_007'  # 地域コード
    for sr in reader.shapeRecords():
        encoded_sr_record = [r.decode('shift-jis').encode('utf-8') for r in sr.record]
        atr = dict(zip(field_names, encoded_sr_record))
        geom = sr.shape.__geo_interface__

        # ポリゴン融合
        # 諸島のように同一の行政区域でも領域が分断されているものを
        # type=MultiPolygonとしてまとめる
        if atr[uniq_key] in geom_dict.keys():
            if len(geom_dict[atr[uniq_key]]['geometry']['coordinates']) == 1:
                geom_dict[atr[uniq_key]]['geometry']['type'] = 'MultiPolygon'
                tmp = tuple(geom_dict[atr[uniq_key]]['geometry']['coordinates'])
                geom_dict[atr[uniq_key]]['geometry']['coordinates'] = []
                geom_dict[atr[uniq_key]]['geometry']['coordinates'].append(tmp)
            geom_dict[atr[uniq_key]]['geometry']['coordinates'].append(geom['coordinates'])
        else:
            geom['coordinates'] = list(geom['coordinates'])
            geom_dict[atr[uniq_key]] = {'type':'Feature', 'id':atr[uniq_key], 'geometry':geom, 'properties':atr}

    output_buffer = []
    output_buffer_apd = output_buffer.append
    for code, geo_data in geom_dict.items():
        geo_data['geometry']['coordinates'] = tuple(geo_data['geometry']['coordinates'])
        output_buffer_apd(geo_data)  # dict data

    # GeoJSON fileの書き出し
    out_fname = os.path.join(os.path.dirname(shape_fname), os.path.basename(shape_fname).replace('.shp', '.json'))

    with open(out_fname, "w") as geojson:
        json.dump({"type": "FeatureCollection", "features": output_buffer}, geojson, indent=2, ensure_ascii=False)
    msg('file written: %s' % out_fname)
    return out_fname

GIS系のデータの扱い方がよくわかっていないので、めちゃくちゃ効率の悪い書き方をしていそうな気がしています。。
(ポリゴン融合しつつShapefileからGeoJsonに一発で変換出来るツールくらいどこかにありそうですし...)

GeoJsonとして下のようなデータが入ったファイルが生成されます。

{
  "type": "FeatureCollection", 
  "features": [
    {
      "geometry": {
        "type": "Polygon", 
        "coordinates": [
          [
            [
              131.91386401800003, 
              34.57300200100008
            ],

ポリゴン融合

上のコード中に含まれていますので割愛します。
ポリゴン融合とはどういうものかについては、mk-modeさんのBLOGを参考にさせていただきました。

コロプレス図色付け用データの取得

yahoo!さんのgeocitiesの市区町村情報から拝借致しました。
ここでも島根県の市区町村データを取得しました。
データ取得のコードは下のようになります。

import pandas as pd

# ジオシティーズの対象のページ
url = 'http://www.geocities.jp/warera_tikyujin/country_information/japan/pref_information/town200704/preft200704_32.html'

df = pd.read_html(url, header=0, index_col=0)[1]

# データの整形
df.columns = ['city_name', 'population', 'household', 'ground', 'remarks']
df['city_name'] = df.index
# 単位の調整
df['population'] = df['population'].astype(np.int32)/10000

GeoJsonとweb上で取得したデータのヒモ付準備

web上で取得したデータには地域コードがないので付与してやります。
このコードをGeoJsonとデータの対応付けをする際のkeyとします。
keyを市区町村名としても良いかもしれませんが、日本語は何かと扱いが面倒くさいのでコードに変換しておくと無難かと思います。
標準地域コードは総務省統計局のページからcsvで取得可能です。

import os
import json

jsonfile = os.getcwd() + '/../city_code.json'
with open(jsonfile) as f:
    city_code_json = json.loads(f.read(), "utf-8")

def cityname_to_code(city_name):
    try:
        return city_code_json[city_name].encode('utf-8')
    except KeyError:
        return np.nan

# 日本語をキーにすると扱いづらいので、地域コードをキーにする
df.loc[:, 'code'] = df['city_name'].apply(cityname_to_code)
# ヒモ付が出来なかったデータを除外
df = df.dropna(subset=['code'])

MAP表示

全ての準備が整ったので引数を正しくあたえるだけで、あとはfoliumさんがよしなに処理してくれます。

import folium
from gdal_test import convert_shp_to_geojson

# Shapefileの読み込み
filepath = os.getcwd() + '/works/input/N03-140401_32_GML/N03-14_32_140401.shp'
geojson = convert_shp_to_geojson(filepath)

# map作成
m = folium.Map(location=[35.472297, 133.050499], zoom_start=7.0, tiles='Mapbox Bright')
# m.geo_json(geo_path=geojson)
m.geo_json(geo_path=geojson, data=df,
    columns=['code', 'population'],
    # key_on='feature.id',
    key_on='feature.id',
    threshold_scale=[1, 5, 10, 20, 50, 100],
    fill_color='BuPu', reset=True,
    legend_name='population (10k)')

m.create_map(path='city.html')

作業時点現在のfolium(version 0.1.3)の実装では、作成されたcity.htmlをchromeで表示させる場合、クロスドメインチェックにひっかかって表示出来ない可能性が高いです。
chromeで表示がうまくいかない場合はfirefoxなどの別のブラウザで確認してみて下さい。
うまく表示出来れば以下のようなMAPが表示されます。

今回でGIS系のデータの扱い方がなんとなくわかりました。
ただ上の内容、sinhrksさんの記事のほぼ丸パクリになってしまいました。。。
次はもう少し数学的要素のあるものを作ってみたいと思います。

国土数値情報とfoliumで市区町村区切りのコロプレス図を描く