Feature vertices to points ... for Basic ArcGIS Pro
Read the help topics...
FeatureClassToNumPyArray can be used to pull the vertices out of geometry features. If they are singlepart without holes, then you have an exact representation that can be used to recreate the inputs. If that isn't the case then reconstruction takes some work (see https://github.com/Dan-Patterson/numpy_geometry for the `npg` module).
NumPyArrayToFeatureClass will not overwrite an existing feature class, even if the overwriteOutput environment is set to True.
------------------------------
Stepwise... a notebook is attached
# ---- Imports
import arcpy
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu
from numpy.lib.recfunctions import append_fields # like arcpy ExtendTableI will assume that you are using a notebook in ArcGIS Pro. I will demo a simple conversion of polygons to all points. The other options can also be coded for, or simply buy a higher license level of Pro.
Assume that your layer is in a map. Specify the layer name and whether you want to carry over extra table fields.
# Specify the layer to use for the `lyr_name` variable
# This only works for featureclasses... no further error checking... the pressure is on
#
lyr_name = 'sq3' # -- input the featureclass name
extra_fields = True # -- True keeps all extra fields, False otherwise... just some? write your own
Now the big step.
- check that the layer exists in the active map. If it does, carry on, otherwise bail.
- the layer datasource and the spatial reference is obtained,
- a list of fields to dump from the output (lines 11-15)
- an array is created (lines 16-21). This is a numpy array produced the arcpy.da function designed for this purpose.
- some cleanup follows, output fields are specified and appended to the output (if any)
- the output featureclass is created (line 36)
- and added to the map (line 37)
# ---- Project, map and layer collection
# The current project with the active map is checked for the layer name you specified above.
# If it exists, the vertices will be converted to a point featureclass.
#
aprx = arcpy.mp.ArcGISProject("CURRENT")
map_ = aprx.activeMap
lyrs = map_.listLayers(lyr_name)
if lyrs:
in_fc = lyrs[0].dataSource # -- path to the featureclass
SR = arcpy.Describe(in_fc).spatialReference
ignore_flds = ['OBJECTID', 'Shape', 'Shape_Area', 'Shape_Length'] # -- add if you want
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
)
# -- output featureclass checks
out_fc = in_fc + "_pnts"
if arcpy.Exists(out_fc): # -- delete if it exists
arcpy.management.Delete(out_fc)
#
# -- array checks and combining coordinate values
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("\nOutput fields :\n{}".format(final_arr.dtype.names))
#
# -- produce the output and add to the map.
arcpy.da.NumPyArrayToFeatureClass(final_arr, out_fc, ["SHAPE_X", "SHAPE_Y"], spatial_reference=SR)
map_.addDataFromPath(out_fc) # -- add results to the map
#
else:
print("No layer called {}".format(lyr_name))And the results are
I will spare you the table contents, but they are there.
--------------
If you prefer the python window
import arcpy
import numpy as np
in_fc = r"C:\...path to your layer...\...layer_name..."
keep_flds = ["OID@", "SHAPE@X", "SHAPE@Y"]
SR = arcpy.Describe(in_fc).spatialReference
arr = arcpy.da.FeatureClassToNumPyArray(
in_fc,
field_names=keep_flds,
spatial_reference=SR,
explode_to_points=True
)Which shows the structured array data format. (snipped to simplify)
arr
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] + 1From 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.
