多方面からダウンロードした地物データをPythonでJSON（GeoJSON）に揃えてみた

地理データ処理の勉強中なのですが、ArcGISを持たず、QGISは機能が多すぎて自分のしたいことがどの機能でやれるのかサッパリ分からぬと日々ぼやいています

あくまで備忘録的に、試行錯誤してやれたことを公開して、数年後に振り返って、自分が何に興味を持っていたのか振り返っていきたいと思っています（初学者なので理論的な部分や基礎的な事項の欠落に関してはご容赦ください）

Python触ったことある人向け。

【目標】

①地理上の境界データを空間データとして扱える形にして、本来テーブルデータとして与えられている数値情報をオブジェクトに結び付けたい　（→データベースの構築）

②そのうえで、オブジェクトの包含関係、交差関係、面積や周長や空間距離に関する基本的な演算をできるようになりたい　（→データベースの幾何演算）

③その演算結果と、テーブルデータとの間で処理を行い、オブジェクトの可視化の際にヒートマップやカラーメッシュマップ、等高線マップ、点の分布、マップ上の文字およびグラフ表現、などの方法で表現したい　（→便利な描画フォーマットでの可視化）

【データの入手先】

・国土交通省国土政策局国土情報課GISホームページ：国土数値情報ダウンロードサービス

・e-Stat 政府の総合窓口：地図で見る統計（統計GIS)/境界データダウンロード

基本的に境界データが手に入るのは上の二つ。

・OpenStreetMap

・CSIS

こちらはまだ利用したことがないので、踏み込みません。

【座標系】

①世界測地系（緯度経度）と、②平面直角座標系（XY座標）と③その他（メルカトル図法など）があります。

それぞれに対し、EPSGコードという4桁くらいの数字コードが付与されています。

JGD座標系とEPSGの一覧表 | OpenなGISのこと

GISのための測地成果、測地系、楕円体、投影座標系、EPSGコードのまとめ - 自然環境保全のための周辺技術

e-Statでは①と②の両方が入手できます。

直観的に分かりやすいのは①ですが、GISの先生に聞いたところ、地域分析をする際には②にそろえた方がいいようです（ユニバーサル横メルカトルでもOKらしいが③なので手間が多い）。②は、全国の地方ごとに座標系が異なるので、必要な座標系にあわせて①を変換する必要があります。今回はまだ物理距離（km）に突っ込むわけではないのでとりあえず、全部①に合わせます。

【今回扱うデータ】

〈DID境界データ〉

▶国交省での配布（2015年度版の場合A16-15_01といった名前）

■拡張子：shp・dbf・shx , geojson , xml

■座標系：JGD2011（srsName="JGD2011 / (B, L)" frame="GC / JST"）

・ geojsonはデータベース形式としてかなり整っている。

"features"という名のリストの中に全オブジェクトが格納されおり、各オブジェクトのプロパティとジオメトリを取得することができる。

・xmlにおける緯度経度データは gml:Curve.gml:segments.gml:LineStringSegment.gml:posListに格納されている。gml:Curveそれぞれに対して緯度経度情報リストを取得すればよさそう。

▶e-Statでの配布（2015年度版の場合h27_did_01といった名前）

■拡張子：shp・dbf・shx , gml・xsd

■座標系：

→世界測地系緯度経度(ddmw)でDLした場合（srsName="EPSG:4612" srsDimention="2"）

→世界測地系平面直角座標系(xymw)でDLした場合（srsName="EPSG:24○○" srsDimention="2"）

e-StatではJGD2000を使用されている模様。国交省のほうはJGD2011なので、どちらかに合わせる方がいいようです。

〈小地域（国勢調査）境界データ〉

▶e-Statでの配布（2015年度版の場合h27ka01といった名前）

■配布拡張子：shp・dbf・shx , gml・xsd

■座標系：

→世界測地系緯度経度(ddmw)でDLした場合（srsName="EPSG:4612" srsDimention="2"）

→世界測地系平面直角座標系(xymw)でDLした場合（srsName="EPSG:24○○" srsDimention="2"）

こちらもJGD2000。

・gmlにおける緯度経度データはgml:featuremember.fme:h27ka○○○○○.gml:surfaceProperty.gml:Surface.gml:patches.　gml:PolygonPatch.gml:exterior.gml:LinearRing.gml:posListに格納されている。

▶国交省での配布

今のところ見つかりません

〈都市地域〉

▶国交省での配布（作成2005年、公表2006年が最新。A09-06_01といった名前）

配布拡張子：xmlのみ

座標系：JGD2000

国交省のGISホームページで配布されているデータには、JGD2000（特に緯度経度なのでEPSG:4612）で配布されているものとJGD2011（特に緯度経度なのでEPSG:6668）で配布されているものが混ざっているようです。

・こちらのxmlファイルはgmlという名前空間を利用していない。jps:GM_Curveやjps:GM_Surfaceというオブジェクトをとっています。

（例えば、DID境界データのxmlは、曲線オブジェクトにおけるポリラインの頂点を<gml:posList>の中に一括して、緯度経度の位置数値情報を格納しています。それに連なる形で、サーフェスオブジェクトの境界定義でこの曲線を参照する…というかたちをとっているのでとても見やすい。

一方で、この都市地域のxmlは、手間がひとつ多い定義の仕方をしています。

1. 曲線オブジェクトの始点・終点になる点を定義する

2. 曲線オブジェクトでこの点を始点・終点として参照しつつ、他の位置数値情報をポリラインの頂点点として直接定義して、その点の配列を<jps:GM_PointArray>の中に一括して格納することで曲線を定義

3. この曲線をサーフェスオブジェクトの境界定義の中で参照する

…といった様子。

全部の点オブジェクトを定義してから曲線オブジェクトを書き始め、全曲線オブジェクトを定義してからサーフェスオブジェクトを書き始める…というふうな記述方法が、とても見づらい…全体の構造としては見やすいけど、参照関係が見づらい…）

【Pythonで編集】

今回は、この3つの地物データを、一番読解性のすぐれる国交省配布DID人口データのGeoJSONのデータ形式に統一しようと思います。とりわけ、市街化区域の冗長な境界データを簡潔な形にしたい…

GMLはXMLベースの言語なので、全部XMLです。PythonでXMLを読み取りするxml.etree.ElementTreeモジュールを利用して、データの処理ができます。

〈都市地域〉→ GeoJSON

今回は北陸地域の某県のデータを用います。

メタデータを省けば、本体データの構造としては、オブジェクトの容器の中に、①点、②曲線、③有向曲線、④サーフェス、⑤サーフェスに付与される区分コード、の順に並んでいます。

ksj:OBJ	jps:GM_Point	id="n00001"
	jps:GM_Curve	id="c00001"
	jps:GM_OrientableCurve	id="c00001"
	jps:GM_Surface	id="a00001"
	ksj:BA01	id="BA01_00001"

①のGM_Pointと②のGM_Curveに全位置座標が入っているので、これをまず取得します。そして、位置の文字列データをfloat型に変えたうえで経度・緯度の順に格納し、jsonファイルの形式に変更させます。

GM_Curveの位置座標を規定する部分に関しては、間接指定と直接指定が混在している仕様になっているようなので、この部分だけは注意して抽出します。例えば下図のように、開始点と終着点だけがGM_Pointの緯度経度を参照しています。GM_PointArrayひとつごとに、リング状のサーフェス境界ひとつを表現しているようです。

jps:GM_PointArray	jps:GM_PointArray.column	jps:GM_Position.indirect	GM_PointRef.point idref = "n00001"
	jps:GM_PointArray.column	jps:GM_Position.direct	<DirectPosition.coordinate>緯度経度
	jps:GM_PointArray.column	jps:GM_Position.direct	<DirectPosition.coordinate>緯度経度
	…
	jps:GM_PointArray.column	jps:GM_Position.direct	<DirectPosition.coordinate>緯度経度
	jps:GM_PointArray.column	jps:GM_Position.indirect	GM_PointRef.point idref = "n00001"

以下使用したコードです。

import xml.etree.ElementTree as ET
import json

def main():
    tree = ET.parse('A09-06_16.xml') 
    root = tree.getroot()
    
    posLists = []
    pointList = []
    i=0
    for obj in root[1][0][0][1]:
        if("GM_Point" in obj.tag):
            posLists.append([])
            pointList.append(obj[0][0][0].text)
        if("GM_Curve" in obj.tag):
            for coordinate in obj[4][0][1][0]:
                if not(coordinate[0][0].text == None):
                    posLists[i].append(coordinate[0][0].text)
                else:
                    posLists[i].append(pointList[i])
            i += 1
    
    dic = {}
    dic["type"] = "FeatureCollection"
    dic["features"] = []  

    for posList in posLists:
        f = {}
        f["type"] = "feature"
        f["geometry"] = {"type":"Polygon"}
        f["geometry"]["coordinates"] = [[]]
        for postr in posList:
            dd = postr.split(" ")
            pos = [float(dd[1]),float(dd[0])] 
            #geojsonは経度→緯度の順番
            f["geometry"]["coordinates"][0].append(pos)
        dic["features"].append(f)
        
    with open('A09-06_16.json', 'w') as f:
        json.dump(dic, f, ensure_ascii=False)
    
if __name__ == '__main__':
    main()

※ただこれだけだと、領域を定義する輪状の曲線のオブジェクトしか存在しないため、2つの問題が出てきます。

①サーフェス定義に利用された外周用のGM_Curveと穴部分のGM_Curveの分類情報が抜け落ちている。

②都市地域/市街化区域/市街化調整区域/その他地域の分類情報が抜け落ちている。

実際にこのjsonファイルをそのまま利用する前に、上の2つの処理を必要に応じて追加する必要があります。これはフォーマットに合わせて、手動で加工したほうが早い気がしたので逐次加工しました。

GeoJSON フォーマット仕様

そして、この処理で作成されたJSONファイルはかなり自由にプロパティ情報を統合することができます。

ただし、以下の限界点に関しては、

・情報が古い。国交省の配布しているデータが2006年。東日本大震災前。

・このポリゴンデータで得られるのは、都道府県のマスタープランで定まる都市計画区域界と、区域区分界（線引き）、非線引き地域の用途地域区域の区分界。一方で、地域生活拠点や都心区域の指定、都市機能誘導区域や居住誘導区域といった立地適正化計画による詳細な区分の境界データは存在しない。

これについては、都道府県ページで得られる最新の都道府県マス、都市マス、立地適正化計画などを見て、修正点をアップデートしていくようにして、位置情報の細かいずれは諦める方針でいくしかなさそう…

〈小地域（国勢調査）境界データ〉→GeoJSON

最後に、こっちも同じ要領でやってしまいましょう。某県某市の小地域境界データ（GML形式）を、e-StatでDLします。他の二つに合わせるため、緯度経度の座標系を用います。

# -*- coding: utf-8 -*-

import xml.etree.ElementTree as ET
import json
import codecs

def main():
    tree = ET.parse('h27ka16201.xml') 
    root = tree.getroot()
    
    dic = {}
    dic["type"] = "FeatureCollection"
    dic["features"] = []

    
    for obj in root:
        if(obj.tag == "{http://www.opengis.net/gml}featureMember"):
            f={}
            f["type"] = "feature"
            f["properties"] = {}
            f["geometry"] = {"type":"Polygon"}
            f["geometry"]["coordinates"] = []
            
            idnum = int(obj[0].attrib["{http://www.opengis.net/gml}id"][2:])
            for prop in obj[0]:
                if(prop.tag == "{http://www.safe.com/gml/fme}KEY_CODE"):
                    f["properties"]["KEY_CODE"] = prop.text
                    f["properties"]["市町村名称"] = "△△市"
                    f["properties"]["小地域符号"] = idnum
            for prop in obj[0]:
                if(prop.tag == "{http://www.safe.com/gml/fme}JINKO"):
                    f["properties"]["人口"] = int(prop.text)
            for prop in obj[0]:
                if(prop.tag == "{http://www.safe.com/gml/fme}AREA"):
                    f["properties"]["面積"] = float(prop.text)
            f["properties"]["調査年度"] = 2015
            for prop in obj[0]:
                if(prop.tag == "{http://www.opengis.net/gml}surfaceProperty"):
                    for RingNum in range(len(prop[0][0][0])):
                        ddstr = str(prop[0][0][0][RingNum][0][0].text)
                        c = 0
                        ido=[]
                        keido=[]
                        for dd in ddstr.split(" "):
                            if(c%2 == 0):
                                ido.append(float(dd))
                            if(c%2 == 1):
                                keido.append(float(dd))
                            c +=1
                        f["geometry"]["coordinates"].append([])
                        for x in range(len(ido)):
                            f["geometry"]["coordinates"][RingNum].append([keido[x],ido[x]])
                            #RingNumが1を超える＝interiorBoundaryが存在する
            
            dic["features"].append(f)  
            

    file = codecs.open('h27ka16201.json', 'w','utf-8')
    json.dump(dic, file, ensure_ascii=False, indent=2)

 
if __name__=='__main__':
    main()

※こちらの方では、穴部分のデータの抜け落ちが生まれないように書きました。

コード普段書いてる人に見せたら唖然とされるくらいにはコードが汚いけれども、とりあえずこれで、JSONファイルにそろえることができたので、できます（？）

【〆】

目標①～③のうち、①のデータ変換を２つやってみました。

GeoJSONのいいところは、shpとdbfなら別々に扱っている地物の情報を、プロパティやジオメトリのカテゴリを含む地物オブジェクトでまとめて表現できるところだと思っています。

次以降の流れですが、このオブジェクトは一つのJSONの中に統合できるので、各種演算が容易にできそうです。これを次は、やってみるのと同時に、小地域データに小地域集計のCSVファイルから得られる年齢別人口などの情報を加えたいので、それをやってみます。

精進します！

追記：

結局、修論でGIS関係の知識を使うことになりました。

知識はどこで生きてくるのかわからないから、面白いことをやることには意味があるんですね。

Dossiri - Kamaeru

建築・都市・生活の領域に関する知識の体系化、技術に対する考察、書籍レビュー等