GIS奮闘記

元GISエンジニアの技術紹介ブログ。主にPythonを使用。

スポンサーリンク

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を分けた方がいいと思いますが、とりあえず今回はまとめて出力してみます。

出力するマップ

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

f:id:sanvarie:20160418082715p:plain

サンプルコード

■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()

結果を確認します。

f:id:sanvarie:20160418082914p:plain

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

f:id:sanvarie:20160418083112p:plain

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

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

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

OpenStreetMapを使ってみよう!

さて本日はマニアックなArcPyのプログラムからは離れてOpenStreetMapを紹介しようと思います。実は私も最近知りました。GISに携わる者としてこんな面白いものを知らなかったのは実に恥ずかしい。。。

OpenStreetMapとは

OpenStreetMap(OSM)は、道路地図などの地理情報データを誰でも利用できるよう、フリーの地理情報データを作成することを目的としたプロジェクトです。
誰でも自由に参加して、誰でも自由に編集でき、誰でも自由に利用する事が出来ます。」

要するにウィキペディアの地図版ってとこですかね?!以下サイトから利用することができます。

https://openstreetmap.jp/

地図の見方

以下ボタンをクリックしてください。

f:id:sanvarie:20160414155913p:plain

 

地図が表示されます。なかなか綺麗ですね!

f:id:sanvarie:20160414160032p:plain

 

拡大するとこんな感じです。ここが拡大の限界なのが微妙ですが、なかなかいいですね!

f:id:sanvarie:20160414160145p:plain

 

地図の編集の仕方

以下ボタンをクリックしてください。

f:id:sanvarie:20160414160355p:plain

 

編集したい場所に移動します。

f:id:sanvarie:20160414160532p:plain

 

ポイント、ライン、エリアのいずれかを選択します。

f:id:sanvarie:20160414160624p:plain

 

ポイントを置きます。

f:id:sanvarie:20160414160741p:plain

 

必要な情報を書き込みます。

f:id:sanvarie:20160414160859p:plain

 

保存します。

f:id:sanvarie:20160414160959p:plain

 

Saveして完了です。

f:id:sanvarie:20160414161058p:plain

 

これで自分の作成した情報がOpenStreetMapに反映されます。作図だけではなく、既存の情報を修正することもできます(コンビニが閉店してラーメン屋になった場合、ポイントの情報を変更する)。

地図データの出力

f:id:sanvarie:20160414161908p:plain

 

エクスポートを押すとデータをエクスポートできます。ただ、形式が「.osm」という形式のため、ArcGISなどで読み込むことができません(なんかプラグイン的なのが公開されていましたが、ArcGIS10.2.0用のはありませんでした・・・)。こいつをShapeに変換できたらすごく便利なのですが、それは今後調べてみようと思います。

 

地図を閲覧するだけならGoogleMapの方がいいと思いますが、この取組自体がとても素晴らしいと感じますし、色々メリットも多いと思います(例えば、通常地図には載らないであろう情報を自分や自分の仲間むけに載せることができますね。私が勤める会社の業務的にも使えそうな気がします)。

 

サイトの使い方などを少し勉強する必要がありますが、ぜひ皆さんも活用してみてはいかがでしょうか!

ArcPyで壊れたフィーチャクラスを抽出する方法

さて、本日は「ArcPyで壊れたフィーチャクラスを抽出する方法」です。フィーチャクラスが壊れるって何?!と思う方もいらっしゃるかと思いますが、フィーチャクラスは壊れるんです。本当に貧弱です。

確認方法

属性テーブルのヘッダカラムをダブルクリックしてください。

例えば、こういったフィーチャクラスがある場合
f:id:sanvarie:20160401155648p:plain

フィーチャの件数は861件です。フィーチャクラスが壊れている状態でヘッダカラムをダブルクリックすると
f:id:sanvarie:20160401155850p:plain

フィーチャの件数が737件に減ってしまいました。

もうこうなってしまったら最後、新しくフィーチャクラスを作成するしかありません。また、確認方法ですが、上記の方法でも問題ないのですが、フィーチャクラスが数多くある場合、あまりにも面倒です。なので、Pythonでの抽出方法を考えてみました。

確認方法

arcpy.da.FeatureClassToNumPyArray を使用します。これはFeatureClassをNumPyArrayに変換してくれる便利な関数です。

①この関数の第二引数にOBJECTIDのみを渡す
→正常なフィーチャ数を取得(861件)

②この関数の第二引数にOBJECTID以外も渡す
→異常なフィーチャ数を取得(737件)

この仕組みを利用すれば、壊れたフィーチャクラスを抽出できますね。

環境

Windows7 64bit
ArcGIS10.2.0
Python2.7.3

サンプルソース

import arcpy

arcpy.env.workspace = "" #GDBが存在する一個上のフォルダを指定
broken_list = []

for w in arcpy.ListWorkspaces(workspace_type="FileGDB"):

    arcpy.env.workspace = w

    for f in arcpy.ListFeatureClasses():

        array1 = arcpy.da.FeatureClassToNumPyArray(f,["OID@"],null_value=-9999)
        array2 = arcpy.da.FeatureClassToNumPyArray(f,["OID@","SHAPE@XY"],null_value=-9999)

        if len(array1) != len(array2):
            #GDB、フィーチャクラス、正常な件数、異常な件数を格納
            broken_list.append([w,f,len(array1),len(array2)])

print broken_list

結果はこのようになります。壊れたフィーチャクラスの結果のみ抽出することができました。
f:id:sanvarie:20160401161336p:plain

そもそもこんなしょっちゅうフィーチャクラスが壊れるってどういうことですかESRIさん!と言いたいところですが、言ってもどうしようもありません。とりあえず、運悪く壊れてしまった場合はこのように対象フィーチャクラスを抽出すればいいですね。

簡単ですが、本日は以上です。

【PythonとGIS】ArcPyで頂点の間引き処理を進化させてみた

さて、本日は頂点の間引き処理を行ってみたいと思います。ArcMapの頂点の間引きというジオプロセシングツールを使えばそれを実現することが可能です。頂点の間引き (Generalize)—ヘルプ | ArcGIS for Desktop
ただ、以下のような問題もあります。

  • Standardライセンスからのみ使用可能で、Basicまでしか持っていないという場合は使用不可

今回はBasicライセンスしか持っていないという場合でもそれに対応できるPythonスクリプトを作成してみました。

ソースはGistにも公開しています。https://gist.github.com/sanvarie/51bcdd9e1668081bde75

環境

windows7 64bit
ArcGIS10.2.0
Python2.7.3

事前準備

  1. GDBを作成する(今回は「Generalize.gdb」を作成しました)
  2. フィーチャクラスを作成する(今回は「Line」というラインのフィーチャクラスを作成しました)

処理概要

  1. フィーチャクラスを作成する(間引いたラインを作成するフィーチャクラス)
  2. 始点終点は間引き処理対象外とする
  3. 頂点間の角度と距離を計測して間引き対象/対象外を選別する
  4. 頂点間の角度の差異が1度以下、かつ、頂点間の長さが1m以上を間引き対象とする

対象データ

こういったデータを作成しました。
f:id:sanvarie:20160325092748p:plain

各頂点はこのようになっています。
f:id:sanvarie:20160325093200p:plain

サンプルソース

■generalize.py

# -*- coding: utf-8 -*-
import arcpy
import math

class Generalize(object):

    def __init__(self,workspace,in_feature,out_feature):
        arcpy.env.workspace = workspace
        self.workspace = workspace
        self.in_feature = in_feature
        self.out_feature = out_feature
        self.generalize_list = list()
        self.check_feature_class()

    def check_feature_class(self):
        """out_featureに指定したフィーチャクラスが無い場合、作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.out_feature)

        if len(exist_out_feature) == 0:
            try:
                print u"フィーチャクラスを作成します"
                arcpy.CreateFeatureclass_management(self.workspace, self.out_feature, "POLYLINE")
            except Exception:
                u"フィーチャクラスの作成に失敗しました"

    def set_generalize_list(self):
        """角度と距離から削除対象の頂点か判定してリストに格納する

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離

        """

        for i,pt in enumerate(self.generalize_list):

            #削除対象のポイントの場合
            if pt[5] == 0:

                #角度計算用に対象の頂点とその前後の頂点の角度を取得する
                angle_first = self.generalize_list[i-1][7]
                angle_next = pt[7]
                angle_last = self.generalize_list[i+1][7]

                #対象の頂点とその前後のポイントの角度の差が1以下、かつ、対象の頂点間の長さが1以上の場合
                if math.fabs(math.fabs(angle_first) - math.fabs(angle_next)) <= 1 \
                and math.fabs(math.fabs(angle_next) - math.fabs(angle_last)) <= 1 and pt[8] >= 1:

                    self.generalize_list[i].append(1)

                else:

                    self.generalize_list[i].append(0)
            else:

                self.generalize_list[i].append(0)

    def calc_angle_length(self):
        """各頂点の角度や長さなどを計算する関数し、計算結果をリストに格納する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離
                         [9]:削除対象区分2  0:削除対象外,1:削除 ←[5]はこれを抽出する前段階のもの
        """

        x2,y2,pg = float(),float(),object()
        generalize_list2 = list()

        for i,pt in enumerate(self.generalize_list):

            #始点、終点でもない場合
            if pt[5] == 0:

                #角度を計算
                result = math.atan2(pt[2]-y2,pt[1]-x2)

                #角度をセット
                self.generalize_list[i].append(math.degrees(result))

                #距離をセット
                self.generalize_list[i].append(pt[6].distanceTo(pg))

            #始点の場合
            elif pt[5] == 1:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i+1][2],pt[1]-self.generalize_list[i+1][1])

                #角度をセット
                if math.degrees(result) < 0:
                    self.generalize_list[i].append(math.degrees(result) + 180)
                else:
                    self.generalize_list[i].append(math.degrees(result) - 180)

                #距離は0でOK
                self.generalize_list[i].append(0)

            #終点の場合
            elif pt[5] == 2:
                #角度を計算
                result = math.atan2(pt[2]-generalize_list2[i-1][2],pt[1]-generalize_list2[i-1][1])

                #角度をセット
                if math.degrees(result) < 0:
                    self.generalize_list[i].append(math.degrees(result))
                else:
                    self.generalize_list[i].append(math.degrees(result))

                #距離は0でOK
                self.generalize_list[i].append(0)

            generalize_list2.append(self.generalize_list[i])
            x2,y2,pg = pt[1],pt[2],pt[6]

        #角度と距離から削除対象の頂点か判定
        self.set_generalize_list()

    def get_point_coord(self):
        """ラインの頂点ごとの情報をリストに格納する"""

        lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]

        #頂点ごとにOID、X、Y、始点、終点の情報を格納する
        for line in lines:
            for parts in line[1]:
                for part in parts:
                    self.generalize_list.append([line[0],part.X,part.Y,line[1].firstPoint.X,line[1].lastPoint.X])

    def calc_first_last(self):
        """各頂点が始点なのか終点なのかそれ以外なのか(削除対象)を判定する関数 """

        for i,pt in enumerate(self.generalize_list):

            #対象のポイントのXが始点の場合
            if pt[1] == pt[3]:

                #始点で削除非対象
                self.generalize_list[i].append(1)

            elif pt[1] == pt[4]:

                #終点で削除非対象
                self.generalize_list[i].append(2)

            else:

                #削除対象
                self.generalize_list[i].append(0)

    def calc_point_distance(self):
        """各種関数を呼び出し、self.generalize_listに要素を追加していく"""

        #ラインの頂点ごとの情報を取得
        self.get_point_coord()

        #各頂点の始点終点それ以外を判定
        self.calc_first_last()

        #頂点のジオメトリ情報を取得
        self.create_point_geometry()

        #頂点の角度と長さの計算
        self.calc_angle_length()

        #削除対象の頂点を間引いてラインを作図する
        self.generalize_line()

    def generalize_line(self):
        """頂点を間引いてラインを作図する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:始点X
                         [4]:終点X
                         [5]:削除対象区分 0:削除対象,1:始点で削除非対称,2:終点で削除非対称
                         [6]:ポイントジオメトリ
                         [7]:角度
                         [8]:距離
                         [9]:削除対象区分2  0:削除対象外,1:削除 ←[5]はこれを抽出する前段階のもの
        """

        oid = int()
        array = arcpy.Array()

        with arcpy.da.InsertCursor(self.out_feature, ["OID@", "SHAPE@"]) as cur:

            for i,pt in enumerate(self.generalize_list):

                #OBJECTIDが次に行ったタイミングでラインを作図する。
                if pt[0] != oid and i != 0:
                    cur.insertRow((oid,arcpy.Polyline(array)))
                    array = arcpy.Array()
                    array.add(arcpy.Point(pt[1], pt[2]))
                else:
                    #削除対象の頂点ではない場合のみ頂点の情報を格納していく
                    if pt[9] != 1:
                        array.add(arcpy.Point(pt[1], pt[2]))

                oid = pt[0]

            cur.insertRow((oid,arcpy.Polyline(array)))

    def create_point_geometry(self):
        """ポイントジオメトリを作成してリストに格納"""
        point = arcpy.Point()

        for i,pt in enumerate(self.generalize_list):
            point.X = pt[1]
            point.Y = pt[2]
            self.generalize_list[i].append(arcpy.PointGeometry(point))

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\Generalize.gdb"

    #間引き元のフィーチャクラス名
    in_feature = "Line"

    #間引き先のフィーチャクラス名
    out_feature = "LineGeneralize"

    #インスタンス化
    generalize = Generalize(workspace,in_feature,out_feature)

    #頂点を間引く
    generalize.calc_point_distance()

結果を確認すると問題なく、頂点が間引きされていることがわかりました。
f:id:sanvarie:20160325100013p:plain

課題

  1. データ量が多い場合は少し時間がかかってしまう
  2. 属性が存在する場合は、間引いた先のフィーチャクラスにはその属性は引き継がれない(generalize_lineメソッドを修正すればいける)

課題は少しありますが、角度や距離をうまく使用すれば、ArcMapの機能よりも柔軟な使用ができそうですね。簡単ではありますが、本日は以上です。

【PythonとGIS】ArcPyでラインの端点の角度を取得する

さて、本日はラインの端点の角度を取得しようと思います。ラインの端点にはシンボルを置くことが多いので、端点の角度を取得して、そこにのっかるポイントに角度を与える場合などに使えるかと思います。

環境

Windows7 64bit
ArcGIS10.2
Python2.7.3

事前準備

  1. 対象GDBの作成(今回はGeneralize.gdbを作成しました)
  2. フィーチャクラスの作成(今回はLineを作成しました)

プログラムで取得できる値

対象のフィーチャクラスにおいて、端点の角度を含む以下情報をListで返します。端点以外の頂点の情報は格納しません。

[0]:OBJECTID
[1]:X
[2]:Y
[3]:頂点区分 1:始点 2:終点
[4]:角度

対象データ

以下のようなサンプルデータを作成しました。
f:id:sanvarie:20160322151530p:plain

左から順に作図しました。
f:id:sanvarie:20160322151650p:plain

サンプル

■calcangle.py

# -*- coding: utf-8 -*-
import arcpy
import math

class CalcAngle(object):

    def __init__(self,workspace,in_feature):
        arcpy.env.workspace = workspace
        self.in_feature = in_feature
        self.generalize_list = list()

    def calc_angle(self):
        """各頂点の角度を計算する関数

        self.generalize_list --ラインの始点終点などの情報
                         [0]:OBJECTID
                         [1]:X
                         [2]:Y
                         [3]:頂点区分 1:始点,2:終点
                         [4]:角度
        """

        for i,pt in enumerate(self.generalize_list):

            #始点の場合
            if pt[3] == 1:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i+1][2],pt[1]-self.generalize_list[i+1][1])

                #角度をセット
                if math.degrees(result) < 0:

                    self.generalize_list[i].append(abs(math.degrees(result)))

                else:

                    self.generalize_list[i].append(abs(math.degrees(result) - 360))

            #終点の場合
            elif pt[3] == 2:

                #角度を計算
                result = math.atan2(pt[2]-self.generalize_list[i-1][2],pt[1]-self.generalize_list[i-1][1])

                #角度をセット
                if math.degrees(result) < 0:

                    self.generalize_list[i].append(abs(math.degrees(result)))

                else:

                    self.generalize_list[i].append(abs(math.degrees(result) -360))

            else:
                self.generalize_list[i].append(0)

    def get_point_coord(self):
        """ラインの頂点ごとの情報をリストに格納する"""

        lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]

        #頂点ごとにOID、X、Y、始点、終点の情報を格納する
        for line in lines:
            for parts in line[1]:
                for part in parts:
                    #始点の場合
                    if part.X == line[1].firstPoint.X and part.Y == line[1].firstPoint.Y:
                        self.generalize_list.append([line[0],part.X,part.Y,1])
                    #終点の場合
                    elif part.X == line[1].lastPoint.X and part.Y == line[1].lastPoint.Y:
                        self.generalize_list.append([line[0],part.X,part.Y,2])

    def calc_point_distance(self):
        """各種関数を呼び出し、self.generalize_listに要素を追加していく"""

        #始点終点の情報を取得
        self.get_point_coord()

        #頂点の角度の計算
        self.calc_angle()

        return self.generalize_list

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\Generalize.gdb"

    #対象のフィーチャクラス名
    in_feature = "Line"

    #インスタンス化
    generalize = CalcAngle(workspace,in_feature)

    #角度計算
    angle_list = generalize.calc_point_distance()

    print angle_list

実行結果はこのようになりました。
f:id:sanvarie:20160322152039p:plain

OBJECTIDが1のデータで確認すると始点(以下画像の緑丸)の角度=156.1363677663898、終点(以下画像の赤丸)の角度=336.1363677663898 で問題なさそうですね。
f:id:sanvarie:20160322152247p:plain

これを使用すれば端点の角度を利用してそこに接するポイントに角度を付与することができるかと思います。また、これを応用して、各頂点ごとの角度などを取得してもおもしろそうですね。

簡単ではありますが、本日は以上です。

【PythonとGIS】ArcPyでライン接合を行ってみた

さて、本日はArcPyでライン接合を行ってみようと思います。背景などでけっこうラインがぶちぶち切れてしまっていることが多いので、この機能があればラインのフィーチャ数を減らすことができますね。Advancedライセンスがあれば、ArcMapで実行できるのですが(うまくいくかは知らない)、いかんせんAdvancedはお高いですよね(けっこういい車が買えるお値段です)。そんなものを購入できるはずもなく、「使えないなら作ろう!」と思い立った次第であります。ArcMapでのライン接合についてはこのサイトに記載してあります。ラインの接合 (Unsplit Line)—データ管理ツールボックス | ArcGIS Desktop

今回もソースをGistに公開しました。https://gist.github.com/sanvarie/410ea0c1bab366f85a68

環境

Windows7 64bit
ArcGIS10.2
Python2.7.3

事前準備

  1. ファイルジオデータベースの作成(今回はArcPyAdvanced.gdbという名前にしました)
  2. 接合対象のフィーチャクラスの作成(今回はLineSplitという名前にしました)
  3. 接合先のフィーチャクラスの作成(今回はLineUnsplitという名前にしました)←接合先のフィーチャクラスが存在しなかったら作成するようにしたので、作らなくても大丈夫です。
  4. Pandasのインストール

接合対象/対象外データ

以下参照
f:id:sanvarie:20160319170244p:plain

サンプルソース

ライン接合処理のサンプルです。以下に処理の概要を記載します。

  1. 接合対象ラインを抽出(つながりのあるラインを抽出)
  2. 接合対象ラインの中でループしているラインを判定
  3. 接合対象ラインの中で端のラインとそのラインの端点(他のラインと接していない端点)を抽出(ループしているラインは除く)
  4. 接合対象ラインの端点から作図順を決める
  5. 接合対象ラインを作図
  6. 接合対象外ラインを作図

■unsplitline.py

# -*- coding: utf-8 -*-
import arcpy
from collections import Counter
import os
import traceback
import pandas as pd

class UnsplitLine():

    column_unsplit_length = "UNSPLIT_LENGTH" #作成するカラム名

    def __init__(self,workspace,in_feature,out_feature):
        arcpy.env.workspace = workspace
        self.workspace = workspace
        self.in_feature = in_feature
        self.out_feature = out_feature
        self.line_list = []
        self.match_point_list = []
        self.finish_line_list = []
        self.divide_line_list = []
        self.list_draw_point_sort = []
        self.finish_oid = []
        self.count = UnsplitLine.counter(self)
        self.check_feature_class()
        self.check_columns()

    def check_feature_class(self):
        """out_featureに指定したフィーチャクラスが無い場合、作成します"""

        exist_out_feature = arcpy.ListFeatureClasses(self.out_feature)

        if len(exist_out_feature) == 0:
            print u"フィーチャクラスを作成します"
            arcpy.CreateFeatureclass_management(self.workspace, self.out_feature, "POLYLINE")

    def check_columns(self):
        """out_featureに指定したフィーチャクラスに「UNSPLIT_LENGTH」カラムが存在しない場合、カラムを作成します"""
        field_list = []
        fields = arcpy.ListFields(self.out_feature)
        for field in fields:
            if field.type != "Geometry":
                field_list.append(field.name)

        if not UnsplitLine.column_unsplit_length in field_list:
            print u"カラムを作成します"
            arcpy.AddField_management(self.out_feature, UnsplitLine.column_unsplit_length, "Double")

    def export_upsplit_lines(self):
        """ライン接合後、Shape_lengthやアイテム数をテキスト出力します"""
        try:
            #フィーチャ数のLengthの集計
            moto_line_length,saki_line_length = self.count_length()

            #フィーチャ数のカウント
            result_in_feature,result_out_feature = self.count_feature()

            file_dir = self.workspace[:self.workspace.rfind("\\")]

            with open(os.path.join(file_dir,"export_upsplit_lines.txt"), "w") as f:

                message = u"接合元フィーチャの総延長=" + str(moto_line_length)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合先フィーチャの総延長=" + str(saki_line_length)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合元フィーチャ数=" + str(result_in_feature)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合先フィーチャ数=" + str(result_out_feature)
                f.write(message.encode("Shift-JIS"))
                f.write("\n")

                message = u"接合したラインのOID"
                f.write(message.encode("Shift-JIS"))
                f.write("\n")
                for line in self.line_list:
                    for l in line:
                        f.write(str(l))
                        f.write(" ")
                    f.write("\n")

        except IOError:
            print u"ファイルの出力に失敗しました。"
        except TypeError:
            print u"ファイルに出力する内容が不正です。"

    def count_feature(self):
        """接合元と接合先のフィーチャクラスのフィーチャ数をカウントします"""
        result_in_feature = arcpy.GetCount_management(self.in_feature)
        result_out_feature = arcpy.GetCount_management(self.out_feature)

        return result_in_feature,result_out_feature

    def count_length(self):
        """接合元と接合先のフィーチャクラスのShape_Lengthを集計します"""
        moto_line_length = float()
        saki_line_length = float()

        for line in [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]:
            moto_line_length += line[1].length

        for line in [row for row in arcpy.da.SearchCursor(self.out_feature, ["OID@","SHAPE@"])]:
            saki_line_length += line[1].length

        return moto_line_length,saki_line_length

    def set_join_oid(self,oid_list):
        """oid_listに格納されているOIDを分解します。[1,2,3]の場合、1,2,3に分解します。"""
        oid_list = set(oid_list)
        oid_list = list(oid_list)
        oid = ",".join([str(x) for x in oid_list])
        return oid

    def copy_feature(self):
        """接合対象ではないラインの作図を行います"""

        try:
            oid = self.set_join_oid(self.finish_line_list)

            with  arcpy.da.InsertCursor(self.out_feature, ["OID@","SHAPE@",UnsplitLine.column_unsplit_length]) as cur:
                for line in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"],"OBJECTID not in (%s)" % oid) :
                    length = self.get_shape_length(line[0])
                    for part in line[1]:
                        array = arcpy.Array()
                        for pnt in part:
                            if pnt:
                                array.add(arcpy.Point(pnt.X, pnt.Y))

                    cur.insertRow((line[0],arcpy.Polyline(array),length))

        except TypeError:
            print traceback.format_exc()

    def get_lines(self,oid=0):
        """接合元のラインを取得する関数です

        returns:
        self.moto_lines -- 接合元のラインの全フィーチャ
        line            -- ラインオブジェクト
        """

        if oid == 0:
            self.moto_lines = [row for row in arcpy.da.SearchCursor(self.in_feature, ["OID@","SHAPE@"])]
            return self.moto_lines
        else:
            line = [row for row in arcpy.da.SearchCursor(in_feature, ["OID@","SHAPE@"],"OBJECTID = %s" % oid)]
            line = [x for x in line for x in x]
            return line

    def search_line_combination(self,line2):
        """接合するラインの組み合わせを取得する関数です

            例:OBJECTID=1,2と4,5,6と10,11,12,13が接続対象の場合
               [[1,2],[4,5,6],[10,11,12,13]]というリストを作成します。

        arguments:
        line2 -- search_lineでのループから渡されるLineオブジェクト

        """

        point_list = []  #端点格納用リスト
        oid = line2[0]

        try:
            for line in self.moto_lines:
                if line[0] != oid:
                    #組み合わせ登録が完了したリストと端点が3か所重なっているライン(接合対象外)のリスト以外
                    if not line[0] in self.finish_line_list:
                        if (line[1].firstPoint.X == line2[1].lastPoint.X and line[1].firstPoint.Y == line2[1].lastPoint.Y) or \
                           (line[1].lastPoint.X == line2[1].lastPoint.X and line[1].lastPoint.Y == line2[1].lastPoint.Y) or \
                           (line[1].lastPoint.X == line2[1].firstPoint.X and line[1].lastPoint.Y == line2[1].firstPoint.Y) or \
                           (line[1].firstPoint.X == line2[1].firstPoint.X and line[1].firstPoint.Y == line2[1].firstPoint.Y):

                            point_list.append(line[0])
                            line_copy = line

            #隣り合っているラインが無い、もしくは、接合対象ラインの組み合わせリストができた場合のみこの分岐に入る
            if len(point_list) == 0:

                #接合対象ラインの組み合わせリストが存在する場合
                if len(self.match_point_list) != 0:
                    self.match_point_list = set(self.match_point_list) #上記処理だとOIDが重複してしまう可能性があるので、重複削除
                    self.line_list.append(list(self.match_point_list))
                    self.match_point_list = []

            #一番端のラインの場合はこの分岐に入る
            elif len(point_list) == 1:
                point_list.append(oid)
                self.finish_line_list.extend(point_list)
                self.match_point_list.extend(point_list)

                #再検索
                self.search_line_combination(line_copy)

            #両端にラインが存在する場合はこの分岐に入る(この分岐では組み合わせ登録は行わない)
            elif len(point_list) == 2:

                oid_list = []

                #端点が三つ以上重なっているか確認
                ret,point_list = self.check_line_combination(oid,point_list)

                #端点が三つ以上重なっている場合
                if ret == -1:
                    self.divide_line_list.extend(point_list)
                else:
                    oid_list.append(oid)
                    self.match_point_list.extend(oid_list)

                    #ループしている場合はこの分岐に入る
                    if len(self.match_point_list) > len(point_list):
                        self.finish_line_list.extend(self.match_point_list)

                    #重複していない要素を抽出
                    l = list(x for x in point_list if x not in self.match_point_list)

                    if len(l) == 1:
                        line = self.get_lines(l[0])
                        self.search_line_combination(line)
                    else:
                        self.search_line_combination(line_copy)

        except Exception:
            print traceback.format_exc()

    def check_line_combination(self,oid,point_list):
        """接合対象のラインにおいて、端点が三つ以上重なっているか確認する関数です

        arguments:
        oid -- search_lineでのループから渡されるOBJECTID
        point_list -- search_line_combinationで抽出した端点のリスト

        returns:
        1 or -1 -- 1:端点が三つ重なっていない -1:端点が三つ重なっている
        point_list -- point_list

        """

        point_list.append(oid)
        point_list = set(point_list)
        point_list = list(point_list)
        oid = ",".join([str(x) for x in point_list])

        devide_list = []

        #対象のOID(複数)で検索
        for row in arcpy.da.SearchCursor(in_feature, ["OID@", "SHAPE@"],"OBJECTID in (%s)" % oid):

            devide_list.append([row[0],row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y])

        #データフレームに格納
        df_devide = pd.DataFrame(devide_list)
        df_devide.columns = ["OID","FirstX","FirstY","LastX","LastY"]

        #始点と終点でグルーピング
        dic_first = df_devide.groupby("FirstX").groups
        dic_last = df_devide.groupby("LastX").groups
        list_first = [x for _,x in dic_first.iteritems()]
        list_last = [x for _,x in dic_last.iteritems()]

        #始点と終点で同じ座標が三つ以上ある場合は処理対象外とする(すべて始点、すべて終点の場合)
        for f,l in zip(list_first,list_last):
            if len(f) >= 3 or len(l) >= 3:
                return -1,point_list

        #始点と終点で同じ座標が三つ以上ある場合は処理対象外とする(始点と終点がそれぞれ3点に集中している場合)
        divide_line_list = []
        for key,row in df_devide.iterrows():
            for key2,row2 in df_devide.iterrows():
                if key !=  key2:
                    if row.FirstX == row2.LastX or row.FirstX == row2.FirstX:
                        divide_line_list.append(row.OID)
            if len(divide_line_list) >= 2:
                return -1,point_list

            divide_line_list = []

        return 1,point_list

    def search_line(self):
        """search_line_combinationでの結果を返す関数です

        returns:
        self.line_list -- 接合対象ラインのリスト(多次元)

        例:OBJECTID=1,2と4,5,6と10,11,12,13が接続対象の場合
        [[1,2],[4,5,6],[10,11,12,13]]というリストが返されます。

        """

        lines = self.get_lines()

        for line in lines:

            self.search_line_combination(line)

        return self.line_list

    def pop_list(self,pop_list):
        """リストから要素をポップする関数です

        arguments:
        pop_list -- ポップ対象のリスト

        """
        try:
            for p in pop_list:
                i = pop_list.index(p)
                self.moto_lines.pop(i)
        except ValueError:
            print traceback.format_exc()
            raise Exception

    def get_coord(self,oid):
        """接合対象ラインにおいて、各頂点の座標を取得する関数です

        arguments:
        oid -- OBJECTID → 接合対象のOBJECTIDが1,2,3 の場合、oid = 1,2,3 という形になります

        returns:
        coord_list -- 接合対象ラインの、各頂点の座標(多次元)
                     [0] = OBJECTID
                     [1] = ポイントごとに振ったID
                     [2] = X
                     [3] = Y
                     [4] = 始点X
                     [5] = 始点Y
                     [6] = 終点X
                     [7] = 終点Y

        loop_list -- 接合対象ラインがループしているのかどうかの判定用リスト(多次元)
                    [0] = OBJECTID
                    [1] = 始点X
                    [2] = 始点Y
                    [3] = 終点X
                    [4] = 終点Y

        """

        coord_list = []
        loop_list = []

        #対象のOID(複数)で検索
        for row in arcpy.da.SearchCursor(self.in_feature, ["OID@", "SHAPE@"],"OBJECTID in (%s)" % oid):

            loop_list.append([row[0],row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y]) #ループしたラインの判定用に使用する

            for part in row[1]:

                for i,pnt in enumerate(part):

                    if pnt:
                        #ポイントにIDを振る
                        pnt.ID = i

                        #ラインのOIDと座標を格納
                        coord_list.append([row[0],pnt.ID,pnt.X,pnt.Y,row[1].firstPoint.X,row[1].firstPoint.Y,row[1].lastPoint.X,row[1].lastPoint.Y])

        return coord_list,loop_list

    def check_geometry(self,df):
        """接合対象ラインにおいて、各頂点が始点なのか終点なのか判定する関数
           判定した結果を、データフレームに結合します。

        arguments:
        df -- Pandasのデータフレーム
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]


        returns:
        df_merge -- Pandasのデータフレーム
                    columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
                    PointType 1:始点 2:終点 None:それ以外

        """

        try:
            list_point_type = []

            #各頂点が始点なのか終点なのかを確認
            for _ ,row in df.iterrows():
                for _ ,row2 in df.iterrows():
                    if row.OID == row2.OID and row.ID == row2.ID:
                        if row.X == row2.FirstX and row.Y == row2.FirstY:
                            #始点として登録
                            list_point_type.append([row2.OID,row2.ID,"1"])
                        elif row.X == row2.LastX and row.Y == row2.LastY:
                            #終点として登録
                            list_point_type.append([row2.OID,row2.ID,"2"])

            #リストをデータフレームに格納
            df_point_type = pd.DataFrame(list_point_type)
            df_point_type.columns = ["OID","ID","PointType"]
            df_point_type = pd.DataFrame(df_point_type.drop_duplicates(["OID","PointType"]))

            #データフレームを結合
            df_merge = pd.merge(df, df_point_type, on=["OID","ID"],how='left')

            return df_merge
        except Exception:
            print traceback.format_exc()

    def check_first_last(self,df):
        """接合対象ラインにおいて、端のラインと端が始点なのか終点なのかを判定する関数です

        arguments:
        df -- Pandasのデータフレーム。
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
              PointType 1:始点 2:終点 None:それ以外

        returns:
        list_first_last -- 端のラインとPointTypeを格納したリスト
                           [0]:OBJECTID
                           [1]:頂点ごとに振ったID
                           [2]:PointType 1:始点 2:終点
        """

        list_first_last = []
        list_edge = []

        #ラインのつながりが二本の場合
        if len(df.groupby("OID").groups.keys()) == 2:
            for _ , row in df.iterrows():
                for _ , row2 in df.iterrows():
                    if row.OID !=  row2.OID:
                        if row.PointType == "1":
                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_first_last.append([row.OID,row.ID,"1"])

                        elif row.PointType == "2":

                            if (row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                               (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY):

                                list_first_last.append([row.OID,row.ID,"2"])

                        else:

                            continue

            return list_first_last

        else:
            for _ , row in df.iterrows():
                for _ , row2 in df.iterrows():
                    if row.OID != row2.OID:
                        if row.PointType == "1":
                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"1"])

                            if (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY) or \
                               (row.FirstX == row2.LastX and row.FirstY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"2"])

                        elif row.PointType == "2":
                            if (row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                               (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY):

                                list_edge.append([row.OID,row.ID,"2"])

                            if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                               (row.LastX == row2.LastX and row.LastY == row2.LastY):

                                list_edge.append([row.OID,row.ID,"1"])

                #重複を削除([[1,2,3],[1,2,3]]があったら[[1,2,3]]にする)
                c = Counter(tuple(x) for x in list_edge)
                list_edge = [list(k) for k,v in c.items()]

                if len(list_edge) >= 2:
                    list_edge = []

                elif len(list_edge) == 1:
                    list_first_last.append([list_edge[0][0],list_edge[0][1],list_edge[0][2]])

                    return list_first_last

    def counter(self):
        """頂点の作図順用のカウンタ"""
        n = [0]
        def inc():
            n[0] += 1
            return n[0]
        return inc

    def search_dataframe(self,df,oid,line_id):
        """データフレームのレコードを抽出する用の関数"""
        df = pd.DataFrame(df[df.OID == oid])
        df = pd.DataFrame(df[df.ID == line_id])

        return df

    def draw_point_sort(self,df,oid,line_id,point_type):
        """接合対象ラインにおいて、各頂点の作図順を設定する関数です
      self.list_draw_point_sort に作図順を格納して後続処理で使用します

        arguments:
        df -- Pandasのデータフレーム。
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointType"]
              PointType 1:始点 2:終点 None:それ以外
        oid -- OBJECTID(端のライン)
        line_id -- 各頂点に割り当てたID(端のライン)
        point_type -- 1:始点 2:終点 ←左記以外の値(None)は受け取らない想定(端のラインのどちらかの端点が入っているため)

        """

        #データフレームから対象のデータを抽出する(必ず一件になるはず)
        df_draw_point_sort = self.search_dataframe(df,oid,line_id)

        #対象のデータのX,Y,ポイントタイプ(始点or終点)
        first_x = df_draw_point_sort["FirstX"].values[0]
        last_x = df_draw_point_sort["LastX"].values[0]
        first_y = df_draw_point_sort["FirstY"].values[0]
        last_y = df_draw_point_sort["LastY"].values[0]

        if point_type == "":
            point_type = df_draw_point_sort["PointType"].values[0]

        #始点の場合
        if point_type == "1":
            #データフレームから対象のデータを抽出する(必ず複数件になるはず)
            df_asc = pd.DataFrame(df[df.OID == oid])

            #ポイントIDの昇順にソート(つまり、始点側から順にソートする)
            df_asc = df_asc.sort(['ID'], ascending=[True])

            #対象のラインの作図順を決定する
            for _ ,row in df_asc.iterrows():
                self.list_draw_point_sort.append([row.OID,row.ID,self.count()])

        #終点の場合
        elif point_type == "2":
            #データフレームから対象のデータを抽出する(必ず複数件になるはず)
            df_desc = pd.DataFrame(df[df.OID == oid])

            #ポイントIDの降順にソート(つまり、終点側から順にソートする)
            df_desc = df_desc.sort(['ID'], ascending=[False])

            #対象のラインの作図順を決定する
            for _ ,row in df_desc.iterrows():
                self.list_draw_point_sort.append([row.OID,row.ID,self.count()])

        #対象のラインに接続するラインで再度処理を行う
        for _ ,row in df.iterrows():
            if row.OID != oid:
                #対象のラインを始点から作図する場合は隣のラインは必ず対象のラインの終点に一致する
                if point_type == "1":

                    if row.X == last_x and row.Y == last_y:

                        if not row.OID in self.finish_oid:

                            self.finish_oid.append(oid)
                            self.draw_point_sort(df,row.OID,row.ID,"") #再検索(更に隣のラインを検索)

                #対象のラインを終点から作図する場合は隣のラインは必ず対象のラインの始点に一致する
                elif point_type == "2":
                    if row.X == first_x and row.Y == first_y:

                        if not row.OID in self.finish_oid:
                            self.finish_oid.append(oid)
                            self.draw_point_sort(df,row.OID,row.ID,"") #再検索(更に隣のラインを検索)

    def search_loop_line(self,df_loop):
        """接合対象ラインがループしているラインなのか判定する関数

        arguments:
        df_loop -- Pandasのデータフレーム。
                  columns = ["OID","FirstX","FirstY","LastX","LastY"]

        returns:
        1 or -1-- 1:ループしたラインではない -1:ループしたライン

        """

        #self.finish_line_listに格納されていないものはここで格納
        for _,row in df_loop.iterrows():
            if not row.OID in self.finish_line_list:
                self.finish_line_list.append(row.OID)

        #接合対象のラインが二本の場合
        if len(df_loop) == 2:

            for _,row in df_loop.iterrows():
                for _,row2 in df_loop.iterrows():
                    if row.OID !=  row2.OID:
                        if ((row.FirstX == row2.LastX and row.FirstY == row2.LastY) or \
                           (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY)) and \
                           ((row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                           (row.LastX == row2.LastX and row.LastY == row2.LastY)):

                            break

                #ループしたラインではない場合
                else:
                    return 1

            #ループしたラインの場合(breakしたらこっちに入る)
            return -1

        elif len(df_loop) >= 3:
            oid_list = []

            #ループしているラインを判定(接合対象ラインにおいて、全始点終点が他ラインの始点終点と一致した場合はループしている)
            for _,row in df_loop.iterrows():
                for _,row2 in df_loop.iterrows():
                    if row.OID !=  row2.OID:
                        if (row.FirstX == row2.FirstX and row.FirstY == row2.FirstY) or \
                           (row.FirstX == row2.LastX and row.FirstY == row2.LastY):

                            oid_list.append(1)

                        if (row.LastX == row2.FirstX and row.LastY == row2.FirstY) or \
                           (row.LastX == row2.LastX and row.LastY == row2.LastY):

                            oid_list.append(2)

                #ループしてないラインの場合はこの分岐に入る
                if len(oid_list) != 2:
                    return 1

                oid_list = []

            #ループしたラインの場合
            return -1

    def get_shape_length(self,oid):
        """接合対象ラインのShape_lengthの合計を取得する関数"""

        total_length = float()
        for line in arcpy.da.SearchCursor(self.in_feature, ["SHAPE@"],"OBJECTID in (%s)" % oid) :
             total_length += line[0].length

        return total_length

    def upsplit_line(self,cur,df,oid):
        """ラインを接合する関数

        arguments:
        cur -- カーソル
        df -- Pandasのデータフレーム
              columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY","PointSort"]
        oid -- OBJECTID(1,2,3が接合対象の場合oidには1,2,3 が格納されている)

        """

        try:
            array = arcpy.Array()

            length = self.get_shape_length(oid)

            for d in df.index:
                array.add(arcpy.Point(df.ix[d][2], df.ix[d][3]))
            cur.insertRow((arcpy.Polyline(array),length))
        except Exception:
            print traceback.format_exc()

    def create_line(self,lines):
        """各種関数を呼び出す処理。ここで、ラインの接合や接合対象ではないラインのコピー関数を呼び出す。

        arguments:
        lines -- search_lineで格納したリスト(多次元)

                 例:OBJECTID=1,2,3 また、4,5 というラインが接合対象の場合、
                    [[1,2,3],[4,5]] という要素が格納されている

        """

        try:
            if len(lines) == 0:
                print u"対象データがありません"
                raise Exception

            with  arcpy.da.InsertCursor(self.out_feature, ["SHAPE@",UnsplitLine.column_unsplit_length]) as cur:

                #接合対象ラインのリストでループ
                for line in lines:

                    #接合対象のラインのリストを分解(例 [1,2,3]→1,2,3)
                    oid = self.set_join_oid(line)

                    #接合対象ラインの座標をリスト化する
                    coord_list,loop_list = self.get_coord(oid)

                    df = pd.DataFrame(loop_list)
                    df.columns = ["OID","FirstX","FirstY","LastX","LastY"]

                    #対象のラインがループしているか判定
                    result = self.search_loop_line(df)

                    #ループしたラインがない場合
                    if result == 1:

                        #データフレーム作成
                        df_coord = pd.DataFrame(coord_list)
                        df_coord.columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]

                        #各頂点を始点、終点、それ以外に仕分け
                        df_checked = self.check_geometry(df_coord)

                        #各ラインのつながりにおいて、端のラインと端点を抽出する
                        list_checked = self.check_first_last(df_checked)

                        #上記で端のラインが判定できなかった場合、エラーとする
                        if not list_checked:
                            print u"端点が判断できないラインが存在しました" + " " + "OID=" + oid
                            raise Exception

                        #接合対象のラインにおいて、端のライン
                        oid_one,line_id,point_type = list_checked[0][0],list_checked[0][1],list_checked[0][2]

                        #隣りあっているラインを抽出して作図順を決めていく
                        self.draw_point_sort(df_checked,oid_one,line_id,point_type)

                        #リストをデータフレーム(各頂点の作図順を格納している)に格納
                        df_point_sort = pd.DataFrame(self.list_draw_point_sort)
                        df_point_sort.columns = ["OID","ID","PointSort"]

                        #データフレームを結合
                        df_checked = pd.merge(df_checked, df_point_sort, on=["OID","ID"],how="left")

                        #重複を削除
                        df_checked = pd.DataFrame(df_checked.drop_duplicates(["X","Y"]))

                        #各頂点の作図順でソート
                        df_checked = df_checked.sort(["PointSort"], ascending=[True])

                        #ラインを作図
                        self.upsplit_line(cur,df_checked,oid)

                        self.list_draw_point_sort = []

                    #ループしたラインがある場合
                    else:

                        #データフレーム作成
                        df_coord = pd.DataFrame(coord_list)
                        df_coord.columns = ["OID","ID","X","Y","FirstX","FirstY","LastX","LastY"]

                        #各頂点を始点、終点、それ以外に仕分け
                        df_checked = self.check_geometry(df_coord)

                        #隣りあっているラインを抽出して作図順を決めていく
                        self.draw_point_sort(df_checked,coord_list[0][0],coord_list[0][1],"")

                        #リストをデータフレーム(各頂点の作図順を格納している)に格納
                        df_point_sort = pd.DataFrame(self.list_draw_point_sort)
                        df_point_sort.columns = ["OID","ID","PointSort"]

                        #データフレームを結合
                        df_checked = pd.merge(df_checked, df_point_sort, on=["OID","ID"],how="left")

                        #各頂点の作図順でソート
                        df_checked = df_checked.sort(["PointSort"], ascending=[True])

                        #ラインを作図
                        self.upsplit_line(cur,df_checked,oid)

                        self.list_draw_point_sort = []

            #非接合対象のラインをコピーする
            self.copy_feature()

            #接合したラインの情報をテキストで出力
            self.export_upsplit_lines()

        except Exception:
            print traceback.format_exc()

if __name__ == '__main__':

    #対象のgdb
    workspace = r"C:\ArcPySample\ArcPyAdvanced.gdb"

    #接合元のフィーチャクラス名
    in_feature = "LineSplit"

    #接合先のフィーチャクラス名
    out_feature = "LineUnsplit"

    #インスタンス化
    unsplitLine = UnsplitLine(workspace,in_feature,out_feature)

    #接合対象のラインを抽出
    lines = unsplitLine.search_line()

    #対象のラインを接合
    unsplitLine.create_line(lines)


このソースをArcMap上で実行した結果が以下になります。
f:id:sanvarie:20160319171524p:plain


ちゃんと作られました。このように一本ずつだったラインが
f:id:sanvarie:20160319171549p:plain


このように接合されました!
f:id:sanvarie:20160319171605p:plain


プログラムの中で「UNSPLIT_LENGTH」というカラムを作成し、接合前のShape_Lengthを格納するようにしました(接合対象のラインが三本あり、それぞれの長さが10,20,30だった場合、60が格納される)。ただ、ちょっと誤差がでてしまっているので、接合結果を比較する場合は、小数点以下をはずすなどしないとだめそうですね・・・
f:id:sanvarie:20160319171709p:plain


出力前後の結果確認用にテキストも出力するようにしています(今回はC:\ArcPySampleの直下に作成)。こちらも微妙に誤差が出ていますね・・・
f:id:sanvarie:20160319171839p:plain

課題

  1. 接合前後でShape_Lengthに微妙に誤差がでてしまっている
  2. 大量のデータを処理する場合、パフォーマンスに難あり
  3. ソースが煩雑になりすぎた(ちょっと色々無理矢理やってしまったところがあります・・・)。共通化など、色々整理が必要ですね。
  4. バグが大量にありそうです・・・

課題は多いものの、とりあえずAdvancedの機能であるライン接合処理ができるようになったことは素直にうれしいですね。自分も色々勉強しながらだったので、今度作るときはもっとうまくプログラムを書くことができそうです。

以上、本日は「【PythonとGIS】ArcPyでライン接合を行ってみた ~Advancedライセンス?!だが断る!~」でした。

ArcPyを使用してラインの頂点数を取得する

さて、本日はArcPyを使用してラインの頂点数を取得してみようと思います。頂点数が多いとパフォーマンスにも影響が出てしまうため、余分な頂点数がありそうなラインを抽出するのに役立つかと思います。

実行環境

事前準備

  1. 「Line」フィーチャクラスを作成します。「頂点数」格納用のカラムも用意します。
  2. 「Line」フィーチャクラスにラインを作図します。

f:id:sanvarie:20160302133304p:plain

f:id:sanvarie:20160302142157p:plain

サンプルコード

ラインの頂点数を取得して属性テーブルに格納するサンプルです。

# -*- coding: utf-8 -*-
import arcpy

#対象のGDB
arcpy.env.workspace = "C:\ArcPySample.gdb"

#対象のフィーチャクラス
infc = "Line"

dictTop = {}
cnt = 0

#ラインフィーチャの頂点数を検索
for row in arcpy.da.SearchCursor(infc, ("OID@","SHAPE@")):
    for part in row[1]:
        for pnt in part:
            if pnt:
                cnt += 1
            else:
                print("error")
                continue
        #ディクショナリにOIDと頂点数を格納
        dictTop[row[0]] = cnt
        cnt = 0
del row

#各フィーチャの頂点数を更新
for i in dictTop:
    cursor = arcpy.da.UpdateCursor(infc, ("OID@","Top"))
    for row in cursor:
        if i == row[0]:
            row[1] = dictTop[i]
            cursor.updateRow(row)
del row

こんな感じになれば成功です。今回は属性に頂点数を格納しましたが、CSV出力したりしてもいいかもしれませんね。
f:id:sanvarie:20160302142614p:plain

簡単ではありますが、本日は以上です。