GIS奮闘記

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

スポンサーリンク

【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の機能よりも柔軟な使用ができそうですね。簡単ではありますが、本日は以上です。