かつてはジオメトリの三角分割が「流行」でした。
3D Analyst拡張機能があれば可能ですが、持っていなくて原理を知り、NumPyとpythonを使えるなら…さあ始めましょう。
1 ジオメトリから配列を作成する（FeatureClassToNumPyArrayが出発点です。以前にその方法のコードを投稿しました）
2 それを三角分割する（私はscipyを使います）
3 各三角形の重心を求める（これもコードを投稿済みです）
4 重心が元の形状内にあるか確認する。
ここで複雑になるのは、多角形は凸または凹で、穴がある場合もない場合もあり、複数のパーツで構成されることもあります。だから時にはジオメトリを最も単純な形に分解して繰り返す必要があります。凸形状は退屈なので無視します。
ここに3つの形状、C、D、Aがあります。Cは穴なしの凹形状、Dは穴ありの凸形状、Aは穴ありの凹形状です。
それらの三角分割されたバージョンです。
そして制約付き三角分割。つまり、元の外側境界や穴の中にない三角形による三角分割です。
すべて問題なさそうです。
三角分割を実行するには、次を使用できます
from scipy.spatial import Delaunay
def triangulate_pnts(pnts):
"""点を三角分割し、三角形を返します。
パラメータ
----------
pnts : array
形状または点群の配列形式のポイント。
geo.shapes または np.ndarray のいずれか。
out : array
三角形のポイント配列。
.. note::
シンプレックスは反時計回りに並んでいますが、この実装では逆になっています。
"""
pnts = np.unique(pnts, axis=0) # 一意なポイントのみ取得
avg = np.mean(pnts, axis=0)
p = pnts - avg
tri = Delaunay(p)
simps = tri.simplices
# -- インデックスホルダー、インデックスで埋め、最初を繰り返しCLへロール
# ポイントを元に戻す
z = np.zeros((len(simps), 4), dtype='int32')
z[:, :3] = simps
z[:, 3] = simps[:, 0]
tmp_ = p[z] + avg
new_pnts= []
for i in tmp_: # 時計回りに並べ替え
if _bit_area_(i) < 0.0: # -- 2025_10_27
i = i[::-1]
new_pnts.append(i)
return new_pnts
def _bit_area_(a):
"""ミニe_area、`areas`と`centroids`で使用されます。
負の面積は穴です。これは`shoelace`式とは意図的に逆になっています。
"""
a = _base_(a)
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 _area_centroid_(a):
r"""単一パート多角形`a`の面積と重心を計算します。
Geo配列のパーツの面積と重心計算にも使用されます。
注意事項
-----
複数パート形状の場合、この構文を使います:
>>> # 穴付き長方形の場合
>>> a0 = np.array([[[0., 0.], [0., 10.], [10., 10.], [10., 0.], [0., 0.]],
[[2., 2.], [8., 2.], [8., 8.], [2., 8.], [2., 2.]]])
>>> [_area_centroid_(i) for i in a0]
>>> [(100.0, array([ 5.00, 5.00])), (-36.0, array([ 5.00, 5.00]))]
"""
a = _base_(a)
x0, y1 = (a.T)[:, 1:]
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0)
e1 = np.einsum('...i,...i->...i', x1, y1)
t = e1 - e0
area = np.sum((e0 - e1) * 0.5)
x_c = np.sum((x1 + x0) * t, axis=0) / (area * 6.0)
y_c = np.sum((y1 + y0) * t, axis=0) / (area * 6.0)
return area, np.asarray([-x_c, -y_c])
こちらがCです。
C
array([[ 0.00, 0.00],
[ 0.00, 10.00],
[ 10.00, 10.00],
[ 10.00, 8.00],
[ 2.00, 8.00],
[ 2.00, 2.00],
[ 10.00, 2.00],
[ 10.00, 0.00],
[ 0.00, 0.00]])
上記で重心・面積・三角分割が求められます。
'どの三角形が元ジオメトリに属するか' を判定するには点内多角形検索を使います。これは以前の私の投稿でも扱いました。
'Point in Polygon ... Geometry Mysteries - Esri Community'
'winding number' アプローチを使っています。
'もしジオメトリオブジェクトの制約付きDelaunay三角分割が必要なら…試してみてください。'
'注意:'
'もっとpython、NumPy、ArcPy関連を見るなら私のgithubサイトへどうぞ。'
'Dan Patterson on github'
