はじめに
G空間情報センターから点群データをダウンロードする時、Webページに表示されるインデックスマップ上でメッシュをポチポチっとクリックしてダウンロードしたい場所を選択したりします。
そこで〇〇市のエリアのメッシュを全部選択してダウンロードするにはどうすれば良いか悩んだことはありませんか？私はありませんでした！昨日までは…
さて、のっぴきならない理由で特定の市町村の範囲と重なる点群データをダウンロードしなくてはならなくなった時にはどうすれば良いでしょうか？
点群データの分割範囲
点群データのファイル名は大抵 08OC7572.zip という形式になっています。これはデータの分割単位が国土基本図の図郭に基づいているため図郭コードがファイル名になっているものです。
08は平面直角座標の8系、OCはレベル50000、75はレベル5000、72はレベル500のコードに対応します。
国土基本図の図郭を生成
図郭のポリゴンおよびコードを生成するツールには以下のようなものがあります。
Pythonパッケージ
japan_basic_section (by のこのこさん）
平面直角座標系の番号を指定することで GeoDataFrameを生成することが出来ます。
from japan_basic_section.grid import Grid
from japan_basic_section.origin_coords import get_coord_info
system_number = 8
info = get_coord_info(system_number)
g = Grid(system_number, 50000)
grid_50000 = g.make_grid()
crs = info[1]["JGD2011"]
grid_50000 = grid_50000.set_crs(crs)
生成された GeoDataFrameの index に 図郭コード が入っています。
さて、一つの平面直角座標系の全体をグリッドに区切ると図郭数はいくつになるでしょうか。
| 地図情報レベル | 図郭数 | 範囲（東西x南北） |
|---|---|---|
| 50000 | 160 | 40km x 30km |
| 5000 | 16000 | 4km x 3km |
| 500 | 1600000 | 400m x 300m |
レベル500ではグリッドの数が多いので DataFrame の行数も多くなります。生成にも時間がかかるのでご注意。
QGISプラグイン
Japanese Grid Mesh（by MIERUNE)
こちらはプラグインなので QGISの GUIでパラメータを指定できます。「メッシュの作成範囲」を指定できるので必要最小限の図郭で済ませることが出来て有り難いです。また地域メッシュの作成機能も備えています。
行政区域ポリゴンデータの入手
行政区域のポリゴンデータは、例によって国土数値情報ダウンロードサイトから入手します。
静岡県点群データに合わせるため、静岡県の行政区域ポリゴンを入手します。Virtual Shizuoka の計測年度に合わせて、とりあえず2020年度の行政区域データ N03-20200101_22_GML.zip をダウンロード。
2024年以降だと浜松市の区が合併して変化しているので要注意。
file_path = Path("data/国土数値情報/行政区域/N03-20200101_22_GML/N03-20_22_200101.shp")
df_admin_area = gpd.read_file(file_path, encoding='shift_jis')
# JGD2011, 平面直角座標8系に設定
df_admin_area.to_crs(epsg=6676, inplace=True)
df_admin_area['N03_003'] = df_admin_area['N03_003'].fillna('')
属性としては以下のような情報が含まれています。郡名については None になっている行にfillna()を使って空文字列を入れています。
| 属性名 | 属性情報 |
|---|---|
| N03_001 | 都道府県名 |
| N03_002 | 北海道の振興局名 |
| N03_003 | 郡名 |
| N03_004 | 市区町村名 |
| N03_005 | 政令指定都市の行政区名 |
| N03_006 | 行政区コード |
plot()するとこのように描画されました。
ポリゴンの結合
シェープファイルを読み込んだ状態では、地物の数が1753個あります。同じ市町村が複数のポリゴンに分かれているためです。これらをdissolve()を使って集約、結合してみます。結合のキーとしては（郡名, 市区町村名）を指定します。
df_admin_area_dissolved = df_admin_area.dissolve(by=['N03_003', 'N03_004'])
結合後の地物の行数は43となります。
図郭ポリゴンの読み込み
今回はQGISプラグインで生成した図郭ポリゴンを利用します。
分かりやすくする為にレベル5000のデータを使っていますが、点群データのファイル名に対応させるにはレベル500の図郭コードが必要です。
file_path_grid_5000 = Path("data/国土基本図図郭/Grid_5000_shizuoka.gpkg")
df_grid_5000 = gpd.read_file(file_path_grid_5000)
QGISプラグインで生成したグリッドの数は2600個。静岡県全体を覆う図郭（下図の赤色部分）に限定すると 778個のグリッドがありました。
行政区域ポリゴンと図郭を重ねる
特定の市町村の行政区域ポリゴンに重なる部分がある図郭ポリゴンを抽出します。空間演算の交差(intersect)を計算します。
df_city = df_admin_area_dissolved.loc[('','磐田市')]
mask = df_grid_5000.geometry.intersects(df_city.geometry)
ax = df_admin_area_dissolved.plot(edgecolor='blue', alpha=0.4)
df_grid_5000[mask].plot(ax=ax, color='red', alpha=0.4)
intersectsによって、ポリゴン同士が交差している部分があればTrueが返ってきます。
Trueに対応する行を抽出すると、指定した市町村に対応するグリッドが抽出できました。
磐田市を覆う図郭はレベル5000の場合 29個となりました。
df_grid_5000[mask]['MapName'].to_list()
[
'08OC06','08OC14','08OC15','08OC16','08OC24',
'08OC25','08OC26','08OC34','08OC35','08OC44',
'08OC45','08OC54','08OC55','08OC56','08OC64',
'08OC65','08OC66','08OC73','08OC74','08OC75',
'08OC76','08OC83','08OC84','08OC85','08OC86',
'08OC93','08OC94','08OC95','08OC96']
To be continued...
レベル500の図郭ポリゴンを使って前述と同様に処理を行うと1514個の図郭が重なりました。
ではこの図郭コードを使って点群データをダウンロードしてみましょう（多分続く…）