ポリゴン ジオメトリの例
複数の三角形があり、一部には穴が開いています。緑色の数字はポリゴン外周リングの ID 値を示します。赤色の数字は、外周リングと内周リングに付けられた連続するポイント ID 値です。
ID 値が次のようになっている場所があることは、見た目からすぐにわかります。
- 一度だけ現れるもの（線分の接続点）、
- 2つの番号があるもの（ポリゴンの始点と終点）、
- 関連付けられた ID が3つ、4つ、あるいは5つあるもの（さまざまな形でポリゴンが接する場所）。
arcpy から numpy へ
FeatureClassToNumPyArray を使ってジオメトリを NumPy 配列形式で出力する場合、いくつかの方法があります。私の場合は、オブジェクト ID と X、Y 座標の値を別々に出力するのが好みです。属性（テキストやその他の数値）を扱う必要がある場合は、別途出力します。これは再構築した出力に簡単に結合できます。
基本的な構文は次のとおりです。
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # 必要
import arcpy
from arcpy.da import FeatureClassToNumPyArray
in_fc = r"C:\path_to_gdb\the_geodatabase.gdb\featureclass_name"
flds = ['OID@', 'SHAPE@X', 'SHAPE@Y'] # -- 使用するフィールド
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
もちろん、処理をわかりやすくするため、ローカルに保存されたデータを使っているはずです。用途に合わせて、上記の3行目を変更してください。
4行目では、OBJECTID フィールドとポリゴン頂点の X、Y 値を取り込むために使用するフィールド トークンを指定しています。
私が使っているルールは次のとおりです。
- 入力フィーチャクラスに OGC ジオメトリ検証を適用し、穴のあるポリゴンでは、穴（内周リング）が親となる外周リングの後に続く順序になっていることを確認します。穴のジオメトリが外周リングより前にあるユースケースは思い当たりません。これにより、ポイント イン ポリゴン判定で穴のリングを外周リングに割り当てる手間を省けます。
- 投影座標系のデータを使い、少なくとも正角図法（UTM、MTM）、または場合によっては正積図法を使用します。ジオメトリ属性の判定には役に立たないウェブ向け投影法は、絶対に使いません。
- 外周リングは時計回り、穴（内周リング）は反時計回りです。
- ポリゴンの始点と終点は X、Y 値が重複しますが、ID 番号は異なります。
- さらに、ESRI が指定するその他すべての要件も満たします。
関連情報はジオメトリのチェックとジオメトリの修復を参照してください。
基本的な配列出力
では、配列の生成に戻りましょう。
# ---- 上記から続く
ids = out['OID@'].astype('int') # -- 値を単純な numpy ndarray として取得
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- 構造化配列を ndarray に変換
LL = np.min(xys, axis=0)
arr = xys - LL
#
ids
[1 1 1 1 1 1 1 1 1 ... 5 5 5 5]
arr
[[ 10.0, 0.0]
[ 20.0, 20.0]
[ 30.0, 0.0]
[ 10.0, 0.0]
[ 17.5, 5.0]
[ 22.0, 5.0]
[ 22.0, 10.0]
[ 17.5, 10.0]
[ 17.5, 5.0]
# ... 省略
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
2行目では、派生した構造化配列からオブジェクト ID の値を取り出し、単純な整数にしています。
3行目では同様に、X、Y 座標を double 型の numpy ndarray として取り出します（stuは最初のコードブロックでインポートし、名前を変更した numpy 関数です）。
4行目では、座標を扱いやすい平面空間に移動しています。すべての座標値から X、Y の最小値を差し引きます。この例では、座標 (0.0, 0.0) は実際には (300,000.00, 5,000,000.00) です。XY 平面を移動させることで、読みやすく、計算もしやすくなります。計算が終わったら座標を元に戻せばよいので、問題ありません。
7、8行目は、9個のポイントで外周リングとその穴を構成する最初の形状 (1) と、最後の形状 (5) に関連付けられた4個のポイントの ID 値です。10～24行目はポイントの座標です。11行目と14行目は外周リング（三角形）の始点と終点の座標で、15～19行目は四角形の穴を構成するポイントです。
ここからできること
この時点で、ポリゴン ジオメトリを表す複数のポイントが得られています。次のことができます。
- ポイントを移動する（X、Y 値に加算または減算するだけです）
- リング上のポイント数を増やす（密化する）、または減らす（単純化する）（そのための numpy 関数があります）。
- など
ただし、個々のポリゴンを扱う必要がある場合は、ポリゴンを分割しなければなりません。ここでは、ポリゴン自体や内周リング、外周リングを分離するためのいくつかの方法を紹介します。
def array_to_rings(arr, from_to_ids=False):
"""配列を外周リングと内周リングに分割します。
パラメーター
----------
arr : array_like
x、y 座標を含む Nx2 の配列／リスト／geo 配列。1つ以上の外周リング
および内周リング（穴）を含む場合があります。
戻り値
-------
配列のリスト。内周リング／穴がある場合は、入れ子になることがあります。
"""
# --
def _area_(a):
"""`_bit_area_` を変更したもの。面積が負の場合は穴です。`shoelace` 式とは
逆向きに計算します。
"""
x0, y1 = (a.T)[:, 1:] # 外積を次のように設定
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # 2024-03-28 変更
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""反時計回りかどうかを判定します。shoelace 式とは逆向きに計算します。"""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nNx2 の np.ndarray が必要です。"
if not isinstance(arr, np.ndarray):
print(msg)
return None
if hasattr(arr, "IFT"):
arr = np.copy(arr, subok=False)
if len(arr.shape) != 2:
print(msg)
return None
#
N = len(arr)
fst = 0
out = []
n = 0
c_o = 0 # -- 外周リングのカウンター
outer_ = []
inner_ = []
frto_ = [] # -- 入力配列 `arr` を分割するための開始・終了 ID
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\n一致するポイントがありません。入力配列に不備があります。")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # 部分配列のサイズ
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # 既存の外周リングに穴のリングを追加
out[-1].append(sub)
else: # 万一、外周リングより先に穴が見つかった場合
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- 配列のリスト、または何らかの形式の配列のリストのリストを返す
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
出力は次のとおりです。
array_to_rings(A5)
([[array([[ 10.000, 0.000],
[ 20.000, 20.000],
[ 30.000, 0.000],
[ 10.000, 0.000]]),
array([[ 17.500, 5.000],
[ 22.000, 5.000],
[ 22.000, 10.000],
[ 17.500, 10.000],
[ 17.500, 5.000]])],
[... 省略],
[... 省略],
array([[ 30.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 0.000],
[ 30.000, 0.000]]),
array([[ 40.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 20.000],
[ 40.000, 0.000]])]
上の出力では、最初の形状が2つの配列のリスト（外周リングと穴）になっています。最後の2つの配列は、単純な三角形である形状4と5のものです。
この時点で、穴を取り除いた出力を生成することもできます。
out = []
for i in arr:
if isinstance(i, list):
out.append(i[0])
else:
out.append(i)
out
[array([[ 10.000, 0.000],
[ 20.000, 20.000],
[ 30.000, 0.000],
[ 10.000, 0.000]]),
array([[ 30.000, 0.000],
[ 20.000, 20.000],
[ 40.000, 20.000],
[ 30.000, 0.000]]),
array([[ 10.000, 0.000],
[ 0.000, 20.000],
[ 20.000, 20.000],
[ 10.000, 0.000]]),
array([[ 30.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 0.000],
[ 30.000, 0.000]]),
array([[ 40.000, 0.000],
[ 35.000, 10.000],
[ 40.000, 20.000],
[ 40.000, 0.000]])]
その後、配列を arcpy のポリゴンに変換できます。
def array_poly(arrs, p_type=None, sr=None):
"""Geo 配列または ndarray からポリゴン フィーチャを作成します。
`geometry_fc` から使用するか、単独で使用できます。
パラメーター
----------
arr : arrays
ポイント配列。
p_type : text
POLYGON または POLYLINE
sr : spatial reference
空間参照オブジェクト、名前、または ID。
注記
-----
配列データからポリラインまたはポリゴン フィーチャを作成できます。
フィーチャは、内周リングの有無を問わずマルチパートにできます。
外周リングは時計回り、内周リング（穴）は
反時計回りに並べます。
"""
def _arr_poly_(arr, SR, as_type):
"""nan 値が現れる位置で配列をスライスして分割します。"""
aa = [Point(*pairs) for pairs in arr]
if as_type.upper() == 'POLYGON':
poly = Polygon(Array(aa), SR)
elif as_type.upper() == 'POLYLINE':
poly = Polyline(Array(aa), SR)
return poly
# --
polys = []
if isinstance(arrs, np.ndarray):
if len(arrs.shape) == 2:
arrs = [arrs]
elif arrs.dtype == 'O':
arrs = arrs
for a in arrs:
p = _arr_poly_(a, sr, p_type)
polys.append(p)
return polys
そして、上記を使ってフィーチャクラスに変換します。
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""`fc_geometry` で作成した配列リストからポリゴン フィーチャを作成します。
パラメーター
----------
a : array or list of arrays
通常は ``pnts_arr`` で作成されるオブジェクト配列の場合もあります。
p_type : string
POLYGON などの大文字表記のジオメトリ タイプ。
gdb : text
ジオデータベースのパスと名前。
fname : テキスト
フィーチャクラス名。
sr : 空間参照
名前またはオブジェクト
戻り値
-------
シングルパートおよび／またはマルチパートのフィーチャクラス。
備考
-----
処理は ``array_poly`` によって実行されます。
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # array_poly を呼び出す
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # 旧形式は in_memory
CreateFeatureclass(wkspace, fname, p_type, spatial_reference=sr)
AddField(fname, 'ID_arr', 'LONG')
with InsertCursor(fname, ['SHAPE@', 'ID_arr']) as cur:
for row in out:
cur.insertRow(row)
CopyFeatures(fname, name)
return None
これは、「ポリゴン パートの削除」ツールと同等です。不思議なことに、ArcGIS Pro ではこのツールを使うのに Advanced ライセンスが必要です。穴を削除／塗りつぶすだけなのに、なぜ Advanced ライセンスが必要なのでしょうか。
上記の出力タイプを POLYLINE に変更すると、「ポリゴンからライン」と同じように、ポリゴンをポリラインに変換できます。このツールにも Advanced ライセンスが必要です。
これは、FeatureClassToNumPyArray を使ってジオメトリを処理する方法を示す簡単な例です。
その他のツールとコードはこちらからご覧いただけます：NumPy GeometryおよびArcGIS Pro 向けツール。
