ArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみた

さて、本日はArcPyでオープンストリートマップのデータ（.osm)をShapeに変換してみようと思います。最近はまっているオープンストリートマップですが、出力できるデータ形式が.osm(実質xml)というものなので、なんだかなぁという感じです。とりあえずShapeで出力してArcMapでも読み込めるようにしようと思います。

仕様

JA:OSM XML - OpenStreetMap Wiki　に仕様が載っているのですが、大事なところだけをピックアップします。

データ形式

• ノード ブロック。属性として緯度経度を持つ

→各ノードの属性（名称とか）

<node id="1420826646" visible="true" version="1" changeset="9191805" timestamp="2011-09-02T10:28:00Z" user="watao" uid="393606" lat="35.5526160" lon="139.6462920">
  <tag k="amenity" v="police"/>
  <tag k="KSJ2:AdminArea" v="神奈川県横浜市港北区"/>
  <tag k="KSJ2:ADS" v="日吉本町1-23-21"/>
  <tag k="KSJ2:PubFacAdmin" v="都道府県"/>
  <tag k="name" v="港北警察署日吉駅前交番"/>
  <tag k="note" v="National-Land Numerical Information (Public Facility) 2006, MLIT Japan"/>
  <tag k="note:ja" v="国土数値情報（公共施設データ）平成19年　国土交通省"/>
  <tag k="source" v="KSJ2"/>
  <tag k="source_ref" v="http://nlftp.mlit.go.jp/ksj/jpgis/datalist/KsjTmplt-P02-v2_0.html"/>
 </node>

• ウェイ ブロック

→各ウェイごとに、そのノードへの参照（ノードへの参照から緯度経度を取得できるみたい）
→各ウェイの属性（名称とか）

<way id="378405415" visible="true" version="3" changeset="35164763" timestamp="2015-11-08T08:44:40Z" user="Madoka" uid="421690">
  <nd ref="3818193608"/>
  <nd ref="3818193607"/>
  <nd ref="3818193610"/>
  <nd ref="3824079880"/>
  <nd ref="3818193605"/>
  <nd ref="3818193606"/>
  <nd ref="3818193608"/>
  <tag k="amenity" v="cafe"/>
  <tag k="indoor" v="room"/>
  <tag k="level" v="0"/>
  <tag k="name" v="タリーズコーヒー　慶應日吉店"/>
  <tag k="opening_hours" v="Mo-Su 07:30-20:00"/>
  <tag k="room" v="yes"/>
  <tag k="source" v="survey"/>
 </way>

• リレーション ブロック

→各リレーションごとに、そのメンバーへの参照
→各リレーションの属性（名称とか）

<relation id="5647223" visible="true" version="3" changeset="36006131" timestamp="2015-12-17T10:28:22Z" user="nakajyou" uid="2435742">
  <member type="way" ref="377939570" role="part"/>
  <member type="way" ref="377939569" role="part"/>
  <member type="way" ref="159628078" role="part"/>
  <member type="way" ref="377939571" role="part"/>
  <member type="way" ref="377939573" role="part"/>
  <member type="way" ref="377939574" role="part"/>
  <member type="way" ref="377939572" role="part"/>
  <member type="way" ref="377939568" role="part"/>
  <member type="way" ref="377939567" role="part"/>
  <member type="relation" ref="5748542" role=""/>
  <member type="relation" ref="5647285" role=""/>
  <tag k="name" v="協生館"/>
  <tag k="type" v="building"/>
 </relation>

ノード＝ポイント、ウェイ＝ライン、ポリゴンですね。ウェイに関しては始点終点が一致している場合、ポリゴン、そうでない場合、ラインと判断できるっぽいです。リレーションというのは、各アイテムの関連性のことでしょうか（例えば、駅ビルと各テナントの関連性とか）？とりあえず、ここは不要と判断したので、本記事では割愛いたします。また、取得する属性は名称だけにしました（細かい仕様が載ってないので、カラムをどうするか判断がつきませんでした・・・）。

出力するShape

• ポイント
• ライン
• ポリゴン

どれも属性は名称のみです。ポイントのカテゴリ（例えば、スーパーマーケットとか）とかでShapeを分けた方がいいと思いますが、とりあえず今回はまとめて出力してみます。

出力するマップ

東急東横線の日吉駅前
↑とてもいいところです

サンプルコード

■exportosm.py

# -*- coding: utf-8 -*-
import arcpy
import os
from xml.etree import ElementTree

class OSM(object):

    def __init__(self,osm_path,out_dir):
        arcpy.env.workspace = out_dir
        self.osm = osm_path
        self.out_dir = out_dir
        self.file_name_point = "osm_point.shp"
        self.file_name_line = "osm_line.shp"
        self.file_name_polygon = "osm_polygon.shp"
        self.out_file_point = os.path.join(out_dir,self.file_name_point)
        self.out_file_line = os.path.join(out_dir,self.file_name_line)
        self.out_file_polygon = os.path.join(out_dir,self.file_name_polygon)
        self.spatial_reference = arcpy.SpatialReference(4612)
        self.node_dict = {}
        self.way_dict = {}

    def create_shape(self):
        """Shapeファイルを作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_point)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_point.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_point,
                                                "POINT",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_point) #カラム追加

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_line)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_line.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_line,
                                                "POLYLINE",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_line) #カラム追加

        exist_out_feature = arcpy.ListFeatureClasses(self.file_name_polygon)

        if len(exist_out_feature) == 0: #指定したディレクトリにosm_polygon.shpが存在しない場合は作成
            arcpy.CreateFeatureclass_management(self.out_dir,
                                                self.file_name_polygon,
                                                "POLYGON",spatial_reference = self.spatial_reference)

        self.add_columns(self.file_name_polygon) #カラム追加

    def add_columns(self,file_path):
        """Shapeファイルにカラムを追加します"""

        field_list = []
        fields = arcpy.ListFields(file_path)
        for field in fields:
            #if field.type != "Geometry":
            field_list.append(field.name)

        if not "name" in field_list: #nameカラムがなかったら作成
            arcpy.AddField_management(file_path, "name", "text")


    def read_osm(self):
        """osmファイルを読み込んで、必要な情報を取得します"""

        tree = ElementTree.parse(self.osm)

        for node in tree.findall(".//node"): #nodeを検索

            for e in list(node):
                if e.attrib["k"] == "name": #名称が存在した場合
                    self.node_dict[node.attrib["id"]] = {"name":e.attrib["v"],"lat":node.attrib["lat"],"lon":node.attrib["lon"]}
                    break

            else: #名称が無い場合はこっち
                self.node_dict[node.attrib["id"]] = {"name":"","lat":node.attrib["lat"],"lon":node.attrib["lon"]}

        tmp_way_dict = {}
        way_list = []
        name = ""

        for way in tree.findall(".//way"): #wayを検索

            for e in list(way):

                if e.attrib.has_key("ref"): #nodeidを取得

                    way_list.append(e.attrib["ref"])

                if e.attrib.has_key("k"):

                    if e.attrib["k"] == "name": #名称がある場合

                        name = e.attrib["v"]

                else:

                    name = ""

            t = set() #重複チェック(wayがlineかpolygonなのか)
            duplicate_list = [x for x in way_list if x in t or t.add(x)]

            if len(duplicate_list) > 0: #重複があったらポリゴンとみなす
                tmp_way_dict[way.attrib["id"]] = {"node_id" :way_list,"name":name,"type":"polygon"}
            else:
                tmp_way_dict[way.attrib["id"]] = {"node_id" :way_list,"name":name,"type":"line"}

            way_list = []

        node_coord = []
        for way in tmp_way_dict.items(): #way_dictに格納したnode_idから実際の座標を取得する

            for c in way[1]["node_id"]:
                node_coord.append([self.node_dict[str(c)]["lat"],self.node_dict[str(c)]["lon"]])

            self.way_dict[way[0]] = {"coord":node_coord,"name":way[1]["name"],"type":way[1]["type"]}

            node_coord = []



    def create_geometry(self):
        """メインメソッド"""

        self.create_shape() #Shape作成
        self.read_osm()     #OSM情報の取得
        self.create_node()  #nodeからポイントを作成
        self.create_way()   #wayからラインとポリゴンを作成

    def create_node(self):
        """nodeからポイントを作成します"""

        array = arcpy.Array()
        with arcpy.da.InsertCursor(self.out_file_point, ["SHAPE@","name"]) as cur:

            for nd in self.node_dict.items():

                point = arcpy.Point()
                point.X = nd[1]["lon"]
                point.Y = nd[1]["lat"]

                pointGeometry = arcpy.PointGeometry(point)

                cur.insertRow((pointGeometry,nd[1]["name"]))

    def create_way(self):
        """wayからラインとポリゴンを作成します"""

        #ラインを作成
        with arcpy.da.InsertCursor(self.out_file_line, ["SHAPE@","name"]) as cur:

            for wd in self.way_dict.items():

                if wd[1]["type"] == "line":

                    array = arcpy.Array()

                    for part in wd[1]["coord"]:

                        array.add(arcpy.Point(part[1],part[0]))

                    cur.insertRow((arcpy.Polyline(array,self.spatial_reference),wd[1]["name"]))


        features = []

        #ポリゴンを作成
        with arcpy.da.InsertCursor(self.out_file_polygon, ["SHAPE@","name"]) as cur:

            for wd in self.way_dict.items():

                if wd[1]["type"] == "polygon":

                    array = arcpy.Array()

                    for part in wd[1]["coord"]:

                        array.add(arcpy.Point(part[1],part[0]))

                    cur.insertRow((arcpy.Polygon(array,self.spatial_reference),wd[1]["name"]))

if __name__ == '__main__':
    osm = OSM("","") #osmのパスとShapeの出力先フォルダを指定
    osm.create_geometry()

結果を確認します。

完全に日吉駅前ですね！素晴らしい！！

位置もぴったりです！ただ、建物の四隅に必ずノードがあることに気づきました。これはオープンストリートマップの仕様なんですかね？！なんか気持ち悪い。あと、プログラム作った後に気づいたのですが、オープンストリートマップを使用する方はみなさんQGISとかを使っているのでしょうか？！（商用のGISには頼らないという思想？！）だとすると、このスクリプトはあまり役に立たないのかなぁと・・・一応、「read_osm」メソッドのところを流用すれば緯度経度や名称を取得することはできますが。。。

とりあえず今回は簡易的な出力を行ってみたのですが、osmファイルの仕様をきちんと理解すればもっと細かい精度でShapeに出すことができそうですね。

簡単ですが、今回は以上です！