ヘルプトピックを読む...
1 Featureclass を NumPy 配列に変換
FeatureClassToNumPyArray はジオメトリフィーチャから頂点を抽出するために使用できます。もしそれらが穴のないシングルパートであれば、入力を再作成するために使用できる正確な表現があります。そうでない場合は再構築に多少の作業が必要です（詳細は https://github.com/Dan-Patterson/numpy_geometry の `npg` モジュールを参照）。
2 NumPy 配列を Featureclass に変換
NumPyArrayToFeatureClass は、overwriteOutput 環境が True に設定されていても既存のフィーチャクラスを上書きしません。
------------------------------
段階的に... ノートブックが添付されています
# ---- インポート
import arcpy
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu
from numpy.lib.recfunctions import append_fields # arcpy の ExtendTable のようなもの
ArcGIS Pro のノートブックを使用していると仮定します。ポリゴンをすべてのポイントに単純に変換するデモを行います。他のオプションもコード化可能ですし、単により高いライセンスレベルの Pro を購入することもできます。
レイヤーがマップ内にあると仮定します。レイヤー名と追加のテーブルフィールドを引き継ぐかどうかを指定してください。
# `lyr_name` 変数に使用するレイヤーを指定
# これはフィーチャクラスのみで動作します... さらなるエラーチェックはありません... プレッシャーがかかっています
#
lyr_name = 'sq3' # -- フィーチャクラス名を入力してください
extra_fields = True # -- True はすべての追加フィールドを保持、False はそうしません... 一部だけ？自分で書いてください
さて、大きなステップです。
- アクティブマップにレイヤーが存在するか確認します。存在すれば続行し、そうでなければ中止します。
- レイヤーデータソースと空間参照が取得され、
- 出力から除外するフィールドのリスト（11-15行目）
- 配列が作成されます（16-21行目）。これはこの目的のために設計された arcpy.da 関数によって生成された numpy 配列です。
- いくつかのクリーンアップが続き、出力フィールドが指定され、出力に追加されます（もしあれば）
- 出力フィーチャクラスが作成されます（36行目）
- そしてマップに追加されます（37行目）
# ---- プロジェクト、マップ、およびレイヤーコレクション
# 現在のプロジェクトとアクティブマップで上記で指定したレイヤー名がチェックされます。
# 存在すれば、頂点はポイントフィーチャクラスに変換されます。
#
aprx = arcpy.mp.ArcGISProject("CURRENT")
map_ = aprx.activeMap
lyrs = map_.listLayers(lyr_name)
if lyrs:
in_fc = lyrs[0].dataSource # -- フィーチャクラスへのパス
SR = arcpy.Describe(in_fc).spatialReference
ignore_flds = ['OBJECTID', 'Shape', 'Shape_Area', 'Shape_Length'] # -- 必要なら追加してください
keep_flds = ["OID@", "SHAPE@X", "SHAPE@Y"]
if extra_fields:
add_flds = [f.name for f in arcpy.ListFields(in_fc) if f.name not in ignore_flds]
keep_flds = keep_flds + add_flds
arr = arcpy.da.FeatureClassToNumPyArray(
in_fc,
field_names=keep_flds,
spatial_reference=SR,
explode_to_points=True
)
# -- 出力フィーチャクラスのチェック
out_fc = in_fc + "_pnts"
if arcpy.Exists(out_fc): # -- 存在すれば削除する
arcpy.management.Delete(out_fc)
#
# -- 配列のチェックと座標値の結合処理
dt_ = arr.dtype
dt_names = [(i[0].replace("@", "_"), i[1]) for i in arr.dtype.descr]
arr.dtype = dt_names
#
final_arr = append_fields(arr, ["Pnt_x", "Pnt_y"], [arr["SHAPE_X"], arr["SHAPE_Y"]])
print("\n出力フィールド :\n{}".format(final_arr.dtype.names))
array([(1, 300009. , 5000001. ), (1, 300000. , 5000000. ),
(1, 300002. , 5000008. ), (1, 300008. , 5000010. ),
... snip
(7, 300004. , 5000009. ), (7, 300002. , 5000009. ),
(7, 300002. , 5000011. )],
dtype=[('OID@', '<i4'), ('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8')])
You can then split the array into its constituent shapes using the OID@ field. Basically determine the difference in the sequential id values, where they differ from 0, get the positional index to use for splitting.
spl = np.nonzero(np.diff(arr["OID@"]))[0] + 1
From there, you can produce subarrays
subs = np.array_split(arr, spl)
subs[0] # -- the first as an example
array([(1, 300009.000, 5000001.000), (1, 300000.000, 5000000.000),
(1, 300002.000, 5000008.000), (1, 300008.000, 5000010.000),
(1, 300010.000, 5000010.000), (1, 300010.000, 5000008.000),
(1, 300009.000, 5000001.000), (1, 300003.000, 5000003.000),
(1, 300007.000, 5000003.000), (1, 300005.000, 5000007.000),
(1, 300003.000, 5000003.000)],
dtype=[('OID@', '<i4'), ('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8')])
From there you could process to just take the first or last point (one of the main tool options that requires and advanced licence).
frst = [i[0] for i in subs] # let's grab the first point of each shape
frst
[(1, 300009.000, 5000001.000),
(2, 300010.000, 5000008.000),
(3, 300008.000, 5000011.000),
(4, 300006.000, 5000012.000),
(5, 300002.500, 5000013.000),
(6, 300004.000, 5000012.500),
(7, 300002.000, 5000011.000)]
As a teaser, you can even place a point at a certain percentage along the polygon perimeter with a little numpy magic
from numpy.lib.recfunctions import structured_to_unstructured as stu
base = [stu(a[["SHAPE@X", "SHAPE@Y"]]) for a in subs]
halfway = [_percent_along_(a, percent=50) for a in base]
halfway
[array([ 300010.000, 5000009.950]),
array([ 300008.500, 5000012.000]),
array([ 300005.000, 5000010.310]),
array([ 300005.056, 5000014.972]),
array([ 300000.214, 5000010.643]),
array([ 300002.527, 5000012.991]),
array([ 300003.894, 5000009.000])]
and the function for placing the points
def _percent_along_(a, percent=0):
"""Add a point along a poly feature at a distance from the start point."""
# a = _base_(a)
if percent > 1.:
percent /= 100.
dxdy = a[1:, :] - a[:-1, :] # coordinate differences
leng = np.sqrt(np.einsum('ij,ij->i', dxdy, dxdy)) # segment lengths
cumleng = np.concatenate(([0], np.cumsum(leng)))
perleng = cumleng / cumleng[-1]
if percent <= 0: # check for faulty distance or start point
return a[0]
if percent >= perleng[-1]: # check for greater distance than cumulative
return a[-1]
_end_ = np.digitize(percent, perleng)
x1, y1 = a[_end_]
_start_ = _end_ - 1
x0, y0 = a[_start_]
t = percent - perleng[_start_]
xt = x0 * (1. - t) + (x1 * t)
yt = y0 * (1. - t) + (y1 * t)
return np.array([xt, yt])
Or you can 'arcpy' it if you prefer to crank up Pro.
So, basic functionality shouldn't be completely removed by license levels. If you have the need convert features to points, remember that tools are available within.
