YouTube | Facebook | X(Twitter) | RSS

ドーナツ ポリゴンの穴を埋めるジオプロセシング ツール

2018/1/9 (火)

穴の空いたポリゴンをドーナツ ポリゴンと呼びますが、この穴を埋めるツールの紹介です。ラスターベクター変換を行うと不要な微小のドーナツ ポリゴンが生成されたので、それらを除外したいと思い作成しました。

はじめに

ArcMap や ArcGIS Pro の標準 GUI では、1フィーチャずつ図形を編集することでドーナツの穴を埋めることはできますが、フィーチャクラス毎に一括でドーナツ ポリゴンの穴を埋めるツールは用意されていません。

ArcObjects を使えばパート(ポリゴンを構成する要素)であるリング(パートのうち輪になっているもの)が外枠(エクステリア)用か内枠(インテリア)用かを判定するメソッドを使って対象リングを見つけて除外することができるのですが、ArcObjects によるジオプロセシング ツールの作成は容易ではありません。

ArcPy ではジオプロセシング ツールを作成しやすい反面、ポリゴンの構成要素であるリングがインテリアかエクステリアか判定できるようなきめ細かいメソッドは用意されていません。一長一短だと思って調べていたところ、ArcPy でも頑張ればインテリア リングの削除ができる方法を紹介している記事がありました(引用)。

その記事を引用して作成したスクリプト ツールがこれです。ドーナツ ポリゴンの穴を除去してくれます。閾値も設定できるので、一定面積以下の穴だけふさぐこともできます。閾値に 0 を指定するとすべての穴をふさぎます。

処理実行前のフィーチャクラス

処理実行後のフィーチャクラス

ダウンロード

参考サイトのコードを元に ArcGIS Pro 2.0 (Python 3.4) で作成したスクリプト ツールです。これからは極力 ArcGIS Pro を使っていきたいという意図です。ZIPファイルをダウンロードして解凍し、ArcGIS Pro の [カタログ] ウィンドウから *.tbx ファイルを参照し、スクリプト ツールを実行してください。コードはツール内に埋め込んだ状態にしていますがパスワード保護はしていないので取り出すこともできます。

なお、*.tbx ファイル自体は ArcMap (Python 2.7) と互換しない(と思う)のですが、コードは共通なので、ArcMap でスクリプト ツールを作成し、以下のコードをそのまま利用することができます。

ソースコード

#coding:cp932

import arcpy
import os

def RemovePolygonHoles_management(in_fc, threshold=0.0):
    """
    この関数はポリゴン フィーチャクラスから穴を削除します。
    閾値を設定すると閾値より小さい面積の穴のみが削除されます。
    閾値が設定されていない場合はすべての穴が削除されます。
    in_fc はポリゴン フィーチャクラスです。
    threshold は数値です。
    """
    desc = arcpy.Describe(in_fc)
    if desc.dataType !="FeatureClass" and desc.dataType != "ShapeFile":
        print "無効なデータです。入力はポリゴン フィーチャクラスのみサポートします。"
        return
    else:
        if desc.shapeType != "Polygon":
            print "無効なデータです。入力はポリゴン フィーチャクラスのみサポートします。"
            return
    if threshold < 0.0: threshold = 0.0 with arcpy.da.UpdateCursor(in_fc, ["SHAPE@"]) as updateCursor: for updateRow in updateCursor: shape = updateRow[0] new_shape = arcpy.Array() for part in shape: new_part = arcpy.Array() if threshold > 0:
                    #シェープのパートから None ポイントを探します。
                    #ArcPy モジュールでは、None ポイントを使用して外部と内部の頂点を区切ります
                    null_point_index = []
                    for i in range(len(part)):
                        if part[i] == None:
                            null_point_index.append(i)
                    #インテリア頂点が存在する場合、ポリゴンを作成してポリゴンの形状領域を所定の閾値と比較します。
                    #閾値より大きい場合は頂点を保持し、そうでない場合は閉じます。
                    if len(null_point_index) > 0:
                        for k in range(0, null_point_index[0]):
                            new_part.add(part[k])
                        for i in range(len(null_point_index)):
                            pointArray = arcpy.Array()
                            #None  ポイントが最後のポイントかどうかを判定
                            if i+1 < len(null_point_index): for j in range(null_point_index[i] + 1, null_point_index[i+1]): pointArray.add(part[j]) else: for j in range(null_point_index[i] + 1, len(part)): pointArray.add(part[j]) #与えられた閾値に対してシェープの面積をチェックするポリゴンを作成 inner_poly = arcpy.Polygon(pointArray) #閾値より大きい場合は新しいパート Array に追加 if inner_poly.area > threshold:
                                if i+1 < len(null_point_index): for k in range(null_point_index[i], null_point_index[i+1]): new_part.add(part[k]) else: for k in range(null_point_index[i], len(part)): new_part.add(part[k]) new_shape.add(new_part) #インテリアが存在しない場合は全体を追加 else: new_shape.add(part) else: #最初の None ポイント インデックスを取得 first_null_point_index = 0 for i in range(len(part)): if part[i] == None: first_null_point_index = i break if first_null_point_index == 0: new_shape.add(part) else: for j in range(first_null_point_index): new_part.add(part[j]) new_shape.add(new_part) if len(new_shape) > 0:
                new_poly = arcpy.Polygon(new_shape)
                updateRow[0] = new_poly
                updateCursor.updateRow(updateRow)

#スクリプト ツール
try:
    in_fc = arcpy.GetParameterAsText(0)		#入力フィーチャクラス
    out_fc = arcpy.GetParameterAsText(1)	#出力フィーチャクラス
    threshold = arcpy.GetParameterAsText(1)	#許容値

    #フィーチャクラスのコピー
    arcpy.CopyFeatureClass_management(in_fc, out_fc)

    #ドーナツ ポリゴンの穴埋め
    RemovePolygonHoles_management(out_fc, float(threshold))

except arcpy.ExecuteError as e:
    print("ツールでエラーが発生しました。")
    print(str(e).decode("UTF-8"))&a<span data-mce-type="bookmark" style="display: inline-block; width: 0px; overflow: hidden; line-height: 0;" class="mce_SELRES_start"></span>mp;amp;lt;span data-mce-type="bookmark" style="display: inline-block; width: 0px; overflow: hidden; line-height: 0;" class="mce_SELRES_start"></span>

引用

  • この記事を書いた人

羽田 康祐

伊達と酔狂のGISエンジニア。GIS上級技術者、Esri認定インストラクター、CompTIA CTT+ Classroom Trainer、潜水士、PADIダイブマスター、四アマ。WordPress は 2.1 からのユーザーで歴だけは長い。 代表著書『"地図リテラシー入門―地図の正しい読み方・描き方がわかる』 GIS を使った自己紹介はこちら。ESRIジャパン(株)所属、元青山学院大学非常勤講師を兼務。日本地図学会第31期常任委員。発言は個人の見解です。

-プログラミング, ArcGIS
-,