51
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

posted at

updated at

日本の地域メッシュを地図上に可視化して理解する

地域メッシュとは

地域メッシュは、緯度・経度に基づき地域をほぼ同じ形状の網目に分割したものです。1つのメッシュはおおよそ正方形で、様々なサイズで定義されています。
頻繁に使うものを簡単にまとめると、

  • 1次メッシュ: 一辺約80km.
  • 2次メッシュ: 一辺約10km. 1次メッシュを縦横8分割したもの。
  • 3次メッシュ: 一辺約1km. 2次メッシュを縦横10分割したもの。
  • 4次メッシュ: 一辺約500m. 3次メッシュを縦横2分割したもの。
  • 5次メッシュ: 一辺約250m. 4次メッシュを縦横2分割したもの。

参考:総務省統計局 | 地域メッシュ統計について

メッシュにはコード体系が決まっており、緯度経度から計算することが可能です(参考:標準メッシュの体系とコード)。ですが、各メッシュコードの定義は、データとしても提供されています(e-Stat | 地図で見る統計(統計GIS) )。この記事では、これを用いて地域メッシュの可視化することで理解したいと思います。

ソフトウェア・ライブラリ

Python(Anaconda/miniconda distribution)を用います。必要なライブラリは下記コマンドからインストール可能です。

conda install -c conda-forge jupyter geopandas descartes shapely \
  matplotlib pip requests pillow chardet mplleaflet && \
  pip install tilemapbase

下記のライブラリは、今回の目的に特徴的です。

  • geopandas: 地理情報データの処理に用いる
  • mplleaflet: インタラクティブな地図上にプロット
  • tilemapbase: 静的な地図上にプロット

データの取得

今回は、1次メッシュ 5339のデータを使用します。ファイルは、こちらのリンク から直接ダウンロードするか、あるいはe-Stat のページから "M5339" のシェイプファイルをダウンロードします。Zipファイルになっているので、これを展開し、下記の様な構成になっていれば準備完了です。

└── QDDSWQ5339
    ├── MESH05339.dbf
    ├── MESH05339.prj
    ├── MESH05339.shp
    └── MESH05339.shx

データの読み込み

geopandas ライブラリを使うとシェイプファイルの読み込みは非常に楽です。

import geopandas as gpd

x = gpd.read_file("QDDSWQ5339/MESH05339.shp")
print(x.shape)
x.head()

(100800, 8)
     KEY_CODE MESH1_ID MESH2_ID MESH3_ID MESH4_ID MESH5_ID  OBJ_ID  \
0  5339000011     5339       00       00        1        1       1   
1  5339000012     5339       00       00        1        2       2   
2  5339000013     5339       00       00        1        3       3   
3  5339000014     5339       00       00        1        4       4   
4  5339000021     5339       00       00        2        1       5   

                                            geometry  
0  POLYGON ((139.00312 35.33333, 139.00000 35.333...  
1  POLYGON ((139.00625 35.33333, 139.00312 35.333...  
2  POLYGON ((139.00312 35.33542, 139.00000 35.335...  
3  POLYGON ((139.00625 35.33542, 139.00312 35.335...  
4  POLYGON ((139.00937 35.33333, 139.00625 35.333... 
  • データは5339から始まる5次メッシュをすべて含んでいます。
  • KEYCODE: 当該メッシュのコード全体
  • MESHX_ID: X次メッシュ部分のコード。これらを1〜5次までくっつけたものがKEY_CODE
  • OBJ_ID: 行番号で特に意味はない(と思う)
  • geometry: 当該メッシュの地理情報。ポリゴンとして定義されている。

ポリゴンの集約

今回のデータは、5次メッシュのデータなので、たとえば "5339" 1次メッシュの全体そのものはどこにも定義されていません。粒度の粗いメッシュについては、5次のポリゴンを集約することで得られます。これには、geopandas.GeoDataFrame の、dissolveメソッドが便利です。

次の関数は、5次メッシュのデータからより低次のデータへの集約を行います。

def aggregate_mesh(x, level):
    tmp = x.copy()
    code_len = [4, 6, 8, 9, 10][level-1]
    tmp["key"] = tmp["KEY_CODE"].str.slice(0, code_len)
    tmp = tmp[["key", "geometry"]]
    tmp = tmp.dissolve(by="key")
    return tmp
  • code_len: 指定された粒度のコードの長さ
  • tmp["key"] に指定された粒度のコードを指定
  • .dissolve でポリゴンを集約
  • .dissolve は、地理情報については和集合を取ります。今回は使っていませんが、それ以外の列がある場合には groupby.aggregate のような計算を同時に行ってくれます
# aggregate to mesh level 1
mesh1 = aggregate_mesh(x, 1)
mesh1

      geometry
key                                                    
5339  POLYGON ((139.29062 35.33333, 139.28750 35.333...

1次メッシュに集約した結果です。今回は "5339" の1次メッシュのみを使っているので、1行に集約されます。

メッシュ地理情報の可視化(インタラクティブな地図)

地理情報の可視化には、geopandas.plotting モジュールが便利な関数を提供しています。特に、ポリゴンについては、plot_polygon_collection が便利です。文字通り、複数のポリゴンを可視化する関数です。

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
plot_polygon_collection(ax, mesh1.geometry, edgecolor="#121212", facecolor="#333333")
  • plot_polygon_collection の第1引数で描画する先(matplotlibAxes)を指定します。そこで pyplot.subplots 関数からAxes を作成します。
  • 第2引数は geopandas.GeoSeries

mesh1.png

結果は上のようなグラフになります。これだとどこにあるかはわかりません。これを地図上にプロットするために、mplleaflet ライブラリを使います。このライブラリは、matplotlib のグラフをインタラクティブな地図上に描画する機能を提供しています。

import mplleaflet

mplleaflet.display(fig)  # jupyterの場合
# mplleaflet.show(fig)   # その他

leaflet1.png

Jupyterで実行すると上図のように地図がアウトプットに表示されます。通常の地図アプリのように拡大したり移動したりが可能です。
.show の方を使うとブラウザの別ページが開きそこに同様の地図が表示されます。

"5339" メッシュは東京から南関東の一部を含むメッシュということがわかります。ぴったり正方形になっていないのは、海の部分の5次メッシュが定義されていないためでしょう。

同様に、2次メッシュを表示してみます。

mesh2 = aggregate_mesh(x, 2)

fig, ax = plt.subplots()
plot_polygon_collection(ax, mesh2.geometry, edgecolor="#121212", facecolor="#333333")
mplleaflet.display(fig)

leaflet2.png

メッシュ地理情報の可視化(静的な地図)

mplleaflet はとても便利な機能ですが、いくつか弱点もあります。

  • プロットする量が多くなると急激に重くなる
  • 現時点では、文字の描画に対応していない

1つ目については、例えば数千〜万の散布図を与えたりするとなかなか描画できません。
代替案として、静的な地図に描画することを試みます。いくつかのオプションがあるようですが、TileMapBase ライブラリを利用します。このパッケージはOpenStreetMapの地図を画像として取得して、それを背景に描画します。この上に描画を重ねることで matplotlib の機能を自由に使うことができます。
弱点としては、mplleaflet ほどシンタックスが単純ではなく、座標変換をユーザー側が行う必要があります。

少し長くなりますが、地図上にポリゴンを描画し、その中央に文字列を描画する関数です。

import tilemapbase
tilemapbase.init(create=True)

from shapely import geometry

def plot_on_map_with_text(geoseries, texts):
    # find the graph range
    rect = geoseries.total_bounds
    edgex = rect[2] - rect[0]
    edgey = rect[3] - rect[1]
    extent = tilemapbase.Extent.from_lonlat(
        rect[0]-edgex*0.3, rect[2]+edgex*0.3,
        rect[1]-edgey*0.3, rect[3]+edgey*0.3)
    extent = extent.to_aspect(1.0)

    fig, ax = plt.subplots(figsize=(8, 8), dpi=100)
    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)

    t = tilemapbase.tiles.build_OSM()
    plotter = tilemapbase.Plotter(extent, t, width=600)
    plotter.plot(ax, t)

    polygons = []
    centers = []
    bounds = geoseries.bounds
    for i in range(len(bounds)):
        # convert to the plottable scale
        minx, miny = tilemapbase.project(bounds['minx'][i], bounds['miny'][i])
        maxx, maxy = tilemapbase.project(bounds['maxx'][i], bounds['maxy'][i])
        polygons.append(
            geometry.box(minx, miny, maxx, maxy))
        centers.append([(minx + maxx)/2.0, (miny + maxy)/2.0])

    polygons = gpd.GeoSeries(polygons)
    plot_polygon_collection(ax, polygons,
                            edgecolor='#121212', facecolor='#000000', alpha=0.4)
    for center, txt in zip(centers, texts):
        ax.text(center[0], center[1], txt, fontdict={'color':'lightblue'})
    return fig, ax
  • tilemapbase.init は使用前に1度実行が必要で、同じ情報を繰り返し取得するのを避けるため、取得した地図画像のキャッシュを作成します。アクセスを繰り返してOpenStreetMapに迷惑をかけないための配慮のようです。
  • shapely はポリゴンをはじめ地理情報の基礎的な要素の機能を提供しているライブラリで、geopandasも依存しています。ここでは、座標系を変換したポリゴンを定義し直すのに用いています。

2次メッシュを、その下2桁のコードともにプロットします。

plot_on_map_with_text(mesh2.geometry, mesh2.index.str.slice(4, 6))

static2.png

こう見ると、メッシュコードのルールがわかります。南西端を (0, 0)として、北へ行くと1つめ、東へ行くと2つめの数字が大きくなっていく仕組みです。2次メッシュは1次メッシュを8分割するので、0から7までの数字が現れます。

3次以降も同様に描画できます。ただし、グラフが細かくなりすぎるので、特定のエリアに限定します。

mesh3 = aggregate_mesh(x[x.MESH2_ID == "45"], 3)
plot_on_map_with_text(mesh3.geometry, mesh3.index.str.slice(6, 8))

static3.png

コードのルールは変わらず、今度は10分割されています。

mesh4 = aggregate_mesh(x[(x.MESH2_ID == "45") & (x.MESH3_ID == "09")], 4)
plot_on_map_with_text(mesh4.geometry, mesh4.index.str.slice(8, 9))

static4.png
4次メッシュは3次メッシュを縦横2分割しており、南西・南東・北西・北東にそれぞれ 1, 2, 3, 4 のコードを付与しています。

mesh5 = aggregate_mesh(x[(x.MESH2_ID == "45") & (x.MESH3_ID == "09")], 5)
plot_on_map_with_text(mesh5.geometry, mesh5.index.str.slice(8, 10))

static5.png
4つの4次メッシュを5次に分割して表示しています。5次メッシュのコード規則は4次メッシュと同様です。

Register as a new user and use Qiita more conveniently

  1. You can follow users and tags
  2. you can stock useful information
  3. You can make editorial suggestions for articles
What you can do with signing up
51
Help us understand the problem. What are the problem?