はじめに

はじめまして！データサイエンティストの山上です。株式会社GEOTRA(@GEOTRA )にてGIS（Geographic Information System：地理情報システム）データを用いた分析を行っています。
本記事では指定された地域のメッシュコードを取得し、それをGeoDataFrameとして視覚化する方法を紹介します。

やること（背景）

指定したエリア（ポリゴン）に含まれるメッシュコードを取得する
- 指定された地域のメッシュを生成
- メッシュをGeoDataFrameとして可視化
前提：指定エリアのポリゴン情報をデータフレームで所持している
バージョン情報：
- Python version 3.9.5
- jismesh.utils version: 2.1.0
- Shapely version: 1.8.0
- Pandas version: 1.4.2
- GeoPandas version: 0.10.2
- Matplotlib version: 3.5.2
- NumPy version: 1.22.3

コードの概要

python

# ライブラリのインポート
import jismesh.utils as ju
from shapely import wkt
from shapely.ops import unary_union
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
import numpy as np

def make_mesh_list_semiauto(gdf: gpd.GeoDataFrame, tgt_level: int):
    # 例外処理
    if not isinstance(tgt_level, int) or not (1 <= tgt_level <= 6):
        raise ValueError("the level is unsupported.")
    
    # 最小値と最大値を設定し、最小値からstepずつ増加させて最大値を超えないリストを作成する
    def generate_list(min_value, max_value, step):
        # リストの生成
        values = list(np.arange(min_value, max_value, step))
        
        # 最大値がリストに含まれていない場合は追加
        if values[-1] < max_value:
            values.append(max_value)
            
        return values
    
    minx, miny, maxx, maxy = gdf.total_bounds # 対象のPOLYGONの最大最小緯度経度を取得
    tgt_polygon = unary_union(gdf.geometry.to_list()) # 対象のPOLYGONを結合
    
    level = 1
    # 候補となる1次メッシュコードを取得
    lon_list = generate_list(minx, maxx, 1)
    lat_list = generate_list(miny, maxy, 0.6667)
    candidate_meshcodes = list(set([ju.to_meshcode(lat, lon, level) for lat in lat_list for lon in lon_list]))
    
    # 対象のメッシュレベルまで交差判定を繰り返す
    if(tgt_level in [2,3,4,5,6]):
        while True:
            level += 1
            next_candidate_meshcodes = [meshcode for candidate_meshcode in candidate_meshcodes for meshcode in ju.to_intersects(candidate_meshcode, level)]
            lat_sw, lon_sw = ju.to_meshpoint(next_candidate_meshcodes, 0, 0)
            lat_ne, lon_ne = ju.to_meshpoint(next_candidate_meshcodes, 1, 1)
            polygon_list = [Polygon([(lon_sw[i], lat_sw[i]), (lon_sw[i], lat_ne[i]), (lon_ne[i], lat_ne[i]), (lon_ne[i], lat_sw[i])]) for i in range(len(next_candidate_meshcodes))]
            candidate_meshcodes = [next_candidate_meshcodes[i] for i in range(len(next_candidate_meshcodes)) if polygon_list[i].intersects(tgt_polygon)]
            
            if level == tgt_level:
                dict_meshcode2polygon = dict(zip(next_candidate_meshcodes, polygon_list))
                geometry = [dict_meshcode2polygon[meshcode] for meshcode in candidate_meshcodes]
                break
    else:
        if len(candidate_meshcodes) >= 2:
            lat_sw, lon_sw = ju.to_meshpoint(candidate_meshcodes, 0, 0)
            lat_ne, lon_ne = ju.to_meshpoint(candidate_meshcodes, 1, 1)
            polygon_list = [Polygon([(lon_sw[i], lat_sw[i]), (lon_sw[i], lat_ne[i]), (lon_ne[i], lat_ne[i]), (lon_ne[i], lat_sw[i])]) for i in range(len(candidate_meshcodes))]
        else:
            lat_sw, lon_sw = ju.to_meshpoint(candidate_meshcodes[0], 0, 0)
            lat_ne, lon_ne = ju.to_meshpoint(candidate_meshcodes[0], 1, 1)
            polygon_list = [Polygon([(lon_sw, lat_sw), (lon_sw, lat_ne), (lon_ne, lat_ne), (lon_ne, lat_sw)])]    
        
        lv1_candidate_meshcodes = candidate_meshcodes
        candidate_meshcodes = [lv1_candidate_meshcodes[i] for i in range(len(lv1_candidate_meshcodes)) if polygon_list[i].intersects(tgt_polygon)]
        dict_meshcode2polygon = dict(zip(lv1_candidate_meshcodes, polygon_list))
        geometry = [dict_meshcode2polygon[meshcode] for meshcode in candidate_meshcodes]
    
    result_gdf = gpd.GeoDataFrame(
        data = {
            'mesh_code': candidate_meshcodes,
        },
        geometry=geometry
    )
    
    # 結果の可視化
    fig,ax = plt.subplots()
    result_gdf.geometry.plot(ax=ax,alpha=.5)
    gdf.geometry.plot(ax=ax)
    plt.show()
    
    return result_gdf

参考：メッシュレベル

メッシュレベル
- 1次(標準地域メッシュ 80km四方): 1
- 2次(標準地域メッシュ 10km四方): 2
- 3次(標準地域メッシュ 1km四方): 3
- 4次(分割地域メッシュ 500m四方): 4
- 5次(分割地域メッシュ 250m四方): 5
- 6次(分割地域メッシュ 125m四方): 6

事例説明

群馬県を例にとって動かしてみる

指定エリア（ポリゴン）が小さい場合

↓ 125mメッシュの場合

6次メッシュ（標準地域メッシュ 125m四方）

python

areas = #ポリゴンのデータフレーム
areas.geometry = areas.geometry.apply(wkt.loads)
areas = gpd.GeoDataFrame(areas,geometry='geometry',crs='EPSG:4326')
gunnma = make_mesh_list_semiauto(areas, 6) # << 6次メッシュなので6に設定

指定エリア（ポリゴン）が巨大な場合

↓ 1kmメッシュの場合

３次メッシュ（標準地域メッシュ 1km四方）

python

areas = #ポリゴンのデータフレーム
areas.geometry = areas.geometry.apply(wkt.loads)
areas = gpd.GeoDataFrame(areas,geometry='geometry',crs='EPSG:4326')
gunnma = make_mesh_list_semiauto(areas, 3) # << 3次メッシュなので3に設定

参考文献

メッシュデータについて | 国土数値情報ダウンロードサイト

令和二年度　国勢調査　メッシュ定義書 | 政府統計の総合窓口

The Shapely User Manual - Shapely 2.0.4 documentation

Google Earthを用いたポリゴン作成 | Qiita

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

終わりに

最後までお読みくださりありがとうございました。

GEOTRAのQiitaでは、引き続き地理空間情報に関する情報を中心に技術系のテクニックを発信していきます。今後も皆さんのお役に立てるコンテンツを配信できればと思っておりますので、皆様のいいね・フォローをお待ちしています！

また、こちらの記事を読んでGEOTRAに興味を持ってくださった方は、是非Wantedlyから会社紹介記事もご覧ください！GEOTRAは、学生インターンを含め、現在一緒に働く仲間を募集しております。興味をお持ちの方は、Wantedly上の募集ページやストーリー、各種サイトをご確認ください！

企業HP：https://www.geotra.jp
note：https://note.com/2022geotra
Wantedly：https://www.wantedly.com/stories/geotra_introduction
LinkedIn：https://www.linkedin.com/company/geotra/

#Pythonプログラミング, #GeoDataFrame, #地理情報システム(GIS), #メッシュコード
#地図データ, #データ可視化, #ジオメトリ, #Shapely, #GeoPandas, #マッププロット