Escenario de geometría de polígonos
Un montón de triángulos, algunos con huecos. Los números verdes indican los valores de ID del anillo exterior del polígono. Los números rojos son los valores secuenciales de ID de los puntos de los anillos exteriores e interiores.
Debería ser evidente a simple vista que hay ubicaciones con valores de ID que:
- aparecen una sola vez (conexiones de segmentos),
- tienen dos números (inicio y fin del polígono),
- y otras que tienen 3, 4 o incluso 5 ID asociados (donde los polígonos se encuentran con diversas configuraciones).
De arcpy a numpy
Al usar FeatureClassToNumPyArray para exportar la geometría al formato de matriz de NumPy, tienes varias opciones. Prefiero exportar el ID del objeto y los valores de las coordenadas X e Y por separado. Si necesito trabajar con atributos (texto u otros valores numéricos), genero una salida independiente que se puede volver a unir fácilmente a cualquier salida reconstruida.
La sintaxis básica es:
importar numpy como np
de numpy.lib.recfunctions importar structured_to_unstructured como stu # necesario
importar arcpy
de arcpy.da importar FeatureClassToNumPyArray
in_fc = r"C:\path_to_gdb\the_geodatabase.gdb\featureclass_name"
flds = ['OID@', 'SHAPE@X', 'SHAPE@Y'] # -- los campos que se usarán
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Por supuesto, estás trabajando con datos almacenados localmente para simplificar las cosas, así que modifica la línea 3 anterior según tu caso.
La línea 4 especifica los tokens de campo que se usarán para incorporar el campo OBJECTID y los valores X e Y de los vértices del polígono.
Las reglas que sigo:
- Valido la clase de entidad de entrada con la validación de geometría OGC para asegurarme de que los polígonos con huecos estén ordenados de modo que los huecos (anillos interiores) sigan a su anillo exterior correspondiente. No se me ocurre ningún caso en el que la geometría de un hueco aparezca antes que su anillo exterior. Así no hay que asignar los anillos de huecos a su anillo exterior mediante comprobaciones de punto en polígono.
- Datos proyectados y una proyección que sea al menos conforme (UTM, MTM) o, en algunos casos, de áreas equivalentes. Desde luego, no esas proyecciones web, que son inútiles para determinar atributos geométricos.
- Los anillos exteriores están orientados en sentido horario; los huecos (anillos interiores), en sentido antihorario.
- El punto inicial y final de un polígono tienen valores X e Y duplicados, pero un número de ID distinto.
- además de todo lo demás que especifica ESRI.
Lee Comprobar geometría y Reparar geometría para obtener información pertinente.
Salida básica de matriz
Volvamos ahora a la generación de la matriz.
# ---- continuando desde arriba
ids = out['OID@'].astype('int') # -- obtener una matriz numpy sencilla de valores
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- convertir la matriz estructurada en una matriz ndarray
LL = np.min(xys, axis=0)
arr = xys - LL
#
ids
[1 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]
# ... fragmento omitido
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
La línea 2 separa los valores objectid de la matriz estructurada derivada y los convierte en enteros simples.
La línea 3 realiza una operación similar al extraer las coordenadas X e Y como una matriz ndarray double de numpy (stues una importación de numpy con un nombre abreviado en el primer bloque de código)
La línea 4 desplaza las coordenadas a un plano menos engorroso. Resto los valores mínimos de X e Y de todos los valores de coordenadas. En mi ejemplo, la coordenada (0.0, 0.0) corresponde en realidad a (300,000.00, 5,000,000.00). El plano XY se ha desplazado, lo que facilita la lectura y simplifica los cálculos. Cuando termines, vuelve a desplazar los cálculos a sus coordenadas originales, así que no hay ningún problema.
Las líneas 7 y 8 contienen los valores de ID de la primera forma (1), que tiene 9 puntos que forman el anillo exterior y su hueco, y los 4 puntos asociados a la última forma (5). Las líneas 10-24 contienen las coordenadas de los puntos. Las líneas 11 y 14 contienen las coordenadas inicial y final del anillo exterior (triángulo), y las líneas 15-19, los puntos del hueco cuadrado.
Qué se puede hacer ahora
En este punto tienes un conjunto de puntos que representan la geometría de los polígonos. Podrías
- desplazarlos (simplemente sumando o restando valores a las coordenadas X e Y)
- aumentar la densidad de puntos a lo largo de los anillos o simplificarlos (hay funciones de numpy para esto
- y demás
Pero si tienes que trabajar con polígonos individuales, tendrás que separarlos. Aquí hay varias opciones para aislar los polígonos en sí o sus anillos interiores y exteriores.
def array_to_rings(arr, from_to_ids=False):
"""Divide una matriz en anillos exteriores e interiores.
Parámetros
----------
arr : array_like
Una matriz/lista/matriz geográfica Nx2 de coordenadas x e y, que puede contener uno
o más anillos exteriores e interiores (huecos).
Devuelve
-------
Una lista de matrices, posiblemente anidadas si hay anillos interiores/huecos.
"""
# --
def _area_(a):
"""Versión modificada de `_bit_area_`. Las áreas negativas son huecos. Se invierte
respecto a la fórmula del `cordón`.
"""
x0, y1 = (a.T)[:, 1:] # producto cruzado dispuesto así
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # modificado el 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Comprobación de sentido antihorario. Se invierte la fórmula del cordón."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nSe requiere una 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 anillos exteriores
outer_ = []
inner_ = []
frto_ = [] # -- ID de inicio y fin para dividir la matriz 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("\nNo hay ningún punto coincidente; la matriz de entrada no es válida.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # tamaño de la submatriz
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # añadir anillos de huecos a un anillo exterior existente
out[-1].append(sub)
else: # por si acaso se encuentra un hueco antes que un anillo exterior
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- devuelve una lista de matrices o una lista de listas de matrices, según corresponda
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
Y el resultado:
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]])],
[... fragmento omitido],
[... fragmento 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]])]
El resultado anterior muestra que la primera forma es una lista de dos matrices: el anillo exterior y el hueco. Las dos últimas matrices corresponden a las formas 4 y 5, que son triángulos simples.
En este punto podríamos generar una salida que elimine los huecos.
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]])]
Luego puedes convertir las matrices en polígonos de arcpy.
def array_poly(arrs, p_type=None, sr=None):
"""Construye entidades poligonales a partir de matrices geográficas o ndarray.
Se usa con `geometry_fc` o por separado.
Parámetros
----------
arr : matrices
Matriz/matrices de puntos.
p_type : texto
POLYGON o POLYLINE
sr : referencia espacial
Objeto, nombre o ID de referencia espacial.
Notas
-----
Se pueden crear entidades de polilínea o polígono a partir de los datos de la matriz. Las
entidades pueden ser multiparte, con anillos interiores o sin ellos.
Los anillos exteriores están ordenados en sentido horario; los anillos interiores (huecos),
en sentido antihorario.
"""
def _arr_poly_(arr, SR, as_type):
"""Divide la matriz donde aparecen valores nan, separándolos."""
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
Y después, convertirlo en una clase de entidad usando lo anterior.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Crea entidades poligonales a partir de la lista de matrices creada por `fc_geometry`.
Parámetros
----------
a : matriz o lista de matrices
Algunas pueden ser matrices de objetos, normalmente creadas por ``pnts_arr``
p_type : cadena
Tipo de geometría en mayúsculas, por ejemplo, POLYGON.
gdb : texto
Ruta y nombre de la geodatabase.
fname : texto
Nombre de la clase de entidad.
sr : referencia espacial
nombre u objeto
Devuelve
-------
Clases de entidad de una sola parte o de varias partes.
Notas
-----
El trabajo lo realiza ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # llamar a array_poly
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # el nombre heredado es 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
Lo que equivale a la herramienta Eliminar parte de polígono que, curiosamente, requiere una licencia Advanced de Pro. ¿Por qué se necesitaría una licencia avanzada para eliminar o rellenar un hueco?
Al cambiar el tipo de salida anterior a POLYLINE, se convierten los polígonos en polilíneas, igual que con Polígono a línea que también requiere una licencia Advanced.
Este es solo un ejemplo sencillo de cómo usar FeatureClassToNumPyArray para trabajar con geometría.
Hay más herramientas y código disponibles en Geometría NumPy y Herramientas para ArcGIS Pro.
