[gelöst] Sammelalgorithmus in einem Raster

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Antworten
hephaistos
User
Beiträge: 7
Registriert: Donnerstag 3. Juli 2008, 15:24
Wohnort: Zürich, Davos

Montag 11. August 2008, 12:30

hallo zusammen
ich möchte eine Rasteroperation schreiben, die die Summe von Werten von direkt benachbarten Pixlel (acht an der zahl) einer Zelle zusammenzählt. Ist dieser Wert kleiner als ein bestimmtes abbruchkriterium (Holzmenge), soll die Nachbarschaft um Radius(r) 1 Grösser werden (16 involvierte Pixel).

Die Randpixelbehandlung möchte ich so gestalten, dass sobald ein Pixel ausserhalb des Rasters mit in die Berechnung des Zentralpixels kommt, dieses den einen nodatavalue oder einen 0-Wert bekommt.

Meine aktuellen Probleme dabei sind:
- die implementierung der randpixelbehandlung
- sobald ich einen anderen Startpixel angebe (xindex, yindex) stürzt IDLE ab.
-ich verstehe nicht, warum meine range von z.B. xmin bis xmax nur drei stellen besitzt.
in meinem beispiel unten müssten es ja 3 sein.
- einen eleganten weg zu finden, um das Zentralpixel nur einmal in die berechnung miteinzubeziehen.
- Radius r wird nicht grösser, auch wenn actualsum < Holzmenge
Ich habe mir dazu folgendes skript geschrieben:

Code: Alles auswählen

Holzmenge = 40
r = 1
xindex = 2
yindex = 2
xmin = xindex - r
xmax = xindex + r 
ymin = (yindex - r)# + 1
ymax = (yindex + r)# - 1
minilist_zeile = []
minilist_spalte = []
minilist_total = []
actualsum = 0

#zentralpixel addieren, geht nur, wenn man annimmt das zentralpixel < Holzmenge
minilist_zeile.append(array[xindex,yindex])

while actualsum <= Holzmenge:
    #zeile + randpixel jeder zeile
    for xmin in range (xmin, xmax):
        minilist_zeile.append(array[xindex-r, xmin])
        xmin = xmin + 1
      #spalte ohne randpixel jeder zeile
        for ymin in range (ymin, ymax+1):
            minilist_spalte.append(array[xindex, ymin])
            ymin = ymin + 1
    r = r + 1
    minilist_total = minilist_zeile + minilist_spalte
    actualsum = sum(minilist_total)
    
print "akutelle Summe: " + str(actualsum)
print "zeilen: "+str(minilist_zeile)
print "spalten: "+str(minilist_spalte)
#print minilist
print "Der Radius um den Holzvorrat bereitzustellen beträgt: " + str(r-1)

Danke für eure Hilfe...


Gruass
Hephaistos
Zuletzt geändert von hephaistos am Donnerstag 14. August 2008, 15:35, insgesamt 1-mal geändert.
BlackJack

Montag 11. August 2008, 13:42

Bei den Randpixeln musst Du halt einfach prüfen, ob beide Koordinaten innerhalb von `array` liegen. Vielleicht mit einer eigenen Funktion um Werte ab zu fragen, die eben entweder den Wert, oder 0 liefert. Am elegantesten wäre vielleicht eine eigene Klasse, die `__getitem__()` entsprechend implementiert.

Was heisst IDLE stürzt ab? Stürzt es wirklich ab, oder bekommst Du eine Ausnahme? Kannst Du ein komplettes Beispiel geben, dass man nachvollziehen kann?

Was meinst Du mit nur "drei Stellen" und es müssten 3 sein? Dir ist klar, dass die Grenzen bei `range()` und `xrange()` die Untergrenze in der Liste ist aber die Obergrenze nicht?

Das "Zentralixel" ist halt ein Sonderfall, den man am Anfang abfangen kann.

Der Radius `r` wird grösser, Du berücksichtigst das in der Schleife aber nirgends bei der Berechnung der Koordinaten.

Die ganzen `minilist_*`\s sind überflüssig, Du kannst da auch gleich die Werte auf summieren.

Ich bin mir ziemlich sicher, dass die beiden inneren ``for``-Schleifen nicht das tun, was Du möchtest. Das fängt damit an, dass die Laufvariable den gleichen Namen hat, wie die jeweilige `range()`-Untergrenze, geht darüber, dass es ja jeweils *zwei* Randzeilen und -spalten gibt, und endet damit, dass die Zeilen 22 und 26 keinen Effekt haben, weil der jeweilige Name im nächsten Schleifendurchlauf gleich wieder an einen anderen Wert gebunden wird.

Pack das Ganze am besten mal in eine oder mehrere Funktionen, die jeweils Teilaufgaben lösen.
Benutzeravatar
numerix
User
Beiträge: 2696
Registriert: Montag 11. Juni 2007, 15:09

Montag 11. August 2008, 13:46

Da bleibt für meinen Geschmack einiges unklar (oder ich stehe auf dem Schlauch):
hephaistos hat geschrieben:hallo zusammen
ich möchte eine Rasteroperation schreiben, die die Summe von Werten von direkt benachbarten Pixlel (acht an der zahl) einer Zelle zusammenzählt. Ist dieser Wert kleiner als ein bestimmtes abbruchkriterium (Holzmenge), soll die Nachbarschaft um Radius(r) 1 Grösser werden (16 involvierte Pixel).
Was verstehst du unter "Rasteroperation"?
Was ist eine Zelle?
Nachbarschaft? Um 1 größer -> 16 Pixel?
hephaistos hat geschrieben:Die Randpixelbehandlung möchte ich so gestalten, dass sobald ein Pixel ausserhalb des Rasters mit in die Berechnung des Zentralpixels kommt, dieses den einen nodatavalue oder einen 0-Wert bekommt.

Meine aktuellen Probleme dabei sind:
- die implementierung der randpixelbehandlung
- sobald ich einen anderen Startpixel angebe (xindex, yindex) stürzt IDLE ab.
-ich verstehe nicht, warum meine range von z.B. xmin bis xmax nur drei stellen besitzt.
in meinem beispiel unten müssten es ja 3 sein.
- einen eleganten weg zu finden, um das Zentralpixel nur einmal in die berechnung miteinzubeziehen.
- Radius r wird nicht grösser, auch wenn actualsum < Holzmenge
Was heißt "stürzt IDLE ab"?
"warum meine range von z.B. xmin bis xmax nur drei stellen besitzt
in meinem beispiel unten müssten es ja 3 sein." Hää?

Da du von "Radius" sprichst, geht es vermutlich hier um Kreise und um Punkte, die bei bestimmtem Radius noch innerhalb des Kreises liegen? Richtig?
hephaistos
User
Beiträge: 7
Registriert: Donnerstag 3. Juli 2008, 15:24
Wohnort: Zürich, Davos

Montag 11. August 2008, 14:32

Danke für eure prompten antworten

Zu euren Gegenfragen:

unter Rasteroperationen verstehe ich in diesem Fall ein Rechnungen mit Werten in einer Liste, die ich über eine spalten, und zeilenangabe ansprechen kann.
Ursprünglich handelt es sich dabei um ein ESRI Grid, welches ich mit gdal in die angesprochene liste geladen habe.

Eine Zelle ist also eine Art Container, den ich mit einem x und y wert ansprechen kann (xindex, yindex).

Unter Nachbarschaft verstehe ich bei Radius 1 die 8 benachbarten Pixel, bei Radius 2 , kommen dann wie numerix gesagt hat, 16 pixelwerte (16 weitere summanden) dazu.
Da du von "Radius" sprichst, geht es vermutlich hier um Kreise und um Punkte, die bei bestimmtem Radius noch innerhalb des Kreises liegen? Richtig?
genau!
Was heißt "stürzt IDLE ab"?
das heisst, dass wenn ich mein skript laufen lassen (und dies ev. fehlerhaft ist) der pythonprompt verschwindet. alle meine Funktionen kann ich nicht mehr sehen, bzw. editieren. danach muss ich mein idle mit sudo kill beenden, neustarten, die funktionen ausführen etc. das ist ziemlich zeitraubend.

die restlichen hinweise versuche ich umzusetzen. Danke für euren input.

Gruass
hephaistos
BlackJack

Montag 11. August 2008, 18:13

Musst Du echt mit ``sudo``(!) killen? Das wäre sehr merkwürdig.

Und wenn es wirklich um *Kreise* geht, dann ist diese Berechnung von Zeilen und Spalten auch nicht das richtige. Dein Quelltext sieht nämlich so aus, als wenn Du ein Quadrat immer grösser lassen werden wolltest.
hephaistos
User
Beiträge: 7
Registriert: Donnerstag 3. Juli 2008, 15:24
Wohnort: Zürich, Davos

Mittwoch 13. August 2008, 09:24

Guten Tag,

Ich habe an meinem "Algorithmus" ein weitergebastelt:
Zu euren Gegenfragen:
Musst Du echt mit ``sudo``(!) killen? Das wäre sehr merkwürdig.
nein, ich muss nicht mit sudo einen prozess "killen", das war ein fehler, den ich bis jetzt immer gemacht habe.
Und wenn es wirklich um *Kreise* geht, dann ist diese Berechnung von Zeilen und Spalten auch nicht das richtige. Dein Quelltext sieht nämlich so aus, als wenn Du ein Quadrat immer grösser lassen werden wolltest.
genau, es soll ein quadrat werden.
Bei den Randpixeln musst Du halt einfach prüfen, ob beide Koordinaten innerhalb von `array` liegen. Vielleicht mit einer eigenen Funktion um Werte ab zu fragen, die eben entweder den Wert, oder 0 liefert. Am elegantesten wäre vielleicht eine eigene Klasse, die `__getitem__()` entsprechend implementiert.
leider habe ich diesen hinweis nicht wirklich verstanden. Darum meine Frage 1: ich weuss nicht wie ich überprüfen soll, ob ein wert, welcher aufsummiert wird, sich ausserhalb des rasters befindet. Hier mein Versuch:
Der xindex,yindex sind die koordinaten des pixel für welches die umbebenden Werte aufsummiert werden.
xmin,xmax grenzt aufsummierenden rahmen ein.

Code: Alles auswählen

def Holzernter(ZuErntendeMenge, xindex, yindex):
    summe = 0
    summe = summe + array[xindex, yindex]
    r = 1
    
    while summe <= ZuErntendeMenge:
        print "r = " + str(r)
        print "summe " + str(summe) + " zem: " + str(ZuErntendeMenge) + "\n"
        xmin = xindex - r
        xmax = xindex + r 
        ymin = yindex - r
        ymax = yindex + r  
        #zeile min, zeile max, zeile min usw.
        if ((xmin >= 0) and (xmax >= 0)) or ((ymin >= 0) and (ymax >=0)):
            for i in range (xmin, xmax+1, 1):
                summe = summe + array[ymax, i] + array[ymin, i]
                print "y = " + str(i) + " value " + str(array[ymax, i])
                print "y = " + str(i) + " value " + str(array[ymin, i])
                #spalten
            for j in range (ymin+1, ymax, 1):
                print "x = " + str(i) + " value " + str(array[j, xmin])
                print "x = " + str(i) + " value " + str(array[j, xmax])
                summe = summe + array[j, xmin] + array[j, xmax]
            r = r + 1
            print summe
        else:
            print "HIER!"
            summe = summe + 0
            r = r + 1
            print summe
    radius = [r-1]
    return radius

Frage 2:
wie kann ich funktionen die bereits fehlerfrei kompiliert wurden einer anderen Funktion zugänglich machen. mit import funktionsname funktioniert das nicht. Alle Dateien befinden sich im selben verzeichnis. bei einem neustart von idle funktioniert die hauptfunktion erst wieder, wenn alle funktionen, die in der hauptfunktion aufgerufen werden neu kompiliert werden.

Frage 3:
Ich weiss, dass ich mit solchen Fragen hier schon nicht ganz richtig bin. Wo würdet ihr Hilfe, zu solchen, in euren Augen ev. trivialen Probleme suchen.

Gruss
BlackJack

Mittwoch 13. August 2008, 10:26

Ob Koordinaten innerhalb Deines `array` liegen, musst Du halt einfach überprüfen. Und zwar jedes mal und nicht nur für die Eckpunkte, denn es kann ja sein, dass das Quadrad teilweise noch im `array` liegt. Da ich nicht weiss was `array` für einen Typ hat, kann ich nicht sagen wie man die Höhe und Breite ermittelt, aber die brauchst Du natürlich zum Überprüfen der Koordinaten.

Da `r` in beiden Zweigen der ``if``-Abfrage um eins erhöht wird, kann man das auch *einmal* nach der Abfrage machen.

Zeile 29 ist recht sinnfrei. :-)

Warum steckst Du das Ergebnis in eine Liste?
hephaistos
User
Beiträge: 7
Registriert: Donnerstag 3. Juli 2008, 15:24
Wohnort: Zürich, Davos

Mittwoch 13. August 2008, 12:41

Saluti,
nun, meinen "array" kreiere ich mit:

Code: Alles auswählen

def createarrayfromgrid(gridname):
    # Importieren des Datensatzes
    dataset = gdal.Open(gridname, GA_ReadOnly )

    if dataset is None:
        print "Es wurde kein file gefunden, das importiert werden könnte"
    # was sagt mir dieser output
    print dataset

    #einige informationen über das importfile, welches oben spezifiziert wurde
    print 'Driver: ', dataset.GetDriver().ShortName,'/', \
              dataset.GetDriver().LongName
    print 'Size is ',dataset.RasterXSize,'x',dataset.RasterYSize, \
    'x',dataset.RasterCount
    print 'Projection is ',dataset.GetProjection()
        
    geotransform = dataset.GetGeoTransform()
    if not geotransform is None:
        print 'Origin = (',geotransform[0], ',',geotransform[3],')'
        print 'Pixel Size = (',geotransform[1], ',',geotransform[5],')'


    #ein GetRasterBand Objekt wird erzeugt
    band = dataset.GetRasterBand(1)

    print 'Band Type=',gdal.GetDataTypeName(band.DataType)

    #minimum und maximum aus dem file lesen
    min = band.GetMinimum()
    max = band.GetMaximum()
    if min is None or max is None:
        (min,max) = band.ComputeRasterMinMax(1)
        print 'Min=%.3f, Max=%.3f' % (min,max)
    if band.GetOverviewCount() > 0:
        print 'Band has ', band.GetOverviewCount(), ' overviews.'
    if not band.GetRasterColorTable() is None:
        print 'Band has a color table with ', \
        band.GetRasterColorTable().GetCount(), ' entries.'
    array = ()
    array = dataset.ReadAsArray()
    
    return array, dataset.RasterXSize, dataset.RasterYSize, geotransform[1]

danke für die Verbesserungen an meinem Code:
Zeile 29 ist recht sinnfrei. Smile
Sinngemäss möchte ich eben erreichen, dass für jede Zelle ausserhalb des Rasters, die Summe unverändert bleibt, der Radius aber erhöht wird und somit innerhalb des Rasters mehr Rasterwerte summiert werden können.

Danke für deine Hinweise!

Gruass
BlackJack

Mittwoch 13. August 2008, 16:27

Die Frage was `array` für ein Typ ist, habe ich jetzt mal per Suchmaschine geklärt. Hoffe ich. Es scheint ein `numpy.array` zu sein. Die haben ein `shape`-Attribut, das die Dimensionen enthält.

Hier völlig ungetesteter Quelltext:

Code: Alles auswählen

def get_value(array, x, y):
    """Returns value at x, y or 0 if x or y is not within the array."""
    width, height = array.shape
    return array[x, y] if 0 <= x < width and 0 <= y < height else 0


def wood_radius(array, needed_amount, start_x, start_y):
    radius = 0
    total = get_value(array, start_x, start_y)
    max_radius = max(array.shape)
    while total < needed_amount and radius <= max_radius:
        radius += 1
        min_x = start_x - radius
        max_x = start_x + radius
        min_y = start_y - radius
        max_y = start_y + radius
        for i in xrange(min_x, max_x + 1):
            total += get_value(array, i, min_y) + get_value(array, i, max_y)
        for i in xrange(min_y + 1, max_y):
            total += get_value(array, min_x, i) + get_value(array, max_x, i)
    if total < needed_amount:
        raise ValueError('not enough wood')
    return radius
Ich war so frei auch den Fall zu berücksichtigen, falls das gesamte `array` gar nicht genug Holz enthält.
hephaistos
User
Beiträge: 7
Registriert: Donnerstag 3. Juli 2008, 15:24
Wohnort: Zürich, Davos

Donnerstag 14. August 2008, 15:35

Hallo BlackJack.
Vielen Dank für dieses codesnippet. Es hat nun funktioniert. Anfangs lieferte mir die Funktion get_value noch einen falschen Datentyp. In meinem 2 dimensionalen numpy array hat es nämlich floats drin. somit musste ich als returnwert der funktion get_value auch einen float liefern, oder?

Die exception raise ValueError('not enough wood') ist hilfreich, danke!

Gruass
hephaistos
Antworten