geopandasとarcpyを使って基盤地図情報を編集して東京都の市区町村ごとに仕分けた後、建物の面積がどれくらいあるのかを集計してみる

主にgeopandas、そして一部でのみarcpyを使い、基盤地図情報で公開されているデータを加工して東京都の各市区町村ごとの建物面積がどれくらいあるのかを計算してみたいと思います。

手順には多くの無駄が含まれているかもしれませんが、記録としてそのまま残しています。

使用したGISソフト
QGIS…ダウンロードした基盤地図情報ファイルである建物ポリゴン(GML形式)をシェープファイルへ変換する
ArcGIS Pro…建物ポリゴンの自己交差の修復

コードの実行はjupyter notebookで行いました。

手順の概要
①建物を表すポリゴンデータを東京都全域分入手
②ポリゴンデータをシェープファイルに変換し、編集ができる状態にする
　②'ポリゴンデータの自己交差を修復
③位置に基づいたメッシュごとに分かれたシェープファイルを統合し、一つの大きなシェープファイルにする
④東京都の行政界ポリゴンデータを使い、ひとかたまりのシェープファイルにクリッピングを行い市区町村ごとに分ける
⑤建物ポリゴンの座標系を地理座標から平面直角座標系に変更し、面積の計算を行えるようにする(ここでのみarcpyを使いました)
⑥建物ポリゴンのシェープファイルに新しい列を追加し、各々のポリゴンの面積を計算した値を書き込む
⑦建物の面積をすべて足し合わせた値を市区町村ごとに記載したcsvファイルを作って出力する

1.建物を表すポリゴンデータを東京都全域分入手

まず、加工元のデータを基盤地図情報ダウンロードサービスからダウンロードしました。
東京都で建物が存在している全域をカバーしようとすると、「5339」のメッシュの範囲のうち35個ほどのデータをダウンロードする必要がありました。島しょ地域は今回対象としなかったため、ダウンロードする範囲には含んでいません。

2.ポリゴンデータをシェープファイルに変換し、編集ができる状態にする
・「基盤地図情報ビューア」での変換は断念
ここでダウンロードできるファイルは、GMLという専用の形式のようで、シェープファイルなどの一般的な形式に変換するためには「基盤地図情報ビューア」という専用のソフトを用いて閲覧・変換を行うというのが基本的な手法のようです。しかしこのソフトを使ってのシェープファイルへの変換には大変難儀しました。建物のポリゴンデータの容量が非常に大きいためか、このソフトを使って一つずつ変換を行おうとしても作業中に強制終了してしまいます。建物密度の低めな多摩地域はビューアからの変換が可能なものもありましたが、区部で同じことを行おうとするとほぼ失敗しました。
そのためこの方法はあきらめ、別の手段を試すことにしました。

・QGISを使ってGMLファイルを読み込み、シェープファイルとして出力
オープンソースのGISソフト「QGIS」を使えば、GMLファイルを直接読み込み、座標系(WGS84)を後から指定してやればシェープファイルとして問題なく出力することが可能だとわかりました。こちらであれば強制終了もなく、また比較的短時間で処理することができました。

基盤地図情報ダウンロードサイトから入手したファイルを解凍すると直接読み込めました

座標系をWGS84に手動で設定すると正しい位置に表示されます

このレイヤをエクスポートし、その際に形式として「シェープファイル」を選択。これをすべてのファイルに対して行い、建物ポリゴンデータのシェープファイルを作成できました。

2' 建物ポリゴンデータの自己交差を修復
変換して得た建物ポリゴンのシェープファイルには、自己交差と呼ばれるジオメトリ演算の際にエラーのもととなる構造上の問題が多く含まれた状態になってしまっていました。(どのタイミングで発生したのかは定かではないですが、シェープファイルとして出力した際に生じたのかもしれません)自己交差のあるファイルに対しては加工や演算ができないため、3の統合を行う前にこれらの問題を解消する必要がありました。すべてのシェープファイルをいったんarcGIS Proで開き、「ジオメトリの修復」を実行しました。

3.位置に基づいたメッシュごとに分かれたシェープファイルを統合し、一つの大きなシェープファイルにする
先ほどの画像の通り、現時点のシェープファイルは基盤地図情報のメッシュごとに分かれており、その範囲にまたがる複数の市区町村が入り混じっています。これを行政界に基づいて整理しなおすために、いったんすべてのシェープファイルを合体(マージ)させ、一つの大きな「東京都建物ポリゴン」を作成します。
この作業にはgeopandasを用いました。

import geopandas as gpd
import glob

#統合したいシェープファイルを格納したフォルダ内のすべてのファイルのパスをshpfileとして取得
shpfiles = [shpfile for shpfile in glob.glob('C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\*.shp')]

#すべてのシェープファイルをジオデータフレームとして読み込み、bld_mという一つのジオデータフレーム(東京都建物ポリゴン)に統合
bld_m = pd.concat([gpd.read_file(shpfile, encoding='utf-8') for shpfile in shpfiles])

#統合して出来上がった「東京都建物ポリゴン」をシェープファイルとして書き出し
bld_m.to_file(driver= "ESRI Shapefile", filename="C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\merge\\東京都建物ポリゴン.shp", encoding='utf-8')

#終了確認
print("done")

ちなみに、シェープファイルの読み込み・書き出しの際に「encoding='utf-8'」と文字コードの指定を入れ忘れてしまった場合、シェープファイルの属性テーブルのうち日本語の列が文字化けしてしまいました。
マージ作業に1時間ほど要して完了しました。

4.東京都の行政界ポリゴンデータを使い、ひとかたまりのシェープファイルにクリッピングを行い市区町村ごとに分ける
一つに統合した建物ポリゴンのシェープファイルに対して、個々の行政区ごとにクリッピングを行い、切り取った建物ポリゴンを個別のシェープファイルとして出力していきます。こちらでも引き続きgeopandasを使っています。
行政区のポリゴンは、全国市区町村界データをダウンロードしたものを使います。建物ポリゴンのシェープファイルと座標系が異なるため、ArcGISで座標変換を行ってWGS84に合わせています。
上記の行政界ポリゴンは全国の都道府県のデータを含んでいますので、前処理として東京都だけを取り出したものを作りました。

import geopandas as gpd

df = gpd.read_file("行政界.shp", encoding='shift-jis')

# 都道府県の情報を含む列[KEN]から、東京都のみを抽出
df = df[df['KEN'].str.contains('東京都')]

#文字コードはutf-8に変更し、新しく「東京都行政界」シェープファイルを出力
df.to_file('東京都行政界.shp', driver='ESRI Shapefile', encoding='utf-8')
#終了確認
print('done')

この東京都の行政界ポリゴンを使ってクリッピング作業を行いました。

import geopandas as gpd

#東京都建物ポリゴンの読み込み
mergedbld = gpd.read_file("C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\merge\\東京都建物ポリゴン.shp", encoding='utf-8')
#クリッピング用行政界ポリゴンの読み込み
tokyoborder = gpd.read_file("C:\\Users\\hoge\\デスクトップ\\行政界\\東京都行政界.shp", encoding='utf-8')

#行政界シェープファイルを読み込んだジオデータフレームから、市区町村の列「SIKUCHOSON」を取り出し、リストにする。
gyoseilist = tokyoborder["SIKUCHOSON"]

#クリッピング開始の目印
print("clippingstart")

#gyoseilistに格納された市区町村名に従って、東京都行政界ポリゴンからその名前を持つ行だけ取り出して東京都建物ポリゴンをクリップする。そしてその市区町村名で出力を繰り返す。
#出力先のパスはf_nameという変数を用意し、ファイル名のみ市区町村の名前で入れ替わるようにしておく。
for i in gyoseilist:
    f_name = "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\市区町村別ポリゴン\\" + i + ".shp"
    bld_clipped = gpd.clip(mergedbld, tokyoborder[tokyoborder["SIKUCHOSON"] == i], keep_geom_type=False)
#該当する市区町村がなかった場合の処理。クリッピングした建物ポリゴンが存在せず、空のジオデータフレームが出力された場合には空であるというメッセージのみを表示。島しょ地域は今回省いているので、空のデータが出てくるはずです。
    if len(bld_clipped.index) == 0:
        print(i + "is empty!")
#処理が終わった市区町村の建物ポリゴンデータは、その市区町村名を表示して完了の目印にする。
    else:
        bld_clipped.to_file(driver= "ESRI Shapefile", filename= f_name, encoding='utf-8')
        print(i)
#終了確認
print("done")

5.建物ポリゴンの座標系を地理座標から平面直角座標系に変更し、面積の計算を行えるようにする

先ほど出力された市区町村ごとの建物ポリゴンデータは、座標系がWGS84と地理座標のままなので、面積の計算には適していないはずです。面積を正確に算出できる平面直角座標系(JGD2011 9系)に変更しました。
シェープファイルの数が市区町村の数ぶんと大量にあるので、一つ一つGISソフト上に読み込んで変換を行うよりも、pythonの繰り返し処理で一気に変換したいと思いました。
pythonのライブラリを用いてシェープファイルの投影変換を行う方法は他にもあるようですが、ひとまずはarcpyで行うことが最もシンプルではないかと思いこれを使用することにしました。

下記のコードは、arcGIS Proのモデルビルダー機能を使って投影変換作業をpythonコードで出力したものを加工して作成しています。

import arcpy
#出力しなかった島しょ地域の市区町村名は省いたリストを用意しておき、今度はこちらを繰り返し処理のためのリストとして使用します。
gyoseilist = ["千代田区","中央区","港区","新宿区","文京区","台東区","墨田区","江東区","品川区","目黒区","大田区","世田谷区","渋谷区","中野区","杉並区","豊島区","北区","荒川区","板橋区","練馬区","足立区","葛飾区","江戸川区","八王子市","立川市","武蔵野市","三鷹市","青梅市","府中市","昭島市","調布市","町田市","小金井市","小平市","日野市","東村山市","国分寺市","国立市","福生市","狛江市","東大和市","清瀬市","東久留米市","武蔵村山市","多摩市","稲城市","羽村市",
"あきる野市","西東京市","瑞穂町","日の出町","檜原村","奥多摩町"]

for i in gyoseilist:
    print(i)
    input_data = "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\市区町村別ポリゴン\\" + i + ".shp"
    output_data = "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\prj\\" + i + "_prj.shp"
    arcpy.Project_management(in_dataset=input_data, out_dataset=output_data, out_coor_system="PROJCS['JGD_2011_Japan_Zone_9',GEOGCS['GCS_JGD_2011',DATUM['D_JGD_2011',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',139.8333333333333],PARAMETER['Scale_Factor',0.9999],PARAMETER['Latitude_Of_Origin',36.0],UNIT['Meter',1.0]]",
    transform_method=["JGD_2011_To_WGS_1984"], in_coor_system="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", preserve_shape="NO_PRESERVE_SHAPE", max_deviation="", vertical="NO_VERTICAL")
#終了確認
print("done")

6.建物ポリゴンのシェープファイルに新しい列を追加し、各々のポリゴンの面積を計算した値を書き込む

5の作業で平面直角座標系に変換した市区町村別の建物シェープファイルに、一つ一つの建物ポリゴンの面積を表す列を追加しました。

import geopandas as gpd
import pandas as pd
import glob
import os

#平面直角座標系に変換したシェープファイルをフォルダから取り出してリストにする
shpfiles = glob.glob('C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\prj\\*.shp')
#面積を追加した新しいシェープファイルの保存先
newpath = "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\area\\"

#取り出したリストから順番に、[area_m2]という列を新しく追加して面積を入れ、新しいシェープファイルを出力する。
for shpfile in shpfiles:
    geom = gpd.read_file(shpfile)
    geom.insert(10,"area_m2", geom.area)
#出力するファイル名を新しく作成するためにパス名を取得し、拡張子を除いたファイル名のみを抽出する。
    filename = (os.path.basename(shpfile).split('.', 1)[0]) 
    geom.to_file(driver= "ESRI Shapefile", filename= newpath + filename + "_area.shp", encoding='utf-8')
#出力が終わった目印のために名前を表示
    print(filename)

#終了確認
print("done")

7.建物の面積をすべて足し合わせた値を市区町村ごとに記載したcsvファイルを作って出力する

個々の建物ポリゴンの面積を列に加えたので、最後に各市区町村ごとに合計値をまとめたレポートを作ってcsvファイルとして出力しました。これがゴールに設定していた「市区町村ごとの建物の面積」になります。

import geopandas as gpd
import pandas as pd
import os

#手順5の投影変換で使ったものと同じ市区町村名リストを使いました。
gyoseilist = ["千代田区","中央区","港区","新宿区","文京区","台東区","墨田区","江東区","品川区","目黒区","大田区","世田谷区","渋谷区","中野区","杉並区","豊島区","北区","荒川区","板橋区","練馬区","足立区","葛飾区","江戸川区","八王子市","立川市","武蔵野市","三鷹市","青梅市","府中市","昭島市","調布市","町田市","小金井市","小平市","日野市","東村山市","国分寺市","国立市","福生市","狛江市","東大和市","清瀬市","東久留米市","武蔵村山市","多摩市","稲城市","羽村市","あきる野市","西東京市","瑞穂町","日の出町","檜原村","奥多摩町"]

#手順6で作成したファイルを保存したフォルダのパス
shp_path = "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\area\\"
#面積の集計結果の出力先
report_path= "C:\\Users\\hoge\\デスクトップ\\建物シェープ作業フォルダ\\集計\\"

#面積の集計結果を入力するための空のリストを作成
areasum_list = []

for i in gyoseilist:
#市区町村名ごとにシェープファイルを読み込んでジオデータフレームを作成
    filename = shp_path + i + "_prj_area.shp"
    shpfile = gpd.read_file(filename, encoding='utf-8')
#面積の集計
    sum = shpfile["area_m2"].sum()
#先ほどの空のリストに面積の集計結果を入れる
    areasum_list.append(sum)

#面積の集計結果が集まったことの目印にリストを表示
print(areasum_list)

#カラム[市区町村]にgyoseilist,カラム[建物面積]にareasum_listを入れた表を作成
df = pd.DataFrame({ "市区町村" : gyoseilist,"建物面積" : areasum_list})

#csvにして出力
df.to_csv(report_path + "面積集計.csv")

#終了確認
ptint("done")

結果の確認
出力されたレポートを確認しました。
問題なく出力できたようです。建物面積の合計値順に並べてみると思っていた順位と違いましたが、単純に各市区町村の面積に沿った結果になっているだと思います。
行政界ポリゴンを使って市区町村全体の面積を出しておき、全体に占める建物面積の割合も出しておけばよかったな、と感じました。