Boucles "for" imbriquées (python)

[FAT_FOUILLEUL]

Bonjour à tous,

Je travaille actuellement sur un script python pour arcgis, dont la première étape est de calculer pour chaque point la distance du point le plus proche.

Je bute depuis quelques heures sur une erreur dont je ne saisis pas le sens. Le code me paraît logique, mais peut être qu'un œil extérieur pourra m'aider à comprendre ce qui ne marche pas, d'autant plus que mon cerveau n'est plus très frais - et que mes compétences en programmation sont assez limitées -

voila le code donc :

import arcpy
from math import *

# Récupère la classe d'entités ou la couche en entrée
inFeatureClass = arcpy.GetParameterAsText(0)
valeur = arcpy.GetParameterAsText(1)
outFeatureClass = arcpy.GetParameterAsText(2)

# Récupérer les données XY
arcpy.AddXY_management(inFeatureClass)

# Ajouter un champ pour calculer distance au point le plus proche
arcpy.AddField_management(inFeatureClass, "distance", "Float")

#Définir les curseurs
cS = arcpy.UpdateCursor(inFeatureClass)
cU = arcpy.UpdateCursor(inFeatureClass)
rg = arcpy.Row

# effectue une boucle sur chaque des enregistrements
for rg in cS:
    X = rg.getValue("POINT_X")
    Y = rg.getValue("POINT_Y")
    #initialiser les variables de distance
    distMin = valeur
    dist = 0.0
    #Pour chaque enregistrement, on compare la distance avec ses voisins
    for rg in cU:
        X2 = rg.getValue("POINT_X")
        Y2 = rg.getValue("POINT_Y")
        #Calcul de la distance
        dist = sqrt(((X2-X)**2)+((Y2-Y)**2))
        if dist < distMin and dist <> 0:
            distMin = dist
            print distMin
    rg.setValue("distance", distMin)
    cS.updateRow(rg)

Le résultat me donne quelque chose comme ceci (50000 est valeur maximale de recherche en entrée):
alors que normalement la distance devrait être calculée pour chacun des enregistrements. Je me demande donc comment gérer les deux boucles "for" ici imbriquées car j'ai l'impression que la seconde boucle (qui compare chaque point aux autres) n'est exécutée qu'une seule fois alors qu'elle devrait être relancée à chaque passage de la première boucle.

N'étant pas expert en python, je me demande comment gérer cette histoire, si par exemple cela provient d'une bête erreur de syntaxe ou même de mon algorithme. Merci d'avance pour toutes vos suggestions.

[Lud]

for rg in cs:
ici rg prend la valeur une à une des valeurs de cs
for rg in cu
pareil et donc tu écrases le premier rg
du coup quand tu fais le rg.setvalue, c'est bien le rg du cu qui prend la valeur ...
Pas sur que ça règle tes histoires de boucles mais à mon avis y a un soucis là. (pour résoudre cela suffit de mettre for rg2 in cu)

[gene]

Lud a raison mais, puisque vous voulez apprendre Python, il y a aussi un problème de conception, puisque la distance entre le point 1 et le point 2 est la même que celle entre le point 2 et le point 1.
Je m'explique:
avec votre script:

premier point de la liste Cs: 
     parcourir la liste Cu du premier point au dernier
deuxième point de la liste Cs:
     parcourir la liste Cu du premier point au dernier
...
dernier point de la liste Cs:
     parcourir la liste Cu du premier point au dernier

ce qui vous oblige a utiliser

if dist < distMin and dist <> 0:

pour éliminer les cas distance(premier point, premier point) etc.
De plus, vous passez n fois par tous les éléments de la liste Cu (parcourir la liste Cu du premier point au dernier)

alors que ce que vous devez faire c'est:
1) n'utiliser qu'une seule liste (Cs par exemple)
2) ne passer qu'une fois par les éléments de la liste

premier point de la liste Cs: 
     parcourir la liste Cs du deuxième point au dernier
deuxième point de la liste Cs:
     parcourir la liste Cs du troisième  point au dernier (puisque dist 1-2 ou 2-1 est déjà connue)
dernier point de la liste Cs:
     rien, puisque toutes les distances dernier - autres points on déjà été calculées

cela se fait avec le slicing ou cadrage des listes que j'ai expliqué dans http://www.forumsig.org/showthread.php?t=33402

Je n'utilise pas ArcGis et ne suis pas sur Windows mais je connais Python et son utilisation dans QGIS ou GRASS GIS, par exemple:

• et donc pour votre script, la première chose à faire est d'écarter le calcul de la distance de la boucle pour en faire une fonction (gain de performance):

def dist(x1,y1,x2,y2):
      return sqrt(((x2-x1)**2)+((y2-y1)**2))

• ensuite si je choisis une liste du genre de ce qu'ArcPy va vous envoyer

Cs= [[232278.92351962379, 117357.2537612508], [233788.23928931833, 117067.46
513346945], [234609.30706803215, 115703.0436776656], [231204.29069160129, 114797
.45421584888], [229972.68902353055, 115799.63988692604]]

• le script devient donc

# première boucle avec enumerate pour avoir un index i à utiliser dans la deuxième boucle
for i, valeur in enumerate(Cs):
     x1 = valeur[0]
     y1 = valeur[1]
     print "index: ", i, "valeurs: ", x1,y1
     #deuxième boucle avec cadrage c.a.d. valeurs après i + 1 seulement
     for j in cs[i+1:]:
         x2=j[0]
         y2=j[1]
         print "valeurs x,y traitées: ", x2,y2
         # utilisation de la fonction distance
         print  "distance:",dist(x1,y1,x2,y2)

• résultat:

index:  0 valeurs x,y:  232278.92352 117357.253761
valeurs x,y traitées:  233788.239289 117067.465133
distance: 1536.8837111
valeurs x,y traitées:  234609.307068 115703.043678
distance: 2857.81358443
valeurs x,y traitées:  231204.290692 114797.454216
distance: 2776.22215028
valeurs x,y traitées:  229972.689024 115799.639887
distance: 2782.96218667

index:  1 valeurs x,y:  232278.92352 117357.253761
valeurs x,y traitées:  234609.307068 115703.043678
distance: 2857.81358443
valeurs x,y traitées:  231204.290692 114797.454216
distance: 2776.22215028
valeurs x,y traitées:  229972.689024 115799.639887
distance: 2782.96218667

index:  2 valeurs x,y:  232278.92352 117357.253761
valeurs x,y traitées:  231204.290692 114797.454216
distance: 2776.22215028
valeurs x,y traitées:  229972.689024 115799.639887
distance: 2782.96218667
valeurs x,y traitées:  232278.92352 117357.253761
distance: 0.0

index:  3 valeurs x,y:  232278.92352 117357.253761
valeurs x,y traitées:  229972.689024 115799.639887
distance: 2782.96218667

index:  4 valeurs x,y:  232278.92352 117357.253761

Conclusions

1. gain de temps: on ne crée pas 2 listes comme dans votre cas
2. gain de temps: on ne repasse pas n fois par tous les éléments de la liste
3. aucune erreur possible puisque vous n'utiliser qu'une fois les éléments présents dans la seule liste nécessaire
4. et donc aucune erreur du genre de celle que vous avez dans votre script, possible

J'espère que ces quelques remarques vous seront utiles dans votre apprentissage de Python, tâche souvent oubliées par les utilisateurs de ArcGIS...

[FAT_FOUILLEUL]

Je vous remercie pour les réponses même si je me suis un peu cassé les dents toute la journée sur le concept de liste^^

Au final l'algorithme que je cherche ne correspondait pas à ce que tu m'as proposé Gene (même si il marche impeccablement ) car il faut que je compare à chaque point courant les points qui ont déjà été passé en revue. Je ne cherche pas ainsi la distance la plus courte sur l'ensemble des points, mais cette distance point par point. Cela dit ma méthode est probablement assez lourde, mais je préfère toujours écrire quelque chose que je comprend.

Donc au final mon code qui me renvoie enfin les résultats escomptés:

import arcpy
from math import *

# Récupère la classe d'entités ou la couche en entrée
inFeatureClass = arcpy.GetParameterAsText(0)
tolerance = arcpy.GetParameterAsText(1)
outFeatureClass = arcpy.GetParameterAsText(2)

# Récupérer les données XY
arcpy.AddXY_management(inFeatureClass)

# Ajouter un champ pour calculer distance au point le plus proche
arcpy.AddField_management(inFeatureClass, "distance", "Float")

# Initialiser les curseurs
cS = arcpy.UpdateCursor(inFeatureClass)
rg2 = arcpy.Row
Clist = []

# Données contenues dans une liste
for rg2 in cS:
    Clist = Clist + [[rg2.getValue("POINT_X"), rg2.getValue("POINT_Y")]]

# Recherche point par point du point le plus proche
for i, valeur in enumerate(Clist):
    x1 = valeur[0]
    y1 = valeur[1]
    dMin = tolerance
    # Boucle pour comparer les distance entre les points
    for j in Clist:
        x2 = j[0]
        y2 = j[1]
        D = dist(x1,y1,x2,y2)
        if D < dMin and D <> 0:
            dMin = D
    print dMin

# Fonction de distance
def dist(x1,y1,x2,y2):
    return sqrt(((x2-x1)**2)+((y2-y1)**2))

enfin le script n'est pas fini, mais je décrète la journée finie

[gene]

même si je me suis un peu cassé les dents toute la journée sur le concept de liste^^

Mais tu devrais, c'est un des fondements des traitements en Python.

De plus, je vois que tu as repris ce que je t'ai proposé en l'ajustant à ce que tu voulais faire et c'est comme ça qu'on apprend.

Comme je le répète toujours, Python n'est pas qu'un langage de script d'ArcGIS mais c'est un langage de programmation complet dont ArcPy n'est qu'un module parmi des milliers d'autres. Une fois qu'on le maîtrise, le reste coule de source, quelque soit le module.

[Barbapapa]

Si tu veux optimiser ton script, tu peux faire joujou avec les k-d tree, cf
http://en.wikipedia.org/wiki/Kd-tree
Où tu as même quelques bouts de code en python.

Mais si tu as peu de points, l’intérêt est franchement limité.

[FAT_FOUILLEUL]

Je vous remercie pour toutes vos réponses, cela m'a bien aidé. J'ai désormais terminé mon script (outil de clustering selon un paramètre de distance en entrée) mais je rencontre un problème quand je veux l'ajouter dans ma toolbox. Le script donne en effet les résultats escomptés lorsque je l’exécute via la console python, mais en revanche il ne ressort qu'un seul point quand je lance l'outil par la toolbox ; comme si tout mes points avaient été fusionnés.

Exemple en image :

Carte avant le traitement



Carte après traitement effectué sur la console python



Carte obtenue en utilisant la toolbox (le script est le même)



Je me demande donc si cela peut venir des paramètres (GetValueIntoText) en entrée, mais à priori non car le script s’exécute tout de même sans messages d'erreur. J'avoue sécher, donc si quelqu'un peut voir d'où vient cette différence entre les deux rendus, je lui en serais extrêmement reconnaissant.

Pour information, voila le code :

#Import des modules
import arcpy
from math import *

#Paramètres en entrée
inFeatureClass = arcpy.GetParameterAsText(0)
tolerance = arcpy.GetParameterAsText(1)
outFeatureClass = arcpy.GetParameterAsText(2)
chpval = arcpy.GetParameterAsText(3)

#Copie de la couche en sortie
FeatureClass = arcpy.CopyFeatures_management(inFeatureClass, outFeatureClass)

#Ajout des données XY
arcpy.AddXY_management(FeatureClass)

#On initialise les curseurs avant d'importer les données dans le programme
cS = arcpy.UpdateCursor(FeatureClass)
rg2 = arcpy.Row
Clist = []
objId = 1
dMin = 1

#En même temps que l'on importe les données (localisation et valeur) de la table
#celle-ci est effacée pour laisser place aux nouvelles coordonnées
for rg2 in cS:
    Clist = Clist + [[objId, rg2.getValue("POINT_X"), rg2.getValue("POINT_Y"), rg2.getValue(chpval)]]
    objId = objId + 1
    cS.deleteRow(rg2)

#Boucle tant qu'aucune distance entre chaque point est inférieure à la tolérance
while (dMin == tolerance) is False:
        dMin = tolerance
        #On essaye ainsi de récupérer la distance minimale à chaque fois
        for i, valeur in enumerate(Clist):
                x1 = valeur[1]
                y1 = valeur[2]
                #Boucle pour comparer lces distance entre les points
                for j in Clist[i+1:]:
                        x2 = j[1]
                        y2 = j[2]
                        # Appel à la fonction de distance
                        D = dist(x1,y1,x2,y2)
                        # Soit dMin la distance au point le plus proche
                        if D < dMin and D <> 0:
                                dMin = D
                                #ptA et ptC les points correspondants à cette distance
                                ptA = Clist.index(Clist[i])
                                ptC = Clist.index(j)
        #On ajoute une nouvelle liste fruit des deux points fusionnés
        #tandis que ces deux points sont supprimés de la liste
        Clist.append([0, (Clist[ptC][1] + Clist[ptA][1])/2, (Clist[ptC][2] + Clist[ptA][2])/2, (Clist[ptC][3] + Clist[ptA][3])])
        Clist.remove(Clist[ptA])
        #Manip pour se retrouver dans la liste
        if ptC > ptA:
                Clist.remove(Clist[ptC-1])
        else:
                Clist.remove(Clist[ptC])
    
print Clist

#Ajout du curseur pour écrire dans la table
rowInsert = arcpy.InsertCursor(outFeatureClass)

#Ecriture des points à partir des coordonnées générées.
for row in Clist:
    inPoint = arcpy.CreateObject("Point")
    inPoint.X = row[1]
    inPoint.Y = row[2]
    newPoint = rowInsert.newRow()
    newPoint.SHAPE = inPoint
    newPoint.setValue(chpval, row[3])
    rowInsert.insertRow(newPoint)
    
del cS
del rg2
del rowInsert

#La fonction de distance
def dist(x1,y1,x2,y2):
    return sqrt(((x2-x1)**2)+((y2-y1)**2))