#1 Fri 20 July 2012 17:22
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
[python] Polyligne 3D vers polygone 3D
Bonjour à tous,
Mon problème est le suivant : j'ai une polyligne 3D (deux points uniquement) que je souhaite extruder jusqu’à un Zmin renseigné afin de créer un polygone 3D. Actuellement je bloque car dès que je compte le nombre de pièces cela me retourne 0 (voir le résultat ci dessous). J'aurais donc besoin de votre aide pour comprendre d'ou vient mon erreur.
SCRIPT :
Code:
import arcpy rows=arcpy.UpdateCursor("W:/Florian/Contour.shp", "", "", "ZMIN") featureList = [] for row in rows: polyligne = row.shape list_pnt=polyligne.getPart(0) print str(len(list_pnt)) p1 = arcpy.Point() p2 = arcpy.Point() p1s = arcpy.Point() p2s = arcpy.Point() array = arcpy.Array() i = 0 for coordPair in list_pnt: if i==0: i=1 p1.X = coordPair.X p1.Y = coordPair.Y p1.Z = coordPair.Z print p1 else: p2.X = coordPair.X p2.Y = coordPair.Y p2.Z = coordPair.Z print p2 p1s=p1 p1s.Z=row.ZMIN print p1s p2s=p2 p2s.Z=row.ZMIN print p2s array.add(p1) array.add(p1s) array.add(p2s) array.add(p2) array.add(p1) # Create a Polygon object based on the array of points # polygon = arcpy.Polygon(array,None,True) print polygon.partCount print polygon.type # Clear the array for future use # array.removeAll() # Delete cursor and row objects to remove locks on the data # del row del rows
RESULTAT :
Code:
2 147591,5449 6838917,3373 83,4290000000037 NaN 147540,006999999 6838905,8854 81,8790000000008 NaN 147591,5449 6838917,3373 76,6770820172 NaN 147540,006999999 6838905,8854 76,6770820172 NaN 0 polygon
Merci.
Dernière modification par Floflo49fb (Fri 20 July 2012 17:23)
Florian Boret
Dream it, Make it, Share it
Hors ligne
#2 Sat 21 July 2012 12:50
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Bonjour,
Je n'ai pas la réponse direct au problème mais j'ai quelques remarques sur le code. N'ayant pas ArcPy, je n'ai rien vérifié dont il se peut que ces remarques ne soient pas toujours justifiées !
D'abord je pense que toute la partie du code dédiée à la récupération de p1 et p2 pourrait fortement être simplifiée par quelque chose dans ce genre :
Code:
for row in rows: polyligne = row.shape p1, p2 = polyligne.getPart(0)
Éventuellement si l'on veut se prévenir d'une exception sur le nombre de points :
Code:
for row in rows: polyligne = row.shape list_pnt = polyligne.getPart(0) try: p1, p2 = list_pnt except: print "Erreur nombre de points supérieur à 2" break #casse la boucle, on passe au feature suivant
Idéalement, je bouclerai aussi sur les différentes parties des entités (à moins que l'on soit absolument certain qu'il n'y aura jamais d'objet multipart). ça n'alourdira pas vraiment le code et ce sera toujours plus propre que de laisser getpart(0).
Globalement l'utilisation de i pr tester si l'on est au 1er ou au second passage de boucle sur list_pnt n'est pas très élégante. Et la partie du code permettant de construire p1s, p2s et le tableau serait mieux placée en dehors de cette boucle (que se passera-t-il si list_pnt contient plus de 2 points ?).
On devrait d'ailleurs pouvoir l'écrire plus simplement :
Code:
p1s = Point(p1.X, p1.Y, row.ZMIN) p2s = Point(p2.X, p2.Y, row.ZMIN)
La construction du tableau me semble correcte, je me demande si l'on ne peut pas écrire directement :
Code:
array.add(p1, p1s, p2s, p2, pt1)
Le polygone me semble lui aussi correctement construit, je ne comprends pas l'erreur. Que retourne polygon.getPart(0) ?
EDIT
En relisant le code je crois avoir identifié le problème. Le fait d'écrire p1s=p1 est à proscrire, car en Python les noms de variable sont comme des étiquettes sur un espace mémoire. En réalité avec cette écriture p1 et ps1 pointent sur la même valeur ainsi en modifiant l'un (cf. ligne p2s.Z=row.ZMIN) on modifie l'autre. Si l'on affiche les valeurs des 4 points à la fin du script on devrait s'apercevoir que p1=p1s et p2=p2s. Le polygone étant invalide il ne peut être créé.
En basant ton code sur mes suggestions précédentes, le script devrait être fonctionnel.
Dernière modification par dominique.lys (Sat 21 July 2012 19:38)
Hors ligne
#3 Mon 23 July 2012 11:03
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Salut Dominique,
Merci pour ta réponse.
Ce matin j'ai retravaillé sur le script en repartant de tes remarques et voici mon bilan :
- j'ai simplifié la récupération de p1 et p2
- pour les entités multipart je suis sur qu'il n'y en aura pas car j'ai fait un "chopper" égal à 2 sur mes données avec FME (Assure que toutes les entités sortant du Transformer ont un nombre de coordonnées inférieure ou égal au NBRE MAX DE VERTEX)
- pour la construction du tableau, le code : array.add(p1, p1s, p2s, p2, pt1) ne fonctionne pas du coup j'ai gardé mon code
- J'ai enlevé le p1s=p1 et p2s=p2 que j'ai remplacé par :
Code:
p1s = arcpy.Point(p1.X, p1.Y, row.ZMIN) >>> (ne pas oublier le arcpy.)
- le polygon.getPart(0) retourne une erreur : <geoprocessing array object object at 0x0CFEFBF0>
Après avoir fait toutes les modifications j'ai toujours le polygon.partCount(0). Par contre en modifiant une coordonnée par exemple en mettant p1s.X = p1X+0.001 alors j'obtiens un polygon.parCount(1) comme s'il ne pouvait pas créer de polygones verticaux.
Dernière version du code :
Code:
import arcpy rows=arcpy.UpdateCursor("W:/Florian/Contour.shp", "", "", "ZMIN") featureList = [] for row in rows: polyligne = row.shape list_pnt = polyligne.getPart(0) try: p1, p2 = list_pnt except: print "Erreur nombre de points superieur à 2" break #casse la boucle, on passe au feature suivant p1 = arcpy.Point() p2 = arcpy.Point() array = arcpy.Array() i = 0 for coordPair in list_pnt: if i==0: i=1 p1.X = coordPair.X p1.Y = coordPair.Y p1.Z = coordPair.Z print p1 else: p2.X = coordPair.X p2.Y = coordPair.Y p2.Z = coordPair.Z print p2 p1s = arcpy.Point(p1.X, p1.Y, row.ZMIN) p2s = arcpy.Point(p2.X, p2.Y, row.ZMIN) print p1s print p2s array.add(p1) array.add(p1s) array.add(p2s) array.add(p2) array.add(p1) # Create a Polygon object based on the array of points # polygon = arcpy.Polygon(array,None,True,True) print polygon.partCount print polygon.type # Clear the array for future use # array.removeAll() # Delete cursor and row objects to remove locks on the data # del row del rows
Florian Boret
Dream it, Make it, Share it
Hors ligne
#4 Mon 23 July 2012 11:37
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Il faut que tu vires toute la boucle sur list_pnt, elle ne sert à rien:
Code:
import arcpy rows=arcpy.UpdateCursor("W:/Florian/Contour.shp", "", "", "ZMIN") featureList = [] array = arcpy.Array() for row in rows: polyligne = row.shape list_pnt = polyligne.getPart(0) try: p1, p2 = list_pnt except: print "Erreur nombre de points superieur à 2" break #casse la boucle, on passe au feature suivant p1s = arcpy.Point(p1.X, p1.Y, row.ZMIN) p2s = arcpy.Point(p2.X, p2.Y, row.ZMIN) array.add(p1) array.add(p1s) array.add(p2s) array.add(p2) array.add(p1) # Create a Polygon object based on the array of points polygon = arcpy.Polygon(array, None, True, True) print polygon.partCount # Add polygon to list featureList.append(polygon) # Clear the array for future use array.removeAll() # Delete cursor and row objects to remove locks on the data del row del rows
array.add([p1, p1s, p2s, p2, pt1]) ne fonctionne pas non plus?
EDIT attention j'ai tapé trop vite il faut lire p1 à la place de pt1...
Juste après la création du polygone, peux-tu rajouter cette boucle :
Code:
for pt in polygon.getPart(0): print(pt.X, pt.Y, pt.Z)
comme ça on verra vraiment ce qu'il y a dans un polygone.
Quelque part ça ne m'étonnerai pas vraiment qu'un shapefile ne puisse pas gérer un polygone avec 2 points colinéaires en Z, c'est justement une des limites de la 2.5D : lorsqu'on "force" des structures de donnée et algo initialement conçu pour des données 2D dans un repère 3D. Un bon exemple est la triangulation de Delaunay appliquée à la modélisation de terrain --> impossible de trianguler une falaise qui serait parfaitement verticale.
Il faut vérifier ce point dans les spécifications du format.
EDIT
Après un test rapide avec la lib pyShp, j'arrive à créer une polygone dans le plan XZ avec colinéarité en Z, le problème n'est pas là.
Dernière modification par dominique.lys (Mon 23 July 2012 14:12)
Hors ligne
#5 Mon 23 July 2012 14:13
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
- enlever la boucle : OK
- array.add([p1, p1s, p2s, p2, p1]) : ne fonctionne pas
- la boucle polygon.getPart(0) ne fonctionne pas non plus ça me retourne un message d'erreur :
Code:
Traceback (most recent call last): File "W:\Florian\travail_sur_polygone_decale_points.py", line 50, in <module> for pt in polygon.getPart(0): File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\arcobjects.py", line 688, in getPart return convertArcObjectToPythonObject(self._arc_object.GetPart(*gp_fixargs(args))) RuntimeError: <unprintable RuntimeError object>
Pour ce qui est de la gestion des entités polygones avec 2 points colinéaires en Z, j'ai fait des tests en fin de matinée :
1 - j'ai généré un polygone 3D avec FME via le Transformer d'extrusion ça fonctionne mais le problème c'est que je n'arrive pas à attribuer un ZMIN avec cet outil.
2- j'ai chargé ce polygone dans ArcGIS qui me permet de le visualiser correctement.
3- j'ai modifié le polygone pour changer un vertex (X Y Z) et j'ai enregistré la modification.
4- j'ai re-modifié le polygone pour remettre les coordonnées de départ du vertex afin de rétablir le polygone avec 2 points colinéaires en Z et la impossible d'enregistrer : "La géométrie de l'entité ne peut être modifiée. La géométrie n'est pas valide.".
Bilan : ArcGIS ne permet pas de créer des entités polygones avec 2 points colinéaires en Z mais on peut charger ces entités si elles sont créées avec d'autres outils
Dernière modification par Floflo49fb (Mon 23 July 2012 14:21)
Florian Boret
Dream it, Make it, Share it
Hors ligne
#6 Mon 23 July 2012 17:06
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Je confirme avec ma version 9.3 et le module arcgisscripting, impossible de créer un polygone du type: [(0,0,0), (1,0,0), (1,0,1), (0,0,1), (0,0,0)]
Si on casse la colinéarité en Z, ça fonctionne: [(0,0,0), (1,0,0), (1,0.1,1), (0,0.1,1), (0,0,0)]
et avec ce type de polygone [(0,0,0), (1,0,0), (1,0,1), (0,0.1,1), (0,0,0)]
on se retrouve avec un triangle --> (1,0,0) et (1,0,1) ont été fusionnés car ils sont colinéaires en Z
si tu peux te permettre d'utiliser pr cette fonction une autre librairie je te conseille http://code.google.com/p/pyshp/
Hors ligne
#7 Mon 23 July 2012 21:11
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Un grand merci pour le coup de mains Dominique, je vais suivre ton conseil je vais passer sur la librairie pyshp.
Je reviendrais peut être vers toi si je bloque (Faut que je t'avoue, je débute sur python) !
Florian Boret
Dream it, Make it, Share it
Hors ligne
#8 Tue 24 July 2012 10:30
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Salut Dominique,
Après avoir fait des essais sur pyshp hier soir et ce matin, je reviens vers toi! J'ai bouquiné la documentation et le site http://geospatialpython.com j'ai compris comment ouvrir un shape lire les coordonnées des entités, voir le nom des champs,... mais la je patine!
Je voulais savoir si t'avais gardé le code de ton test rapide d'hier avec la lib pyShp. Si oui pourrais tu me le transmettre pour que je jettes un oeil au fonctionnement.
Encore merci.
Florian Boret
Dream it, Make it, Share it
Hors ligne
#9 Tue 24 July 2012 12:29
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Je te fais un peu du clé en main, mais bon ça te fera un bon exemple pr tes prochains développements:
Code:
# -*- coding:Latin-1 -*- import shapefile inShpPath="C:\\...\\pline.shp" outShpPath="C:\\...\\3dPoly.shp" fieldName="ZMIN" shp=shapefile.Reader(inShpPath) shapeRec = shp.shapeRecords()#extrait toute les géométries et données attributaires fields = shp.fields #chaque champs est décrit par un tuple de 3 valeurs (nom, type, longueur, précision) #print(fields) fieldsNames=[field[0] for field in fields if field[0] != 'DeletionFlag']#on recupere uniquement le nom des champs en exluant le champs system deletionflag #Note : cette syntaxe s'appel une "list compréhension", elle permet de synthétiser une boucle en 1 ligne #c'est un truc à bien comprendre et maitriser en python #print(fieldsNames) fieldIdx=fieldsNames.index(fieldName)#on cherche l'index du camps zmin #print(fieldIdx) newShp = shapefile.Writer(shapefile.POLYGONZ) newShp.fields = list(fields)#on calque la définition des champs du 1er shp for idx, entity in enumerate(shapeRec):#on parcours les entités geom=entity.shape.points #la geometrie (liste de points) values=entity.record #les données attributaires de cette entité try: p1, p2 = geom #print(p1, p2) except: print("Erreur nombre de points superieur à 2") break p1s = (p1[0], p1[1], values[fieldIdx]) p2s = (p2[0], p2[1], values[fieldIdx]) #print(p1s, p2s) polygon=[p1, p1s, p2s, p2, p1] #ajout du polygone et de ses données newShp.poly([polygon], shapeType=15)#list de list car polygon peut être multipart, shapetype15=POLYGONZ newShp.records.extend([values])#on copie les données initiales #Saving newShp.save(outShpPath)
Dernière modification par dominique.lys (Tue 24 July 2012 12:30)
Hors ligne
#10 Tue 24 July 2012 15:19
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Merci Dominique pour le script et toutes les explications jointes, c'est parfait! C'est clair que c'est une bonne référence pour d'autres scripts surtout que la librairie pyshp me parait assez simple à utiliser.
J'ai quand même modifié un peu ton code car on ne récupérait pas la valeur Z des points de la ligne, du coup je le remets en dessous pour ceux qui voudraient s'en servir.
Code:
import shapefile inShpPath="W:/Florian/Contours_toit" outShpPath="W:/Florian/Murs" fieldName="ZMIN" shp=shapefile.Reader(inShpPath) shapeRec = shp.shapeRecords()#extrait toute les géométries et données attributaires fields = shp.fields #chaque champs est décrit par un tuple de 3 valeurs (nom, type, longueur, précision) #print(fields) fieldsNames=[field[0] for field in fields if field[0] != 'DeletionFlag']#on recupere uniquement le nom des champs en exluant le champs system deletionflag #Note : cette syntaxe s'appel une "list compréhension", elle permet de synthétiser une boucle en 1 ligne #c'est un truc à bien comprendre et maitriser en python #print(fieldsNames) fieldIdx=fieldsNames.index(fieldName)#on cherche l'index du camps zmin #print(fieldIdx) newShp = shapefile.Writer(shapefile.POLYGONZ) newShp.fields = list(fields)#on calque la définition des champs du 1er shp for idx, entity in enumerate(shapeRec):#on parcours les entités geom=entity.shape.points #la geometrie (liste de points) z=entity.shape.z values=entity.record #les données attributaires de cette entité try: p1, p2 = geom z1, z2 = z #print(p1, p2) except: #print("Erreur nombre de points superieur à 2") break p1b = (p1[0], p1[1],z1) p2b = (p2[0], p2[1],z2) #print(p1b, p2b) p1s = (p1[0], p1[1], values[fieldIdx]) p2s = (p2[0], p2[1], values[fieldIdx]) #print(p1s, p2s) polygon=[p1b, p1s, p2s, p2b, p1b] #ajout du polygone et de ses données newShp.poly([polygon], shapeType=15)#list de list car polygon peut être multipart, shapetype15=POLYGONZ newShp.records.extend([values])#on copie les données initiales #Saving newShp.save(outShpPath)
PS : je viens de le lancer sur un fichier de lignes plus gros. Le résultat dans l'après midi...
Dernière modification par Floflo49fb (Tue 24 July 2012 17:10)
Florian Boret
Dream it, Make it, Share it
Hors ligne
#11 Tue 24 July 2012 16:53
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Ah oui j'avais oublié cette histoire pr récupérer les z, mais attention il me semble qu'avec shp.shapes()[0].z[0] tu récupères à chaque fois les z des points de la première polyligne de la couche. Tous tes points p1b p2b auront le même z.
Essaye peut-être plutôt comme ça :
Code:
geom=entity.shape try: p1, p2 = geom.points p1.append(geom.z[0]) p2.append(geom.z[1])
Hors ligne
#12 Tue 24 July 2012 17:13
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Oui je m'en suis rendu compte après avoir généré 1200 polygones!
Du coup, j'ai re-modifié le script au dessus mais je n'ai pas fait exactement ce que tu me proposes.
Je viens de re-générer les 1200 polygones et c'est ok, script validé!
Florian Boret
Dream it, Make it, Share it
Hors ligne
#13 Tue 24 July 2012 18:33
- dominique.lys
- Participant assidu
- Date d'inscription: 5 Oct 2006
- Messages: 473
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Content que ça fonctionne.
Avec le recul je me dis que tu aurais tout aussi bien pu rester sur ArcPy : en biaisant la colinéarité Z avec un petit nombre l'impact restait négligeable. J'ai un peu vendu ma soupe là.
Même si je travail beaucoup avec les outils ESRI, j’apprécie de pouvoir facilement bidouiller des shapes avec pyShp. C'est simple et efficace et petit à petit on se peut se constituer une bonne bibliothèque de fonctions réutilisables quelles que soient les solutions professionnelles choisies.
Python dispose de nombreux autres modules relatifs aux SIG dont les utilisateurs d'ArcGIS auraient torts de se priver. Comme par exemple Shapely qui peut facilement être couplé à pyShp et qui permet notamment de calculer des relations topologiques (intersection, union ...).
Si cela t'intéresse on trouve beaucoup de bonnes infos sur http://www.portailsig.org/
Bon courage pr la suite de ton projet.
Hors ligne
#14 Tue 24 July 2012 22:12
- Floflo49fb
- Participant assidu
- Lieu: Montpellier
- Date d'inscription: 29 Aug 2009
- Messages: 250
- Site web
Re: [python] Polyligne 3D vers polygone 3D
Je suis content du résultat, ça faisait plusieurs semaines que je cherchais désespérément une solution avec ArcGIS et FME!
C'est vrai que j'aurais pu rester sur Arcpy, j'y ai aussi pensé mais il aurait fallu mettre un décalage de 3 mm en Z aux polygones pour être dans les tolérances Esri d'après mes tests. Et pour quelqu'un qui a une formation de géomètre, 3 mm de décalage c'est déjà 3mm d'erreur!
Merci.
Florian Boret
Dream it, Make it, Share it
Hors ligne