Schleifenproblem

[Janjan]

Ja moin,

langsam bin ich mit meinem Latein am ende... habe selber in arcgis(geografische Informationssoftware) polylinien digitalisiert und möchte Entferungen der einzelnen Koordinaten einer Polylinie berechnen. Basiert alles auf dem sog. WGS84-Ellipsoid.

Was ich nicht verstehe, warum die unten markierte Schleife exakt immer ZWEIMAL durchlaufen wird und der nächste Datensatz(es sind mehr als zwei) nich auch berechnet wird, kann da kein logikfehler erkennen... die folgenden polylinien in dem shapefile sind haargenau vom selben typ wie die ersten beiden.

Code: Alles auswählen

# -*- coding: cp1252 -*-
import arcgisscripting, math
gp = arcgisscripting.create(9.3)
gp.Workspace = r"C:\Users\XXX\Desktop\pipeline_shape\\"
fc = "eigene_pipes.shp"
desc = gp.Describe(fc)
rows = gp.SearchCursor(fc)
row = rows.Next()
feat = row.shape

line_x = [] #Neue Liste für X-Koordinaten
line_y = [] #Neue Liste für Y-Koordinaten
line_xy = [line_x,line_y]
distanzen = []

def polyline():      #Funktion um die einzelnen Koordinaten(X,Y) zu bekommen und zu speichern
    a = 0
    while a < feat.PartCount:
        array = feat.GetPart(a)
        array.Reset
        pnt = array.Next()
        #print pnt
        while pnt:
            print pnt.x, pnt.y
            line_x.append(pnt.x)
            line_y.append(pnt.y)
            pnt = array.Next()
        a = a + 1

def berechnung(): #Formel für Entferungsmessung basierend auf dem WGS84-Ellipsoid
    pi = 3.14159265
    f = 0.003352811
    a = 6378.137
    F = (line_x[0]+line_x[1])/2.0
    G = (line_x[0]-line_x[1])/2.0
    l = (line_y[0]-line_y[1])/2.0
    #print "F: "+str(F),"G: "+str(G),"l: "+str(l)

    #Umrechnung in Bogenmaß (Radiant),Python erwartet sin,cos wert als radiant/bogenmaß
    F = (pi/180.0)*F
    G = (pi/180.0)*G
    l = (pi/180.0)*l
    #print "F: "+str(F),"G: "+str(G),"l: "+str(l)
    S = (math.sin(G)*math.sin(G))*(math.cos(l)*math.cos(l)) + (math.cos(F)*math.cos(F))*(math.sin(l)*math.sin(l))
    #print "S: "+str(S)
    C = (math.cos(G)*math.cos(G))*(math.cos(l)*math.cos(l)) + (math.sin(F)*math.sin(F))*(math.sin(l)*math.sin(l))
    #print "C: "+str(C)
    w = math.atan(math.sqrt(S/C))
    #print "w: "+str(w)
    D = 2.0*w*a
    #print "D: "+str(D)
    R = (math.sqrt(S*C))/w
    #print "R: "+str(R)
    H1 = (3.0*R-1)/(2.0*C)
    #print "H1: "+str(H1)
    H2 = (3.0*R+1)/(2.0*S)
    #print "H2: "+str(H2)
    dist = D*(1.0 + f*H1*((math.sin(F)*math.sin(F)))*(math.cos(G)*math.cos(G))) - (f*H2*(math.cos(F)*math.cos(F))*(math.sin(G)*math.sin(G)))
    print "Distanz: "+str(round(dist))+" km"
    distanzen.append(dist)
    

for xy in line_xy:
      b = 0
      polyline()
      while b == 0:
          berechnung()
          line_x.remove(line_x[0])    #erste indexposition löschen, damit nächster abstand berechnet werden kann
          line_y.remove(line_y[0])
          if len(line_x) == 1:
              line_x.remove(line_x[0])   #Liste leeren für nächste Polyline
              line_y.remove(line_y[0])
              b = b + 1
            
              row = rows.Next()
      

distanzen.sort()
print distanzen[0:]
distanzen.reverse()
print distanzen[0:]

Vielleicht sieht ja jemand was ich nicht sehe. wie gesagt läuft zweimal durch, die Koordinaten sind korrekt, verglichen mit Ansicht in der Software.

Ausgabe ungefähr so:
x,y
x,y
x,y
x,y
Distanz: 123.44 km
Distanz: 432.22 km
Distanz: 12.44 km
Distanz: 65.22 km
[12.44,65.22,123.44,432.22]
.
.
.
.

[0x1cedd1ce]

Setz bitte deinen Code in Python-Tags, ansonsten da bleiben dann auch die Einrückungen erhalten. Aber auch ohne Einrückungen habe ich deinen Fehler gefunden.
Du hast die Liste line_xy, diese Liste enthält die Listen line_x und line_y, das sind 2 Stück. Daher wird die Schleife nur zweimal durchlaufen.
Nutz doch erst einmal deine polyline-Funktion um eine Liste mit allen Polylines zu erstellen, die kannst du dann in der For-Schleife verwenden. Dort brauchst du dann auch nicht mehr Polyline aufrufen. Das würde auch diese grauenhafte prozedurale Programmierung aufhübschen. (Grauenhaft weil Prozedural)

[BlackJack]

@Janjan: Um das mit dem prozedural == grauenhaft noch einmal zu präzisieren: Ganz kommt man um prozedurale Programmierung nicht herum in Python, aber man sollte in Funktionen nicht auf globalen Datenstrukturen operieren. Das macht das Programm unübersichtlich und damit fehleranfälliger.

Werte, die man in Funktionen verwendet (ausser Konstanten) sollten als Argumente an die Funktionen übergeben werden und das Ergebnis sollte als Rückgabewert von der Funktion an den Aufrufer zurück gegeben werden. So sieht man besser was eine Funktion benötigt und liefert und kann sie auch einfacher einzeln mit ein paar Werten testen, was die Fehlersuche erleichtert und die Chance der Wiederverwendbarkeit erhöht.

Ist `Reset` auf dem was `GetPart` von einem Feature zurück gibt tatsächlich nur ein Attribut, oder ist das eine Methode die man auch *aufrufen* sollte?

Eine ``while``-Schleife in der Code für die Abbruchbedingung nur eine Variable hochzählt schreibt man besser/einfacher als ``for``-Schleife. Die ``while``-Schleifen über die Iteratoren kann man mit der `iter()`-Funktion ebenfalls in einfachere ``for``-Schleifen umwandeln.

Ich habe auch nicht ganz verstanden warum Du die Koordinaten auf zwei Listen aufteilst, denn sie sollen doch wohl gemeinsam verarbeitet werden!? Also sollte man sie vielleicht auch besser zusammen in einer Liste speichern.

Die Kommentare zu den Funktionen sollten nicht "inline" stehen, sondern als DocStrings unter der Kopfzeile der Funktion.

Und Abkürzungen sollte man möglichst vermeiden, wenn sie nicht allgemein bekannt sind. `feature` statt `feat` oder `point` statt `pnt` machen ein Programm leichter lesbar und verständlicher.

Die `polyline()`-Funktion könnte dann zum Beispiel so aussehen:

Code: Alles auswählen

def get_polyline(feature):
    """Funktion um die einzelnen Koordinaten (X,Y) von einem Feature
    abzufragen.
    
    Gibt eine Liste mit den Koordinaten als Tupel zurück.
    """
    result = list()
    for i in xrange(feature.PartCount):
        points = feature.GetPart(i)
        points.Reset()
        result.extend((p.x, p.y) for p in iter(points.Next, None))
    return result

Eine Konstante für π sowie Umrechungsfunktionen zwischen Grad und Bogenmass sind schon im `math`-Modul enthalten.

Statt ``str(round(dist))`` solltest Du besser mit Zeichenkettenformatierung arbeiten.

Was Du da in der Hauptschleife unten machst ist mir auf den ersten Blick nicht ganz klar, da möchte ich auch nicht weiter drüber nachdenken, denn eine Liste verarbeiten, in dem man mit `remove()` immer am Anfang Elemente entfernt ist keine gute Idee. Selbst wenn man `pop()` dafür verwenden würde, wird es vom Vorgehen her nicht besser. Hier sollte man eine Schleife schreiben die immer zwei aufeinander folgende Listenelemente nimmt und als Argumente an eine Funktion übergibt, welche die Distanz berechnen und die Liste selbst sollte man dabei nicht verändern.

Das Anfordern des Iterators über die Zeilen am *Anfang* des Quelltextes um den dann am Ende, nach zwei längeren Funktionsdefinitionen, endlich zu benutzen, trägt nicht zur Übersichtlichkeit bei. Am besten lässt Du allen Quelltext auf Modulebene, der zusammen die Hauptfunktion ausmacht, in einer Funktion verschwinden.

[Janjan]

das mit dem reset ist mir selber nicht ganz klar, habe den polyline code ausm netz, ich selber stecke auch noch in den anfängen. aber ja habt beide natürlich recht. math.pi hab gar nich dran gedacht...

Besten Dank für Eure SCHNELLEN Tipps!!

[Lenzer]

Habt ihr die Lösung gefunden? Ich finde das Forum an sich sehr interessant, da habe ich mir selber schon viele Lösungen herausgeholt

[Janjan]

Code: Alles auswählen

# -*- coding: cp1252 -*-
import arcgisscripting, math
gp = arcgisscripting.create(9.3)
gp.OverwriteOutput = True
gp.Workspace = r"C:\Users\XXX\Desktop\pipeline_shape\\"
featureclass = "eigene_pipes.shp"
path = r"C:\Users\XXX\Desktop\pipeline_shape\eigene_pipes.shp"
rows = gp.SearchCursor(featureclass,"","","Shape")
row = rows.Next()
feature = row.shape

polyline = []
distanzen = []

def berechnung(): #Formel für Entferungsmessung basierend auf dem WGS84-Ellipsoid
    f = 0.003352811
    a = 6378.137
    F = (polyline[0][0] + polyline[1][0])/2.0
    G = (polyline[0][0] - polyline[1][0])/2.0
    l = (polyline[0][1] - polyline[1][1])/2.0
    #print "F: "+str(F),"G: "+str(G),"l: "+str(l)

    #Umrechnung in Bogenmaß (Radiant),Python erwartet sin,cos wert als radiant/bogenmaß
    F = (math.pi/180.0)*F
    G = (math.pi/180.0)*G
    l = (math.pi/180.0)*l
    #print "F: "+str(F),"G: "+str(G),"l: "+str(l)
    S = (math.sin(G)*math.sin(G))*(math.cos(l)*math.cos(l)) + (math.cos(F)*math.cos(F))*(math.sin(l)*math.sin(l))
    #print "S: "+str(S)
    C = (math.cos(G)*math.cos(G))*(math.cos(l)*math.cos(l)) + (math.sin(F)*math.sin(F))*(math.sin(l)*math.sin(l))
    #print "C: "+str(C)
    w = math.atan(math.sqrt(S/C))
    #print "w: "+str(w)
    D = 2.0*w*a
    #print "D: "+str(D)
    R = (math.sqrt(S*C))/w
    #print "R: "+str(R)
    H1 = (3.0*R-1)/(2.0*C)
    #print "H1: "+str(H1)
    H2 = (3.0*R+1)/(2.0*S)
    #print "H2: "+str(H2)
    dist = D*(1.0 + f*H1*((math.sin(F)*math.sin(F)))*(math.cos(G)*math.cos(G))) - (f*H2*(math.cos(F)*math.cos(F))*(math.sin(G)*math.sin(G)))
    #print "Distanz: "+str(round(dist))+" km"
    distanzen.append(dist)
    return dist

def get_polyline(featureclass):
    """Funktion um die einzelnen Koordinaten (X,Y) von einem Feature
    abzufragen.
   
    Gibt eine Liste mit den Koordinaten als Tupel zurück.
    """
    for i in xrange(feature.PartCount):
        points = feature.GetPart(i)
        points.Reset()
        polyline.extend((p.x, p.y) for p in iter(points.Next, None))
        polyline.append("polyline_end")
        row = rows.Next()
    return polyline

polyline_copy = []
'''
Kopie der Liste "polyline" über welche die Schleife läuft
'''
count_rows = int(gp.GetCount_management(gp.describe(featureclass).catalogpath).getoutput(0))
print count_rows
'''
Öffnen des Shapefiles um die vorhandenen Reihen auszuzählen
'''
for row in range(0,count_rows):
    #row = rows.Next()
    get_polyline(featureclass)
    polyline_copy.extend(polyline)
    for xy in polyline_copy[:]:
        if polyline[1] == "polyline_end":
            polyline.remove(polyline[1])
            polyline.remove(polyline[0])
            distanzen.append("polyline_end")
            break
        else:
            berechnung()
            polyline.remove(polyline[0])

print polyline[:], distanzen[:]

Ja hiermit gehts, ist jetzt noch ohne Module und noch "dreckige" Verarbeitung der polyline liste, aber erstmal das was ich ursprünglich wollte.

Ergebnis ist dann:
5
[] [442.09738758010224, 140.11909616984204, 263.87739870197476, 240.11303429593511, 180.62064972160584, 382.81689563507547, 384.33671455231962, 'polyline_end', 1008.1779130704485, 1367.5596733401865, 'polyline_end', 1086.2681880631606, 877.80512125477537, 2383.7159795399152, 399.27538793245236, 'polyline_end', 1706.9842726917416, 1998.070083022114, 386.64282921804647, 1113.0946209158358, 814.99958104890277, 2880.1615310393586, 'polyline_end', 936.36879948287719, 1230.2239082469459, 272.38153374403493, 801.43678146492368, 'polyline_end']