Leia os tópicos de ajuda...
1 Featureclass Para Array NumPya0
FeatureClassToNumPyArray pode ser usado para extrair os vértices de feições geométricas. Se forem singlepart sem buracos, então você tem uma representação exata que pode ser usada para recriar as entradas. Se não for o caso, a reconstrução requer algum trabalho (veja https://github.com/Dan-Patterson/numpy_geometry para o módulo `npg`).
2 Array NumPy para Featureclassa0
NumPyArrayToFeatureClass não sobrescreverá uma feature class existente, mesmo se o ambiente overwriteOutput estiver definido como True.
------------------------------
Passo a passo... um notebook está anexado
# ---- Importações
import arcpy
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu
from numpy.lib.recfunctions import append_fields # como arcpy ExtendTable
Vou assumir que você está usando um notebook no ArcGIS Pro. Vou demonstrar uma conversão simples de polígonos para todos os pontos. As outras opções também podem ser codificadas, ou simplesmente compre um nível de licença mais alto do Pro.
a0
Assuma que sua camada está em um mapa. Especifique o nome da camada e se deseja carregar campos extras da tabela.
# Especifique a camada a usar para a variável `lyr_name`
# Isso funciona apenas para featureclasses... sem verificação adicional de erros... a pressão está alta
#
lyr_name = 'sq3' # -- insira o nome da featureclass
extra_fields = True # -- True mantém todos os campos extras, False caso contrário... apenas alguns? escreva o seu próprio
a0
Agora o grande passo.
- verifique se a camada existe no mapa ativo. Se existir, continue, caso contrário saia.
- o datasource da camada e a referência espacial são obtidos,
- uma lista de campos para descartar da saída (linhas 11-15)
- um array é criado (linhas 16-21). Este é um array numpy produzido pela função arcpy.da projetada para esse propósito.
- alguma limpeza segue, campos de saída são especificados e adicionados à saída (se houver)
- a featureclass de saída é criada (linha 36)
- e adicionada ao mapa (linha 37)
# ---- Projeto, mapa e coleção de camadas
# O projeto atual com o mapa ativo é verificado para o nome da camada que você especificou acima.
# Se existir, os vértices serão convertidos em uma featureclass de pontos.
#
aprx = arcpy.mp.ArcGISProject("CURRENT")
map_ = aprx.activeMap
lyrs = map_.listLayers(lyr_name)
if lyrs:
in_fc = lyrs[0].dataSource # -- caminho para a featureclass
SR = arcpy.Describe(in_fc).spatialReference
ignore_flds = ['OBJECTID', 'Shape', 'Shape_Area', 'Shape_Length'] # -- adicione se quiser
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
)
# -- verificações da featureclass de saída
out_fc = in_fc + "_pnts"
if arcpy.Exists(out_fc): # -- delete se existir
arcpy.management.Delete(out_fc)
#
# -- verificações do array e combinação dos valores das coordenadas
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("\nCampos de saída :\n{}".format(final_arr.dtype.names))
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] + 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.