Scénario de géométrie polygonale
Un ensemble de triangles, dont certains comportent des trous. Les chiffres verts indiquent les identifiants des anneaux extérieurs des polygones. Les chiffres rouges indiquent les identifiants séquentiels des points des anneaux extérieurs et intérieurs.
Il devrait être évident à l’œil qu’il existe des emplacements dont les identifiants :
- n’apparaissent qu’une seule fois (jonctions de segments),
- apparaissent deux fois (début et fin du polygone),
- et d’autres sont associés à 3, 4, voire 5 identifiants (là où des polygones se rejoignent selon différentes configurations).
D’arcpy à numpy
Lorsque vous utilisez FeatureClassToNumPyArray pour convertir la géométrie au format de tableau NumPy, plusieurs options s’offrent à vous. Je préfère extraire séparément l’identifiant de l’objet et les valeurs des coordonnées X et Y. Si je dois travailler avec des attributs (texte ou autres valeurs numériques), je produis un résultat distinct, qui peut facilement être joint à tout résultat reconstruit.
La syntaxe de base est la suivante :
import numpy as np
from numpy.lib.recfunctions import structured_to_unstructured as stu # nécessaire
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'] # -- les champs à utiliser
out = FeatureClassToNumPyArray(in_fc, flds, explode_to_points=True)
Bien sûr, vous travaillez avec des données stockées localement pour simplifier les choses ; modifiez donc la ligne 3 ci-dessus selon votre cas d’utilisation.
La ligne 4 indique les jetons de champ à utiliser pour récupérer le champ OBJECTID ainsi que les valeurs X et Y des sommets du polygone.
Les règles que j’applique :
- Je valide la classe d’entités en entrée à l’aide de la validation géométrique OGC afin de m’assurer que les polygones comportant des trous sont ordonnés de sorte que les trous (anneaux intérieurs) suivent leur anneau extérieur parent. Je ne vois vraiment aucun cas d’utilisation où la géométrie d’un trou précéderait son anneau extérieur. Cela évite d’avoir à associer les anneaux des trous à leur anneau extérieur à l’aide de tests point dans polygone.
- Des données projetées et une projection au moins conforme (UTM, MTM) ou, dans certains cas, équivalente. Surtout pas ces projections Web inutiles pour déterminer les attributs géométriques.
- Les anneaux extérieurs sont orientés dans le sens horaire, les trous (anneaux intérieurs) dans le sens antihoraire.
- Le point de début et le point de fin d’un polygone ont les mêmes valeurs X et Y, mais des numéros d’identification différents.
- et tout le reste précisé par ESRI.
Consultez Vérifier la géométrie et Réparer la géométrie pour obtenir des renseignements pertinents.
Résultat de tableau de base
Revenons maintenant à la création du tableau.
# ---- suite du code ci-dessus
ids = out['OID@'].astype('int') # -- obtenir un simple ndarray NumPy de valeurs
xys = stu(out[['SHAPE@X', 'SHAPE@Y']]) # -- convertir le tableau structuré en 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]
# ... extrait omis
[ 40.0, 0.0]
[ 35.0, 10.0]
[ 40.0, 20.0]
[ 40.0, 0.0]]
La ligne 2 extrait les valeurs objectid du tableau structuré obtenu et les convertit en entiers simples.
La ligne 3 effectue une opération similaire en récupérant les coordonnées X et Y sous la forme d’un ndarray NumPy de nombres réels double précision (stuest un import NumPy renommé dans le premier bloc de code)
La ligne 4 déplace les coordonnées vers un plan moins contraignant. Je soustrais les valeurs X et Y minimales de toutes les coordonnées. Dans mon exemple, la coordonnée (0.0, 0.0) correspond en réalité à (300,000.00, 5,000,000.00). Le plan XY a été décalé, ce qui facilite la lecture et simplifie les calculs. Une fois terminé, il suffit de rétablir les coordonnées d’origine dans les calculs ; cela ne pose donc aucun problème.
Les lignes 7 et 8 contiennent les valeurs d’identifiant de la première forme (1), qui compte 9 points formant l’anneau extérieur et son trou, ainsi que les 4 points associés à la dernière forme (5). Les lignes 10 à 24 contiennent les coordonnées des points. Les lignes 11 et 14 indiquent les coordonnées de début et de fin de l’anneau extérieur (triangle), et les lignes 15 à 19 celles des points du trou carré.
Que peut-on faire maintenant ?
À ce stade, vous disposez d’un ensemble de points représentant la géométrie des polygones. Vous pourriez
- les déplacer (il suffit d’ajouter ou de soustraire des valeurs aux coordonnées X et Y)
- densifier le nombre de points le long des anneaux ou les simplifier (il existe des fonctions NumPy pour cela
- etc.
Mais si vous devez travailler avec des polygones individuels, il faut les séparer. Plusieurs options permettent alors d’isoler les polygones eux-mêmes et/ou leurs anneaux intérieurs et extérieurs.
def array_to_rings(arr, from_to_ids=False):
"""Répartit un tableau entre anneaux extérieurs et intérieurs.
Paramètres
----------
arr : array_like
Tableau/liste/tableau géographique Nx2 de coordonnées x, y, pouvant contenir un
ou plusieurs anneaux extérieurs et intérieurs (trous).
Valeurs renvoyées
----------------
Liste de tableaux, éventuellement imbriqués en présence d’anneaux intérieurs/trous.
"""
# --
def _area_(a):
"""Version modifiée de `_bit_area_`. Les aires négatives correspondent aux trous. Formule inversée par rapport à
la formule du « lacet ».
"""
x0, y1 = (a.T)[:, 1:] # produit vectoriel défini comme suit
x1, y0 = (a.T)[:, :-1]
e0 = np.einsum('...i,...i->...i', x0, y0) # modifié le 2024-03-28
e1 = np.einsum('...i,...i->...i', x1, y1)
return np.sum((e0 - e1) * 0.5)
def _is_ccw_(a):
"""Vérifie le sens antihoraire. Formule inversée par rapport à celle du lacet."""
return 0 if _area_(a) > 0. else 1
# --
#
msg = "\nUn np.ndarray Nx2 est requis."
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 # -- compteur d’anneaux extérieurs
outer_ = []
inner_ = []
frto_ = [] # -- identifiants de début et de fin pour séparer le tableau d’entrée `arr`
while n < N:
ids = np.nonzero((arr[fst] == arr[fst:]).all(-1))[0]
if len(ids) >= 2:
pair = ids[:2]
else:
print("\nAucun point correspondant ; le tableau d’entrée est défectueux.")
return None
fst, lst = pair + n
sub = arr[fst:lst + 1]
n = n + len(sub) # taille du sous-tableau
res = _is_ccw_(sub)
if res != 1:
outer_.append(sub)
out.append([sub])
c_o += 1
else:
inner_.append(sub)
if len(out) > 0: # ajouter les anneaux des trous à un anneau extérieur existant
out[-1].append(sub)
else: # au cas où un trou apparaîtrait avant un anneau extérieur
out.append([sub])
frto_.append([c_o, fst, lst + 1])
fst = lst + 1
# -- renvoyer une liste de tableaux ou, sous une forme ou une autre, une liste de listes de tableaux
out = [i[0] if len(i) == 1 else i for i in out]
if from_to_ids:
return out, frto_
return out
Et voici le résultat :
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]])],
[... extrait omis],
[... extrait omis],
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]])]
Le résultat ci-dessus montre que la première forme est une liste de deux tableaux : l’anneau extérieur et le trou. Les deux derniers tableaux correspondent aux formes 4 et 5, qui sont de simples triangles.
Nous pourrions à ce stade produire un résultat sans les trous.
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]])]
Vous pouvez ensuite convertir les tableaux en polygones arcpy.
def array_poly(arrs, p_type=None, sr=None):
"""Assemble des entités polygonales à partir de tableaux Geo ou ndarray.
Utilisée par `geometry_fc` ou pouvant être utilisée séparément.
Paramètres
----------
arr : tableaux
Tableau(x) de points.
p_type : texte
POLYGON ou POLYLINE
sr : référence spatiale
Objet, nom ou identifiant de référence spatiale.
Remarques
---------
Des entités polylignes ou polygones peuvent être créées à partir des données du tableau. Les
entités peuvent être multiparties, avec ou sans anneaux intérieurs.
Les anneaux extérieurs sont ordonnés dans le sens horaire ; les anneaux intérieurs (trous) sont ordonnés
dans le sens antihoraire.
"""
def _arr_poly_(arr, SR, as_type):
"""Découpe le tableau aux endroits où apparaissent des valeurs nan, en les séparant."""
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
Puis créer une classe d’entités à l’aide de la fonction ci-dessus.
def geometry_fc(a, p_type=None, gdb=None, fname=None, sr=None):
"""Crée des entités polygonales à partir de la liste de tableaux créée par `fc_geometry`.
Paramètres
----------
a : tableau ou liste de tableaux
Certains peuvent être des tableaux objet, généralement créés par ``pnts_arr``
p_type : chaîne
Type de géométrie en majuscules, par exemple POLYGON.
gdb : texte
Chemin et nom de la géodatabase.
fname : texte
Nom de la classe d’entités.
sr : référence spatiale
nom ou objet
Renvoie
-------
Classes d’entités à partie unique et/ou à parties multiples.
Remarques
-----
Le traitement est effectué par ``array_poly``.
"""
if p_type is None:
p_type = "POLYGON"
out = array_poly(a, p_type.upper(), sr=sr, IFT=IFT) # appeler array_poly
name = gdb + "/" + fname
wkspace = arcpy.env.workspace = r'memory' # l’ancien format est 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
Ce qui équivaut à l’outil Eliminate Polygon Part qui, curieusement, nécessite une licence Advanced dans Pro. Pourquoi donc faudrait-il une licence Advanced pour supprimer ou remplir un trou ?
En remplaçant le type de sortie ci-dessus par POLYLINE, vous convertissez les polygones en polylignes, comme avec Polygon to Line qui nécessite également une licence Advanced.
Voici un exemple simple d’utilisation de FeatureClassToNumPyArray pour travailler avec la géométrie.
Vous trouverez d’autres outils et du code sur NumPy Geometry et Tools for ArcGIS Pro.
