opendata
python2.7

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

More than 3 years have passed since last update.

オープンデータでゴニョゴニョする準備として国土数値情報を使って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が表示されます。

島根画像.jpg

今回でGIS系のデータの扱い方がなんとなくわかりました。

ただ上の内容、sinhrksさんの記事のほぼ丸パクリになってしまいました。。。

次はもう少し数学的要素のあるものを作ってみたいと思います。