banniere

Le portail francophone de la géomatique


Toujours pas inscrit ? Mot de passe oublié ?
Nom d'utilisateur    Mot de passe              Toujours pas inscrit ?   Mot de passe oublié ?

#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

 

Pied de page des forums

Powered by FluxBB