Zenn
🗾

地域メッシュコードからポリゴンを作成する

2023/11/08に公開
Python
GIS
GeoJSON
tech

はじめに

この記事では総務省統計局が公開する地域メッシュ統計の特質・沿革 に基づき、地域メッシュコードからポリゴンを作成するPythonプログラムについて説明します。

地域メッシュコードについて

地域メッシュコードとは日本地図に重畳されたメッシュの格子ごとの地域区画を一意に定める要素番号であり、4桁または6~11桁の数字から構成されています。
上記の資料のP18~20に定義に関する記載があるため、気になる方はそちらをお読みください。

ポリゴンについて

ポリゴンは主要な地理データ型の1つで、地球上のある領域を表します。
本記事では地理データの情報を持つjsonであるgeojson形式でポリゴンを作成します。
GeoJSONのポリゴンは、経度と緯度の座標の配列を含む「coordinates」キーを持つオブジェクトです。
この座標の配列は、ポリゴンの外周を定義するもので、必ず最初の座標と最後の座標が同じになるように閉じている必要があります。

プログラムはPyPIで公開されており、下記のコマンドでインストールできます

Terminal
pip install jismeshcode

また、GitHubからもソースコードを確認できます

対象読者

GIS学習者

地域メッシュコードからポリゴンを計算する

基準地域メッシュの区分方法について

基準地域メッシュは大きく分けて5つの区分があります。
総務省統計局が提供する地域メッシュ統計とはから区分に関する表を抜粋します。

区画の種類 区分方法 緯度の間隔 経度の間隔 一辺の長さ 地図との関係
第１次地域区画 全国の地域を偶数緯度及びその間隔(120分）を３等分した緯度における緯線並びに１度ごとの経線とによって分割してできる区域 40分 １度 約80km 20万分の１地勢図の１図葉の区画
第２次地域区画 第１次地域区画を緯線方向及び経線方向に８等分してできる区域 5分 7分30秒 約10km ２万５千分の１地勢図の１図葉の区画
基準地域メッシュ（第３次地域区画） 第２次地域区画を緯線方向及び経線方向に10等分してできる区域 30秒 45秒 約1km
２分の１地域メッシュ（分割地域メッシュ） 基準地域メッシュ（第３次地域区画）を緯線方向，経線方向に２等分してできる区域 15秒 22.5秒 約500m
４分の１地域メッシュ （分割地域メッシュ）２分の１地域メッシュを緯線方向、経線方向に２等分してできる区域 7.5秒 11.25秒 約250m
度分秒について

度分秒は地球上の位置を示すために使われる、球体の位置を表す単位です。
緯度と経度を度分秒から十進数に変換するには、次の数式を使用します。
十進数 = 度 + (分 / 60) + (秒 / 3600)

次の図の各区内の数字は第1次地域区画のメッシュコードを表しています。
日本の国土にかかる第１次地域区画（世界測地系）

第1次地域区画では5438(4桁)、第2次地域区画では543823(6桁)…といったように、一片の長さが短くなる(区画が小さくなる)ごとに地域メッシュコードの長さは長くなります。
しかし、地域メッシュコードのつけ方に癖があるため、計算は少々複雑です。
それぞれの区画について該当地域を計算する方法を説明します。

第1次地域区画(1次メッシュ)の計算

それでは1次メッシュから取り組みます。
1次メッシュの地域メッシュコードの付け方は下図を参考にしました。

4桁の数字の左から2番目までの2桁が南端緯度に1.5を掛けた数字、4桁の数字の左から3,4番目の2桁が西端経度の下2桁の数字を意味しています。
また、緯度の間隔が40度、経度の間隔が1度であるため、基準点は次のプログラムで求まります。

meshcode = 5438

# Extract base latitude and longitude from the meshcode
south_latitude = int(meshcode[:2]) / 1.5
west_longitude = 100 + int(meshcode[2:4])

# Initial latitude and longitude intervals for the first order mesh
delta_latitude = 2/3  # 40 minutes in degrees
delta_longitude = 1  # 1 degree

上記で求めた基準点から残りの頂点を計算します。
それぞれの辺の長さの分だけ足すだけです。
4点ではなく5つの点を求めているのは、ポリゴンとして表現するためには最初と最後の点が結ばれていなければならないからです。

coordinates = [
    [
        [west_longitude, south_latitude],  # Bottom left
        [west_longitude + delta_longitude, south_latitude],  # Bottom right
        [west_longitude + delta_longitude, south_latitude + delta_latitude],  # Top right
        [west_longitude, south_latitude + delta_latitude],  # Top left
        [west_longitude, south_latitude]  # Close the loop
    ]
]

geojsonとして表現するには次のように記述します。
propertiesプロパティに関してはあってもなくてもよいです。

properties_data = {"meshcode": meshcode}

# Construct the GeoJSON
geojson = {
    "type": "Feature",
    "geometry": {
        "type": "Polygon",
        "coordinates": coordinates
    },
    "properties": properties_data
}

以上で1次メッシュコードからポリゴンを作成できました

第2次地域区画(2次メッシュ)の計算

続いて2次メッシュコードに取り組みます。
2次メッシュコードの該当区画は1次メッシュコードのそれとは求め方が異なり、1次メッシュコードの区画を縦横方向にそれぞれ8分割しています。

6桁の数字の左から4番目までの4桁が1次メッシュコード、左から5番目の数字が縦の等分の区画、左から6番目の数字が横の等分の区画を表しています。
また、緯度の間隔が5分、経度の間隔が7分30秒であるため、基準点は次のプログラムで求まります。

delta_latitude = 1/12 # 5 minutes in degrees
delta_longitude = 1/8 # 7.5 minutes in degrees
south_latitude += int(meshcode[4]) * delta_latitude
west_longitude += int(meshcode[5]) * delta_longitude

+=演算子による計算がありますが、元の変数は1次メッシュコードの計算で求められた緯度経度です。
これ以降、上位の区画で計算された緯度経度が元の変数に格納されます。
区画の個数の分だけ辺を加算することで格子番号の位置を求めています。
頂点の計算とgeojsonの作成は先ほどと同じですので割愛します。

基準地域メッシュ（第３次地域区画）の計算


2次メッシュコードの該当区画は1次メッシュコードの区画を縦横方向にそれぞれ8分割していましたが、3次メッシュコードの該当区画は2次メッシュコードの区画を縦横方向にそれぞれ10分割しています。
8桁の数字の左から6番目までの6桁が2次メッシュコード、左から7番目の数字が縦の等分の区画、左から8番目の数字が横の等分の区画を表しています。
少しずつ計算が複雑になってきましたね。
緯度の間隔は30秒、経度の間隔は45秒です。

delta_latitude = 1/120 # 30 seconds in degrees
delta_longitude = 3/240 # 45 seconds in degrees
south_latitude += int(meshcode[6]) * delta_latitude
west_longitude += int(meshcode[7]) * delta_longitude

２分の１地域メッシュ（分割地域メッシュ）の計算

4次メッシュコード以降はさらに計算が複雑になります。

9桁の数字の左から8番目までの8桁が3次メッシュコード、左から9番目の数字が基準地域メッシュの各辺を２等分して得られる４個の区画に、南西側、南東側、北西側、北東側の順に１～４の番号を付け、これをそれぞれの区画を示す数字を表しています。
縦横それぞれ2分割じゃダメだったんですかね。
緯度の間隔が15秒、経度の間隔が22.5秒であるため、基準点は次のプログラムで求まります。

delta_latitude = 15/3600 # 15 seconds in degrees
delta_longitude = 22.5/3600 # 22.5 seconds in degrees
south_latitude += [None, 0, 0, delta_latitude, delta_latitude][int(meshcode[8])]
west_longitude += [None, 0, delta_longitude, 0, delta_longitude][int(meshcode[8])]

左から9番目の数字に即した区域をif文を使わずに表現するためにリストを使用しています。
リストの先頭をNoneとしたのは要素を可読性向上のためです。

４分の１地域メッシュ（分割地域メッシュ）の計算


5次メッシュコードは4次メッシュコードの分割方法と同じです。

delta_latitude = 7.5/3600  # 7.5 seconds in degrees
delta_longitude = 11.25/3600  # 11.25 seconds in degrees
south_latitude += [None, 0, 0, delta_latitude, delta_latitude][int(meshcode[9])]
west_longitude += [None, 0, delta_longitude, 0, delta_longitude][int(meshcode[9])]

まとめ

地域メッシュコードからポリゴンを作成する方法について、5次メッシュコードまで説明しました。
ライブラリでは6次メッシュコードや5倍地域メッシュや2倍メッシュについても取り扱っているため、興味のある方はコードをご覧ください。

この記事で紹介したコードをjismeshcodeという名前で作成したことで、国勢調査のデータの可視化がとても楽になりました。
ライブラリ化するにあたり、参考リンクの記事にとてもお世話になりました。
また、作成したライブラリに貢献していただいた方々にも感謝いたします。
ライブラリ作成に関する話は別の記事に書こうと思います。

最後は下記リンクから入手できる2020年国勢調査における東京都近辺の人口及び世帯の4次メッシュのデータから作成した人口（総数）のヒートマップで締めくくりたいと思います。
https://www.e-stat.go.jp/gis/statmap-search?page=1&type=1&toukeiCode=00200521&toukeiYear=2020&aggregateUnit=H&serveyId=H002005112020&statsId=T001101&prefCode=13

こちらの可視化のコードもGitHubに載せています。
jismeshcodeを皆様のGISライフに役立てていただけますと嬉しいです。
ここまで読んでいただきありがとうございました🙌

参考リンク

https://zenn.dev/karaage0703/articles/db8c663640c68b

Discussion