さて、本日は頂点の間引き処理を行ってみたいと思います。ArcMapの頂点の間引きというジオプロセシングツールを使えばそれを実現することが可能です。頂点の間引き (Generalize)—ヘルプ | ArcGIS for Desktop
ただ、以下のような問題もあります。
- Standardライセンスからのみ使用可能で、Basicまでしか持っていないという場合は使用不可
今回はBasicライセンスしか持っていないという場合でもそれに対応できるPythonスクリプトを作成してみました。
ソースはGistにも公開しています。https://gist.github.com/sanvarie/51bcdd9e1668081bde75
環境
windows7 64bit
ArcGIS10.2.0
Python2.7.3
事前準備
- GDBを作成する(今回は「Generalize.gdb」を作成しました)
- フィーチャクラスを作成する(今回は「Line」というラインのフィーチャクラスを作成しました)
処理概要
- フィーチャクラスを作成する(間引いたラインを作成するフィーチャクラス)
- 始点終点は間引き処理対象外とする
- 頂点間の角度と距離を計測して間引き対象/対象外を選別する
- 頂点間の角度の差異が1度以下、かつ、頂点間の長さが1m以上を間引き対象とする
対象データ
こういったデータを作成しました。
各頂点はこのようになっています。
■generalize.py
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]
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)
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))
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@"])]
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):
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):
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__':
workspace = r"C:\ArcPySample\Generalize.gdb"
in_feature = "Line"
out_feature = "LineGeneralize"
generalize = Generalize(workspace,in_feature,out_feature)
generalize.calc_point_distance()
結果を確認すると問題なく、頂点が間引きされていることがわかりました。
課題
- データ量が多い場合は少し時間がかかってしまう
- 属性が存在する場合は、間引いた先のフィーチャクラスにはその属性は引き継がれない(generalize_lineメソッドを修正すればいける)
課題は少しありますが、角度や距離をうまく使用すれば、ArcMapの機能よりも柔軟な使用ができそうですね。簡単ではありますが、本日は以上です。