Le portail francophone de la géomatique

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 :

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 :

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 :

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 :

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)

Il faut que tu vires toute la boucle sur list_pnt, elle ne sert à rien:

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 :

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)

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/

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.

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.

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)

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é!

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.