Pythonとfoliumと作物統計データを使って北海道の市町村別に大根の生産量マップを作成する

できたもの

北海道の市町村別の大根生産量マップ（2017年）:

上の図は画像ですが、実際に生成されるファイルは動的に動かせるマップを含んだhtmlファイルです。マーカーをクリックすると、市町村の名前・大根生産量・ランキングもポップアップされます。
こんな:

文字つぶれてますが、最初にポップアップを開いたところが大樹町で、次にポップアップを開いたところが帯広市です。

はじめに

今年の3月までIT企業にいましたが、6月から北海道の大樹町という地域で大根の選別をしています。人生いろいろあります。

さて、都道府県別に見ると、大根の生産量第1位は北海道だそうです。まあ広いしな。では、北海道の市町村別に見ると、最も大根を生産しているのはどの市町村でしょうか。いま私が住まう大樹町はランキングの上位に入るでしょうか。folium というPythonのモジュールを使って、色付けしたマップを作成し可視化してみます。

ここでは北海道の市町村ごとの大根の生産量マップを作製しますが、都道府県と作物の種類の組み合わせを変えれば、お好きな都道府県・作物で同様のマップを作成することができます。このページの末あたり にいくつかの作例をあわせて付けています。夏休みの自由研究にいかがでしょう。

※ folium等の基本的な情報はあまり書いていません。foliumのクイックスタート や 末尾の参考 に挙げたページなんかをあわせてご覧ください。

利用言語・ライブラリ

• Python: 3.6.6
• folium: 0.9.1
• Pandas: 0.23.1
• json: 0.9.1
• あとExcel（Libreとかでも問題ないです）

ざっくり手順

• 市町村区の境界データと市町村別の大根生産量データをダウンロードする。
• 大根生産量データをExcelで加工する。
• ここからPython:
  • 境界データをjsonモジュールで読み込む。
  • 加工した生産量データをpandasで読み込む。
  • foliumを使って両者を紐づけて地図として表示する。
  • （マーカーを付けたい場合は以下も続けて）
  • 市町村の役場の位置データをダウンロードする。
  • pandasで読み込んで、生産量データと紐づける。
  • foliumを使って地図にマーカーを落とす。

使用したデータ

北海道の市町村区の境界データ と 北海道の市町村別の大根生産量データ と 日本の市町村の役場位置データ を用いました。それぞれ、行政区域コード をキーとして紐づけます¹。

北海道の市町村区の境界データ

入手先: 国土数値情報 行政区域データ → リンク

• 国土交通省が公開している。
• 上のリンク先からzipをダウンロードする。zipを解凍した中に含まれるGeoJSON形式の.geojsonファイルを用いる。
• （メモ）たどりかた:
  • 国土数値情報 ダウンロードサービスのトップページ
  • → "2. 政策区域/行政区域"
  • → "ダウンロードするデータの選択"のセクションで北海道を選んで"次へ"
  • → 好きな年度のデータを選んで"次へ"
  • → アンケートに答えて利用約款に同意したらデータをダウンロードできる

北海道の市町村別の大根生産量データ

入手先: 作物統計調査 → リンク

• 農林水産省が公開している。
• 上のリンク先からExcelファイルをダウンロードする。
• （メモ）たどりかた:
  • 統計情報トップページ
  • → 農林水産統計データ/品目別分類/野菜
  • → 生産/作物統計
  • → 調査一覧/作況調査（野菜）
  • → 市町村別データ/平成29年
  • → 野菜調査/8-1 野菜（根菜類）/北海道

日本の市町村の役場位置データ

入手先: 日本の地方自治体一覧オープンデータ → リンク

• 有志により作成され公開されている。
• 上記リンク先のGithubから"localgovjp-utf8.csv"をダウンロードする。

データの加工

大根統計データを整形する

Excelを使ってざざっとやります。
なお、農林水産省のサイトから入手したデータでは、大根は、春だいこん・夏だいこん・秋冬だいこんの3種に分類されています。今回は夏だいこんの数値のみ用います（大樹町が夏大根しかやってないから）。あしからず。
また、作物統計内で、収穫量が「…」または「x」とされている市町村は欠測（NaN）扱いで行ごと消します²。

整形前イメージ:

整形後イメージ:

整形手順:

1. （今回は要らないので）D列以降 "夏だいこん"に関する3列以外を削除する。
2. ヘッダ4行を削除する。
3. 最下部にある"利用上の注意"を削除する。
4. 列名を1行に収める。こんな感じ:
5. 作付面積列・収穫量列・出荷量列の書式設定をユーザ定義から標準に変更する。こうしておかないと、csvにエクスポートしたときに余計なスペースやダブルクォーテーションが入ってしまう。
6. 都道府県コード列と市町村コード列を結合して、"行政区域コード"と名前を付ける。
   単純に結合すると都道府県コード（北海道は"01"）の0が落ちてしまうので、
   =TEXT(都道府県コードのセル,"0#")&市町村コードのセル という関数など使う。
7. 作付面積と収穫量の値が"…"または"x"である行を削除する。
8. 残った市町村を対象に、大根の収穫量のランキングを計算する。あとで使う。=RANK(収穫量のセル,収穫量のセル全範囲絶対参照) など。"ランキング"と名前を付ける。
9. UTF-8でcsvにエクスポートする。

北海道の市町村の境界を表示した地図を出力する

ここからPythonです。
ひとまず、大根のデータを紐づける前に、北海道の市町村区の境界を表示してみます。

hokkaido.py

import folium
import json

# 市町村の境界データの読み込み
f = open(ダウンロードしたgeojsonファイルのパス, 'r', encoding='utf-8_sig')
geojson = json.load(f)
f.close()

# マップの作製
map_center = [43.06417, 141.34694] #札幌
m1 = folium.Map(location=map_center, tiles='cartodbpositron', zoom_start=7)
folium.Choropleth(
    geo_data=geojson,
).add_to(m1)
folium.LayerControl().add_to(m1)
m1.save('hokkaido.html')

あら簡単。
セーブされたhtmlファイルを開くとこんな感じです。地図の移動・ズームなんかも行えます:

実は問題点あり（後述）

上の図、一見のところよさげに見えますが、実は問題点があります。
札幌市が、中央区・北区・…・清田区と、区で分割されている点です。
それの何が問題なのか と 解決方法は 後の節 に書きます。

これ以降、その問題を解消して、札幌市を札幌市として一区画にしている前提とします。
→
問題解決前（札幌市が区で分かれていた） → 問題解決後（市に統合した）

北海道の市町村別に大根のデータを可視化する

すでにデータの加工は済んでるので、紐づけるキー（行政区域コード）をfoliumさんに教えてあげるだけでよしなにしてくれます。

make_prod_map.py

import folium
import json
import pandas as pd

# 市町村の境界データの読み込み
f = open(ダウンロードしたgeojsonファイルのパス, 'r', encoding='utf-8_sig')
geojson = json.load(f)
f.close()

# 大根データの読み込み
prod_df = pd.read_csv(Excelからエクスポートしたcsvのパス, dtype='object') # dtypeをobjectにしておかないと行政区域コード先頭の0が落ちる
prod_df['収穫量'] = prod_df['収穫量'].astype('float')

# マップの作製
map_center = [43.06417, 141.34694] #札幌
m2 = folium.Map(location=map_center, tiles='cartodbpositron', zoom_start=7)
folium.Choropleth(
    name='choropleth',
    geo_data=geojson,
    data=prod_df,
    columns=['行政区域コード', '収穫量'], # dataとして指定したDataframeの キーを第1引数に、表示したい値を第2引数に指定する
    key_on='feature.properties.N03_007', # 境界データ側の行政区域コードにあたる項目のキー名
    fill_color='BuGn',
    nan_fill_color='darkgray', # "geo_data"側（境界データ）には存在するが、"data"側（大根データ）に存在しない市町村はNaN扱いになる。NaNと分かるように灰色で塗る。
    fill_opacity=0.8,
    nan_fill_opacity=0.8,
    line_opacity=0.2,
    legend_name='2017年夏だいこん収穫量（t） ... 灰色は生産なし'
).add_to(m2)
folium.LayerControl().add_to(m2)
m2.save('product_amount_map.html')

セーブされたhtmlファイルを開くとこんな感じです:

んん……なんかすかすか……。北海道には178の市町村があるそうですが、作物統計のデータ上、夏だいこんを生産しているのはそのうち29だけみたいです。すかすかもやむなしか……。生産量が多い地域についても、あんまり法則性は見られないですね。（ちょっとデータを加工したらすこしマシな結果が出ました。次の次の節にて）

地図をリッチにする - マーカーとポップアップの設定

気を取り直してマップを改良しよう。北海道の土地勘が薄いため、上のマップだとどの区域がどの市町村なのかいまいちピンと来ません。
（大根を生産している）市町村にマーカーを配置して、市町村名や生産量の数値をポップアップで表示できるようにします。マーカーを配置する位置は、各市町村の中心的な役場の位置にしましょう。

make_prod_map_with_marker.py

import folium
import json
import pandas as pd

# 市町村の境界データの読み込み
f = open(ダウンロードしたgeojsonファイルのパス, 'r', encoding='utf-8_sig')
geojson = json.load(f)
f.close()

# 大根データの読み込み
prod_df = pd.read_csv(Excelからエクスポートしたcsvのパス, dtype='object') # dtypeをobjectにしておかないと行政区域コード先頭の0が落ちる
prod_df['収穫量'] = prod_df['収穫量'].astype('float')

# 各行政区域の役場データの読み込み、大根データとの結合
local_df = pd.read_csv(ダウンロードした役場位置データのcsvのパス, dtype='object')
local_df.rename(columns={'cid': '行政区域コード'}, inplace=True)
local_df['行政区域コード'] = local_df['行政区域コード'].str.zfill(5) # 行政区域コードの先頭の0が元データでは落ちているので追加する。
prod_df = pd.merge(prod_df, local_df, on='行政区域コード', how='inner') # 生産がNaNの市町村にはマーカーを落とさないという方針なので、inner joinにする。

# マップの作製
map_center = [43.06417, 141.34694] #札幌
m3 = folium.Map(location=map_center, tiles='cartodbpositron', zoom_start=7)
folium.Choropleth(
    name='choropleth',
    geo_data=geojson,
    data=prod_df,
    columns=['行政区域コード', '収穫量'], # dataとして指定したDataframeの キーを第1引数に、表示したい値を第2引数に指定する
    key_on='feature.properties.N03_007', # 境界データ側の行政区域コードにあたる項目のキー名
    fill_color='BuGn',
    nan_fill_color='darkgray',　# "geo_data"側（境界データ）には存在するが、"data"側（大根データ）に存在しない市町村はNaN扱いになる。NaNと分かるように灰色で塗る。
    fill_opacity=0.8,
    nan_fill_opacity=0.8,
    line_opacity=0.2,
    legend_name='2017年夏だいこん収穫量（t） ... 灰色は生産なし'
).add_to(m3)
folium.LayerControl().add_to(m3)

# マーカーの配置
city_count = len(prod_df)
for index, item in prod_df.iterrows():
    # マーカーをクリックすると現れるポップアップの設定
    popup_text = "市町村名: " + item['city'] + "<br/>"
    popup_text = popup_text + "収穫量（t）: " + str(item['収穫量']) + "<br/>"
    popup_text = popup_text + "ランキング: " + str(item['ランキング']) + '位 / ' + str(city_count) + '位'
    html = folium.Html(popup_text, script=True)
    # マーカーの設定
    folium.Marker(
        location=[item['lat'], item['lng']],
        popup=folium.Popup(html, max_width=150, min_width=150),
        icon=folium.Icon(color='green')
    ).add_to(m3)

m3.save('product_amount_map_with_marker.html')

セーブされたhtmlファイルを開くとこんな感じです（冒頭の図の再掲）:

マーカーのサイズがちょっと大きい気がするけど、これを変更する方法は分からなかった。
マーカーをクリックするとこんな感じのポップアップが表示されます:

大樹町は10位かー。ちなみに1位は帯広市でした。

htmlファイルをダウンロードできるようにしました。動かしてみたい方はご自由にどうぞ（ブラウザさえあれば動かせます）: リンク

ちょっと加工する - 生産力だとどうだろう。

大根の生産量を可視化してみましたが、まあなんか、目立って面白い結果にはなりませんでした。

ちょっとデータを加工してみましょう。
作物統計には、収穫量のほかに、作付面積のデータも入っています。当然、作付面積が広いほど生産量も増えるでしょうから、生産量を作付面積で割った単位面積あたりの生産量（いわば生産力と言えるでしょう）を見てみます。

Excelに戻って、生産力を計算してcsvをエクスポートしなおします。使うカラムや色遣いなどをちょっと変えて、make_prod_map.py をもう1度動かします。

こうなりました:

十勝平野や根釧台地の市町村は生産力高めな傾向が見えますね。道東は北海道の中でも夏季が冷涼な地域なので夏だいこんの生育に向いているのかも知れません。生産力でランキングをつけると大樹町は3位だったのでなんか嬉しかったです。

あ、あと余談ですけど、大根の旬って秋から冬な感じがするじゃないですか。夏に大根を作るってもしやハウスとか人工環境下でやってるのかな、と思っていたんですが、農家さんちに来てみたら普通に露地栽培でした。涼しいから自然環境下で育つんすね。

東京23区と政令指定都市はもう1手順 要る

途中の節で挙げた問題点の補足説明です。

なにが問題なのか

実は、市町村の境界データと 作物統計のデータとは、1対1で紐づいてはいません。
たとえば、北海道で言うと、境界データは札幌市の中央区・北区・…・清田区と区で分割されていますが、作物統計のデータは札幌市でひとくくりになっています。札幌市の行政区域コード（01100）と、各区の行政区域コード（たとえば中央区は01101）は別物なので上手く対応付けられません。

全データを確認したわけではないですが、東京都の特別区（23区） と 政令指定都市 について、

• 農林水産省の作物統計は、特別区または市単位でデータを集計している。
• 国土交通省の行政区域データは、区単位で境界データを整備している。

ようです。東京都または政令指定都市を有する都道府県で同じことをしたい方は、以下もご参照ください。

暫定的解決方法

以下の手順で解決しました。やや迂遠かも知れません。

札幌市の境界データを入手する

"札幌市"の境界データがあればなんとかなるでしょう。探してみたら見つかりました:

入手先: 北海道札幌市 (01201A1968) | 歴史的行政区域データセットβ版 → リンク

• ROIS-DS人文学オープンデータ共同利用センター というところが公開している。
• 上のリンク先から、好きな基準年月日のTopoJSONファイルをダウンロードする。最新のものを選べば問題ないでしょう。
• （メモ）たどりかた
  • Geoshapeリポジトリ トップページ
  • → 歴史的行政区域データセットβ版
  • → 市区町村ID インデックス
  • → 01201A1968 北海道札幌市

入手したTopoJSONファイルをGeoJSON形式に変換する

TopoJSONはGeoJSONの友達のようです。お互いに相互に変換できます。国土交通省の行政区域データに組み込んで使いたいので、いまダウンロードしたTopoJSONファイルをGeoJSON形式に変換します。変換用のツールもいろいろあるようですが、ちゃちゃっとやりたいのでウェブサービスを利用しました。

使ったウェブサービス: MyGeoData Converter

ユーザ登録は不要で利用できます。TopoJSONファイルを選択して"Convert Now!"するだけなので、特に操作で困ることはないでしょう。
1つ注意する点として、国土交通省の境界データの測地系（Coordinate System）はJGD2011なので、Output Parameterをそれに合わせておきましょう:

ちなみに、無料会員だと月に3ファイルまでしか変換できないらしいのでお気をつけください。

札幌市の境界データを結合する

手順としては2ステップ:

1. 国土交通省の境界データから要らない区を削除する。今回は、札幌市を構成するすべての区（それらの区の行政区域コードを別途に入手する必要があります。ググればすぐ出ます）。
2. 上のサイトから取得したTopoJSONファイルには行政区域コードが入っていませんでした。行政区域コードを追加した上で（今回は札幌市の行政区域コードを別途に入手する必要があります。ググればすぐ出ます）、国土交通省の境界データに追加します。

こんな感じのコードになります:

# 1. 国土交通省の境界データから要らない区を削除する。
exclude_codes = ['01101', '01102', '01103', '01104', '01105', '01106', '01107', '01108', '01109', '01110'] # 札幌市に属する区の行政区域コード
# print('削除前: ' + str(len(geojson['features'])))
new_features = [feature for feature in geojson['features'] if feature['properties']['N03_007'] not in exclude_codes]
geojson['features'].clear()
geojson['features'] = new_features
# print('削除後: ' + str(len(geojson['features'])))

# 2. 行政区域コードを追加した上で、国土交通省の境界データに追加する。
s = open(上のウェブサービスで変換したGeoJSONファイルのパス, 'r', encoding='utf-8_sig')
sapporo_geojson = json.load(s)
s.close()

for feature in sapporo_geojson['features']:
    feature['properties']['N03_007'] = '01100' # 札幌市の行政区域コード
    geojson['features'].append(feature)

このスニペットを、ここまでのコード（hokkaido.py, make_prod_map.py, make_prod_map_with_marker.py) の "市町村の境界データの読み込み" のあとに追加すれば、それより下のコードを変更することなく、札幌市を1つの区域として扱うことができます。

作例集

※ 市町村の境界データの利用手続き上、300x400より大きい地図を1サイトに6枚以上掲載する場合は 申請が要る みたいなので、画像小さめ。

北海道、ばれいしょ生産量マップ（2017年）:

... 北海道で広く作られてることが分かります。十勝平野と網走・知床あたりが強い。

北海道、甜菜生産量マップ（2017年）:

... ばれいしょほど広い地域では作られてないようだ。しかしこいつも十勝平野と網走・知床あたりが強い。

新潟、米（水稲）生産量マップ（2017年）:

... 3月まで住んでたゆかりある新潟。さすがに全域で生産がある。新潟市が1位なんだ。ちょっと意外。

広島、そば生産量マップ（2017年）:

... 出身県。ひい婆さんが住まう北広島はそばが名産で、ランキング1位じゃないかな、と思ったけど、庄原市に負けてた。海に近いとあんまり作ってないのも分かる。

おわりに

Pythonとfoliumを使って北海道の市町村別に大根の生産量/生産力マップを作成しました。

行政区域と作物統計の組み合わせを変えれば、どの都道府県でも、市町村別/作物別に、生産量や生産力の色付きマップを作製することができるでしょう。
というか、行政区域コードが記載されているデータであれば、作物統計以外でも紐づけは簡単ですので、色んなパターンのマップを作製できそうです。面白いマップを作ったらぜひ教えてください。

参考

• Folium 0.9.1 公式ドキュメント
  • モジュールのリファレンス
  • クイックスタート
• akatak’s blog: 【珈琲たいむ】Pythonでコロプレス図を描く
• サッカーを通じて観戦力と分析力を上げていくブログ: 【データ可視化×Python】foliumを使って福島県の統計データを地図上に可視化してみる（コロプレス図）

1. リンク先の記述のとおり、呼称に揺れがあるようですが、ここでは"行政区域コード"という名前で統一します。 ↩

2. NaN扱いした市町村において、大根を「生産していない」とまで言い切れるかは微妙。「…」は「事実不詳又は調査を欠くもの」を意味して、「x」は「個人又は法人その他の団体に関する秘密を保護するため、統計数値を公表しないもの」を意味するため。http://www.maff.go.jp/j/tokei/kouhyou/sakumotu/sakkyou_yasai/gaiyou/index.html#13 ↩