地域メッシュとは

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

１次メッシュ: 一辺約80km.
２次メッシュ: 一辺約10km. １次メッシュを縦横８分割したもの。
３次メッシュ: 一辺約1km. ２次メッシュを縦横10分割したもの。
４次メッシュ: 一辺約500m. ３次メッシュを縦横2分割したもの。
５次メッシュ: 一辺約250m. ４次メッシュを縦横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: 静的な地図上にプロット

データの取得

今回は、１次メッシュ 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から始まる５次メッシュをすべて含んでいます。
KEYCODE: 当該メッシュのコード全体
MESHX_ID: X次メッシュ部分のコード。これらを１〜５次までくっつけたものがKEY_CODE
OBJ_ID: 行番号で特に意味はない（と思う）
geometry: 当該メッシュの地理情報。ポリゴンとして定義されている。

ポリゴンの集約

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

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

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...

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

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

地理情報の可視化には、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 の第１引数で描画する先（matplotlib の Axes）を指定します。そこで pyplot.subplots 関数からAxes を作成します。
第２引数は geopandas.GeoSeries

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

import mplleaflet

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

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

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

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

mesh2 = aggregate_mesh(x, 2)

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

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

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

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

１つ目については、例えば数千〜万の散布図を与えたりするとなかなか描画できません。
代替案として、静的な地図に描画することを試みます。いくつかのオプションがあるようですが、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 は使用前に１度実行が必要で、同じ情報を繰り返し取得するのを避けるため、取得した地図画像のキャッシュを作成します。アクセスを繰り返してOpenStreetMapに迷惑をかけないための配慮のようです。
shapely はポリゴンをはじめ地理情報の基礎的な要素の機能を提供しているライブラリで、geopandasも依存しています。ここでは、座標系を変換したポリゴンを定義し直すのに用いています。

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

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

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

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

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

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

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))

４次メッシュは３次メッシュを縦横２分割しており、南西・南東・北西・北東にそれぞれ 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))

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

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