GIS奮闘記

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

スポンサーリンク

ArcPyレシピ集② ~ワークスペースのフィーチャクラス名とジオメトリタイプの取得方法~

さて、GWも終わり憂鬱な気分となっておりますが、本日もArcPyレシピを書いてみたいと思います。

本日のお題

おそらく業務で使用するワークスペース(FGDBなど)には多数のフィーチャクラスが存在しているのではないでしょうか。状況によってはそれらの名称とジオメトリタイプを取得したい場合があるかと思います。今回はこれを題材にしたいと思います。

対象のワークスペースとフィーチャクラス

以下のようなサンプルデータを作成しました。

ワークスペース
・WorkSpaceSample.gdb

②フィーチャクラス
・HOUSE・・・ポリゴン
・HOUSE_ANO・・・アノテーション
・ROAD・・・ライン
・SHOP・・・ポイント

これらをディクショナリで取得しようと思います。{'ROAD': 'Polyline', 'SHOP': 'Point'}のような感じですね。

f:id:sanvarie:20160509124855p:plain

環境

ArcGIS10.2
Python2.7.3

サンプル

今回も簡単ですね。

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

def get_all_featureclass_type():

    geo_type_dict = {}
    fc_list = arcpy.ListFeatureClasses()

    for fc in fc_list:

        geo_type_dict[fc] = get_featureclass_type(fc)

    return geo_type_dict

def get_featureclass_type(fc):

    fc_type = arcpy.Describe(fc).FeatureType

    #アノテーションのジオメトリタイプはポリゴンになってしまうのでこの処理を追加
    if fc_type == "Annotation":

        geo_type = "Annotation"

    else:

        geo_type = arcpy.Describe(fc).shapeType

    return geo_type

arcpy.env.workspace = ""  #gdbを指定
print get_all_featureclass_type()

以下のような結果になりました。
f:id:sanvarie:20160509125809p:plain

解説

①まず、「ListFeatureClasses()」で「arcpy.env.workspace」で指定したワークスペースに存在するフィーチャクラスのリストを取得します(上記「②フィーチャクラス」参照)。
②取得したフィーチャクラスごとにループして「get_featureclass_type」を呼び出しています。
③「get_featureclass_type」内でジオメトリタイプを取得しています。

※1.アノテーションのジオメトリタイプはポリゴンとして認識されてしまうので、「Annotation」という文字列を代入しています。
※2.「Describe」とは、地理データのプロパティにアクセスする手段です。例えばフィーチャクラスの「Describe」は以下です。

f:id:sanvarie:20160509130545p:plain

「Describe」はけっこう使える関数なので、ぜひ積極的に使ってみてください!簡単ですが、本日は以上です。

ArcPyレシピ集① ~フィーチャクラスのカラム名の取得方法~

さて、本日から簡単なArcPyのレシピ集を書いてみたいと思います。もっといい方法があるよ!という方はご連絡いただけると助かりますm(_ _)m

本日のお題

フィーチャクラスを扱うにおいて、カラム名をとってきたいことはよくあると思います。今回はこれを題材にしたいと思います。

対象フィーチャクラス

こんな感じのフィーチャクラスを作ってみました。

f:id:sanvarie:20160503160851p:plain

環境

ArcGIS10.2
Python2.7.3

サンプル

超簡単です。フィーチャクラスのカラム名をリストにつっこんでいます

def get_column_name(in_feature):
    """フィーチャクラスのカラム名をリストで返します"""

    fields = arcpy.ListFields(in_feature)
    return [field.name for field in fields]

arcpy.env.workspace = ""      #gdbを指定
columns = get_column_name("") #フィーチャクラスを指定

print columns

以下のような結果になりました。超簡単ですね!

f:id:sanvarie:20160503161140p:plain

解説

解説するほどでもありませんが、arcpy.ListFields()で対象フィーチャクラスのFieldオブジェクトを取得することができます。Fieldオブジェクトのプロパティは以下のようになっています。

f:id:sanvarie:20160503161546p:plain

「name」というプロパティがありますね。今回はそれをリスト内包表記で取得しています。本当にこれだけですね。
ListFieldsだけではなく、ListFeatureClasses(対象のワークスペースに存在するフィーチャクラスを取得するメソッド)とか色々あるので、ぜひ使ってみてください。

このレシピ集、どこまで続けられるか(ブログ主の知識不足)。。簡単ですが、本日は以上です!

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

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

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