Cenário de geometria de polígonos
Um monte de triângulos, alguns com buracos. Os números em verde indicam os valores de ID do anel externo do polígono. Os números em vermelho são os valores sequenciais de ID dos pontos dos anéis externo e interno.
Deve ser visualmente óbvio que há locais com valores de ID que:
- aparecem apenas uma vez (conexões de segmentos),
- têm dois números (início e fim do polígono),
- e outros têm 3, 4 ou até 5 IDs associados (onde os polígonos se encontram em diversas configurações).
Do arcpy ao numpy
Ao usar FeatureClassToNumPyArray para gerar a geometria no formato de array NumPy, há várias opções. Prefiro gerar o ID do objeto e os valores das coordenadas X e Y como uma saída separada. Se preciso trabalhar com atributos (texto ou outros valores numéricos), produzo uma saída separada, que pode ser facilmente associada novamente a qualquer saída reconstruída.
A sintaxe básica é:
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # necessário
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'] # -- os campos a usar
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Claro, você está trabalhando com alguns dados armazenados localmente para simplificar as coisas; então, modifique a linha 3 acima de acordo com seu caso de uso.
A linha 4 especifica os tokens de campo a serem usados para trazer o campo OBJECTID e os valores X e Y dos vértices do polígono.
As regras que sigo:
- Valido a classe de feição de entrada usando a validação de geometria OGC para garantir que os polígonos com buracos estejam ordenados de modo que os buracos (anéis internos) venham depois de seu anel externo correspondente. Não consigo imaginar um caso de uso em que a geometria de um buraco venha antes de seu anel externo. Isso evita ter que associar anéis de buracos ao anel externo por meio de verificações de ponto no polígono.
- Dados projetados e uma projeção que seja, no mínimo, conforme (UTM, MTM) ou, em alguns casos, de área equivalente. Definitivamente, não aquelas projeções web, que são inúteis para determinar atributos de geometria.
- Os anéis externos são orientados no sentido horário; os buracos (anéis internos), no sentido anti-horário.
- Os pontos inicial e final de um polígono têm valores X e Y duplicados, mas um número de ID diferente.
- além de todas as outras coisas especificadas pela ESRI.
Leia Verificar geometria e Reparar geometria para obter informações relevantes.
Saída básica do array
Agora, de volta à geração do array.
# ---- continuando do trecho acima
ids = out['OID@'].astype('int') # -- obtém um ndarray simples do NumPy com os valores
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- converte o array estruturado em ndarray
LL = np.min(xys, axis=0)
arr = xys - LL
#
ids
[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]
# ... trecho omitido
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
A linha 2 separa os valores de objectid do array estruturado derivado, convertendo-os em inteiros simples.
A linha 3 realiza uma operação semelhante, trazendo as coordenadas X e Y como um ndarray double do NumPy (stué uma importação do NumPy com renomeação no primeiro bloco de código)
A linha 4 desloca as coordenadas para um plano menos complicado. Subtraio os valores mínimos de X e Y de todos os valores de coordenadas. No meu exemplo, a coordenada (0.0, 0.0) é, na verdade, (300,000.00, 5,000,000.00). O plano XY foi deslocado, facilitando a leitura e simplificando os cálculos. Ao terminar, basta deslocar os resultados dos cálculos de volta para as coordenadas originais, então não há problema.
As linhas 7 e 8 mostram os valores de ID da primeira forma (1), que tem 9 pontos formando o anel externo e seu buraco, e os 4 pontos associados à última forma (5). As linhas 10 a 24 contêm as coordenadas dos pontos. As linhas 11 e 14 são as coordenadas inicial e final do anel externo (triângulo), e as linhas 15 a 19 são os pontos do buraco quadrado.
O que é possível fazer agora
Neste ponto, você tem vários pontos que representam a geometria dos polígonos. Você pode
- deslocá-los (basta somar/subtrair valores às/das coordenadas X e Y)
- aumentar o número de pontos ao longo dos anéis ou simplificá-los (há funções do NumPy para isso
- etc.
Mas, se precisar trabalhar com polígonos individuais, terá que separá-los. É aqui que há várias opções para isolar os próprios polígonos e/ou seus anéis internos e externos.
def array_to_rings(arr, from_to_ids=False):
"""Particiona um array em anéis externos e internos.
Parâmetros
----------
arr : array_like
Um array/lista/array geográfico Nx2 de coordenadas x, y, que pode conter um
ou mais anéis externos e internos (buracos).
Retorna
-------
Uma lista de arrays, possivelmente aninhados se houver anéis internos/buracos.
"""
# --
def _area_(a):
"""Versão modificada de `_bit_area_`. Áreas negativas são buracos. Invertida em relação
à fórmula do `shoelace`.
"""
x0, y1 = (a.T)[:, 1:] # produto vetorial organizado da seguinte forma
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # modificado em 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Verificação de sentido anti-horário. Invertida em relação à fórmula do shoelace."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nÉ necessário um np.ndarray Nx2."
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 # -- contador de anéis externos
outer_ = []
inner_ = []
frto_ = [] # -- IDs de início e fim para dividir o array de entrada `arr`
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nNenhum ponto correspondente; o array de entrada está com defeito.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # tamanho do subarray
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # adiciona os anéis de buracos a um anel externo existente
out[-1].append(sub)
else: # apenas para o caso de um buraco ser encontrado antes de um anel externo
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- retorna uma lista de arrays ou uma lista de listas de arrays, de alguma forma
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
E a saída:
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]])],
[... trecho omitido],
[... trecho omitido],
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]])]
A saída acima mostra que a primeira forma é uma lista de dois arrays: o anel externo e o buraco. Os dois últimos arrays correspondem às formas 4 e 5, que são triângulos simples.
Neste ponto, poderíamos gerar uma saída que removesse os buracos.
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]])]
Em seguida, você pode converter os arrays em polígonos do arcpy.
def array_poly(arrs, p_type=None, sr=None):
"""Monta feições poligonais a partir de arrays geográficos ou ndarrays.
Usada por `geometry_fc` ou pode ser usada separadamente.
Parâmetros
----------
arr : arrays
Array(s) de pontos.
p_type : texto
POLYGON ou POLYLINE
sr : referência espacial
Objeto, nome ou ID de referência espacial.
Observações
-----
É possível criar feições de linha ou polígono a partir dos dados do array. As
feições podem ser multipartes, com ou sem anéis internos.
Os anéis externos são ordenados no sentido horário; os anéis internos (buracos), no sentido
anti-horário.
"""
def _arr_poly_(arr, SR, as_type):
"""Divide o array nos pontos em que aparecem valores nan, separando-os."""
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
E depois, usando o código acima, em uma classe de feição.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Forma feições poligonais a partir da lista de arrays criada por `fc_geometry`.
Parâmetros
----------
a : array ou lista de arrays
Alguns podem ser arrays de objetos, normalmente criados por ``pnts_arr``
p_type : string
Tipo de geometria em maiúsculas, por exemplo, POLYGON.
gdb : texto
Caminho e nome da geodatabase.
fname: texto
Nome da classe de feições.
sr: referência espacial
nome ou objeto
Retorna
-------
Classes de feições de parte única e/ou multipartes.
Observações
-----
O trabalho é feito por ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # chama array_poly
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # o legado é 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
O que equivale à ferramenta Eliminar parte do polígono que, curiosamente, exige uma licença Advanced do Pro. Por que você precisaria de uma licença avançada para remover/preencher um buraco?
Ao alterar o tipo de saída acima para POLYLINE, você transforma polígonos em linhas, assim como a ferramenta Polígono para linha, que também exige uma licença Advanced.
Este é apenas um exemplo simples de como você pode usar FeatureClassToNumPyArray para trabalhar com geometria.
Mais ferramentas e código disponíveis em Geometria NumPy e Ferramentas para ArcGIS Pro.
