ArcGIS Rechenschritte erst mit Objekt 1 dann mit Objekt 2

[MarcNAV]

Habe ein größeres Programm geschrieben, das im Grunde eigentlich läuft. folgendes Problem habe ich aber: Sagen wir ich habe die Dateien a.asc und b.asc.

wenn ich die Rechenschritte printen lasse siehts in etwa so aus:

Projected: a.asc
Projected: b.asc
Rastered: b.asc
Clipped: b
Zonal Statistics: b


Ich verstehe nicht warum das Programm nach dem Projezieren nur mit b weitermacht. Weiß jemand Rat?

Besser wäre es wohl wenn das Programm erst alle Rechenschritte mit a.asc rechnet und dann mit b.asc.

Code: Alles auswählen

import shutil as st
import arcpy, os, sys, glob
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
from os.path import basename


def main():

#Dateistamm folgen, Ordner aus Datei erzeugen und Datei hinein kopieren

    root_path = 'F:/test/PythonClippingTree'         #!!!!!!!!!!Eingabe
    extension = '.asc'
    for path, _, filenames in os.walk(root_path):
        for filename in filenames:
            if filename.endswith(extension):

                source = os.path.join(path, filename)
                target = source[:-len(extension)]
                os.mkdir(target)
                st.move(source, target)
                target2 = target + "_rastered"
                target3 = target + "_clipped"
                target4 = target + "_zonalstat"
                os.mkdir(target2)
                os.mkdir(target3)
                os.mkdir(target4)

###################################################################
#Projektion

                arcpy.env.workspace = (target)
                arcpy.env.overwriteOutput = True
                # set local variables
                coordinateSystem = 'PROJCS["Bessel_1841_Lambert_Conformal_Conic",GEOGCS["GCS_Bessel_1841",DATUM["D_Bessel_1841",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Standard_Parallel_1",50.0],PARAMETER["Standard_Parallel_2",51.0],PARAMETER["Latitude_Of_Origin",47.0],UNIT["Meter",1.0]]'
                arcpy.DefineProjection_management(filename, coordinateSystem)
                print 'Projected: ' + os.path.basename(filename)

###################################################################
 #.asc in Raster wandeln

    files = glob.glob(target + '/*.asc')
    for f in files:

            inASCII = f
            x = os.path.splitext(os.path.basename(f))[0]
            outRaster = os.path.relpath(target2) + '/%s'%x
            rasterType = ("INTEGER")
            arcpy.ASCIIToRaster_conversion(inASCII, outRaster, rasterType)
            print 'Rastered: ' + os.path.basename(f)

###################################################################
#Clipping:


    # input workspace
    raster_ws = target2
    shp_ws = (r'F:/test/PythonClipping/CalcPythClipZonal/Schablone')     #!!!!!!!!!!Eingabe
    arcpy.env.overwriteOutput = True

    # output workspace
    out_ws = target3


    # get list raster
    arcpy.env.workspace = raster_ws
    rasters = [os.path.join(raster_ws, r) for r in arcpy.ListRasters()]

    # get shapefile
    arcpy.env.workspace = shp_ws
    shapef = (r'F:/test/PythonClipping/CalcPythClipZonal/Schablone/schab.shp')   #!!!!!!!!!!Eingabe


    # get dictionary {raster1 : shapefile.shp, ...}
    clip_dict = dict(zip(rasters, shapef))


    # clip all rasters
    outputList = [] # container
    for raster, shapefile in clip_dict.iteritems():
        out_raster = os.path.join(out_ws, os.path.basename(raster))
        arcpy.Clip_management(raster, "", out_raster, shapef, "", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
        print 'Clipped: ' + os.path.basename(raster)
        outputList.append(out_raster) # apend all the output to the container


###################################################################
#Zonal Statistics:

    # Check out any necessary licenses
    arcpy.CheckOutExtension("spatial")

    # Input data source
    arcpy.env.workspace = target3
    arcpy.env.overwriteOutput = True

    ZoneShapefile = ("F:/test/PythonClipping/BL_Grenzen/l_d_bundneu_polygon.shp")

    # Output File
    OutputFolder = target4

    # Loop through a list of files in the workspace
    RasterFiles = arcpy.ListRasters()

    # Local variables:
    for filename in RasterFiles:
        print("Zonal Statistics: " + filename)
        inRaster = arcpy.env.workspace + "/" + filename
        fileroot = filename
        outRaster = OutputFolder + "/" + fileroot + ".dbf"

        # Process: Zonal Statistics
        arcpy.gp.ZonalStatisticsAsTable(ZoneShapefile, "GRID_CODE", inRaster, outRaster, "DATA", "ALL")

###################################################################
#Ordner alle in target verschieben


    st.move(target2, target)
    st.move(target3, target)
    st.move(target4, target)


if __name__ == '__main__':
    main()

[BlackJack]

@MarcNAV: Weil der entsprechende Code offensichtlich nicht *in* der Schleife, sondern *nach* der Schleife steht. Einrückung hat in Python eine wichtige Bedeutung. Daran sieht der Compiler was zu einer Schleife, einem ``if``/``else``, Funktionen, und anderen Blockstrukturen gehört und was nicht.

Die Kommentare aus '####'…-Zeichen sollte man auch entsprechend einrücken, damit die die Übersicht nicht erschweren. Das ist dem Compiler zwar egal, aber erschwert menschlichen Lesern den Überblick.

Letztendlich ist das auch alles ziemlich lang und viel für eine Funktion. Das sollte man vielleicht sinnvoll auf mehrere Funktionen aufteilen.

[MarcNAV]

stimmt so einfach kann die Lösung manchmal sein. Dachte ich muss da noch komplizierter rangehen mit ner Laufvariablen...... Vielen Dank jetzt funktionierts

[BlackJack]

@MarcNAV: Noch ein paar Anmerkungen zum Quelltext:

`shutil` würde ich nicht mit `st` abkürzen. Der Modulname ist nicht allzu lang und wird auch nicht oft benutzt. Das ist es nicht Wert den Quelltext dafür unverständlicher zu machen weil man gegebenfalls erst nachsehen muss was `st` ist.

Das `sys`-Modul wird importiert, aber nirgends verwendet. Das gleiche gilt für `env` aus `arcpy` und die `basename()`-Funktion aus `os.path`. Auch der Sternchenimport aus `arcpy.sa` scheint nichts zu importieren was irgendwo benutzt wird.

Die `main()`-Funktion ist wie schon gesagt deutlich zu lang und zu verschachtelt um das gut nachvollziehen zu können. Dadurch werden auch einige Zeilen zu lang. Der
Style Guide for Python Code empfiehlt 80 Zeichen pro Zeile.

Zum Aufteilen auf Funktionen könnte man als erstes das finden der Dateinamen in eine eigene Funktion heraus ziehen. Damit bekommt man schon mal zwei Einrückebenen aus der `main()`-Funktion heraus.

Die Durchnummerierten `target*`-Namen sind unschön weil die Nummer dem Leser nichts vermittelt und er erst die Definition suchen muss um zu erfahren was `target3` eigentlich *bedeutet*. Das ist insbesondere in dieser überlangen Funktion nicht gut wo die Namen sehr weit oben definiert, dann aber erst in den Teilschritten und ganz am Ende wieder verwendet werden. Man könnte die Definitionen näher an die Teilschritte verschieben in denen der jeweilige Name als erstes benötigt wird und dann beim Aufteilen in Unterfunktionen mit in die jeweilige Unterfunktion verschieben.

Muss man `arcpy.env.overwriteOutput` tatsächlich immer wieder auf `True` setzen? Reicht es nicht das einmal am Anfang des Programms zu machen? Falls nicht, sollte man es auch nicht ganz am Anfang machen, sondern immer kurz bevor es gebraucht wird.

Diese überlange `coordinateSystem`-Zeichenkette kann man wesentlich lesbarer formatieren.

Bei der Ausgabe am Ende der Projektion ist das `os.path.basename()` überflüssig, denn `filename` enthält ja nur den Dateinamen und nicht auch einen Pfad davor.

An drei Stellen im Quelltext sind überflüssige Klammern um einzelne Werte.

Man sollte immer `os.path.join()` zum Verbinden von Pfadteilen verwenden.

Namen sollten beschreiben was der Wert dahinter bedeutet. `files` und `filenames` sind zwei unterschiedliche Sachen, die man auch richtig benennen sollte um Verwirrungen zu vermeiden. Und einbuchstabige Namen ausserhalb von sehr beschränkten Gültigkeitsbereichen wie „list comprehensions”, Generatorausdrücken, oder ``lambda``-Definitionen sind in der Regel nicht ausreichend.

Die Definition von `inASCII` beim Konvertieren von ASC zu Rasterdateien ist sinnlos.Das gleiche gilt für `fileroot` im letzten Berechnungsschritt.

Was beim Clipping mit `clip_dict` veranstaltet wird ist totaler Unsinn. Die Werte sind nicht das was der Kommentar behauptet und sie werden dann auch überhaupt gar nicht verwendet. Die Werte in `outputList` werden auch nur in der Liste gesammelt, dann aber überhaupt nicht verwendet.

``print`` ist in Python 2 keine Funktion, darum sollte man das auch nicht mit Klammern so schreiben als wäre es eine.

Ich würde dann ungefähr bei so etwas heraus kommen (ungetestet):

Code: Alles auswählen

#http://www.python-forum.de/viewtopic.php?t=34196&p=260303#p260303
import shutil
import glob, os
import arcpy
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension('Spatial')


def iter_filenames(root, extension):
    for path, _, filenames in os.walk(root):
        for filename in filenames:
            if filename.endswith(extension):
                yield path, filename


def project(target, filename):
    arcpy.env.workspace = target
    arcpy.env.overwriteOutput = True
    # set local variables
    coordinate_system = (
        'PROJCS['
            '"Bessel_1841_Lambert_Conformal_Conic",'
            'GEOGCS['
                '"GCS_Bessel_1841",'
                'DATUM['
                    '"D_Bessel_1841",'
                    'SPHEROID["Bessel_1841",6377397.155,299.1528128]'
                '],'
                'PRIMEM["Greenwich",0.0],'
                'UNIT["Degree",0.0174532925199433]'
            '],'
            'PROJECTION["Lambert_Conformal_Conic"],'
            'PARAMETER["False_Easting",0.0],'
            'PARAMETER["False_Northing",0.0],'
            'PARAMETER["Central_Meridian",9.0],'
            'PARAMETER["Standard_Parallel_1",50.0],'
            'PARAMETER["Standard_Parallel_2",51.0],'
            'PARAMETER["Latitude_Of_Origin",47.0],'
            'UNIT["Meter",1.0]'
        ']'
    )
    arcpy.DefineProjection_management(filename, coordinate_system)
    print 'Projected:', filename


def convert_asc_to_raster(target):
    rastered_target = target + '_rastered'
    os.mkdir(rastered_target)
    for asc_file_path in glob.glob(os.path.join(target, '*.asc')):
        asc_file_basename = os.path.splitext(os.path.basename(asc_file_path))[0]
        raster_file_path = os.path.join(
            os.path.relpath(rastered_target), asc_file_basename
        )
        arcpy.ASCIIToRaster_conversion(
            asc_file_path, raster_file_path, 'INTEGER'
        )
        print 'Rastered:', os.path.basename(asc_file_path)
    return rastered_target


def clip(target, rastered_target):
    arcpy.env.workspace = rastered_target
    arcpy.env.overwriteOutput = True

    raster_filenames = arcpy.ListRasters()

    arcpy.env.workspace = 'F:/test/PythonClipping/CalcPythClipZonal/Schablone'

    clipped_target = target + '_clipped'
    os.mkdir(clipped_target)
    # Clip all raster files.
    for raster_filename in raster_filenames:
        arcpy.Clip_management(
            os.path.join(rastered_target, raster_filename),
            '',
            os.path.join(clipped_target, raster_filename),
            'F:/test/PythonClipping/CalcPythClipZonal/Schablone/schab.shp',
            '',
            'ClippingGeometry',
            'NO_MAINTAIN_EXTENT'
        )
        print 'Clipped:', raster_filename

    return clipped_target


def calculate_zonal_statistics(target, clipped_target):
    zonalstat_target = target + '_zonalstat'
    os.mkdir(zonalstat_target)

    # Check out any necessary licenses
    arcpy.CheckOutExtension('spatial')
    arcpy.env.workspace = clipped_target
    arcpy.env.overwriteOutput = True

    for raster_filename in arcpy.ListRasters():
        print 'Zonal Statistics:', raster_filename
        arcpy.gp.ZonalStatisticsAsTable(
            'F:/test/PythonClipping/BL_Grenzen/l_d_bundneu_polygon.shp',
            'GRID_CODE',
            os.path.join(arcpy.env.workspace, raster_filename),
            os.path.join(zonalstat_target, raster_filename + '.dbf'),
            'DATA',
            'ALL'
        )

    return zonalstat_target


def main():
    #
    # Dateistamm folgen, Ordner aus Datei erzeugen und Datei hinein kopieren
    #
    root_path = 'F:/test/PythonClippingTree'
    extension = '.asc'
    for path, filename in iter_filenames(root_path, extension):
        source = os.path.join(path, filename)
        target = source[:-len(extension)]
        os.mkdir(target)
        shutil.move(source, target)

        project(target, filename)
        rastered_target = convert_asc_to_raster(target)
        clipped_target = clip(target, rastered_target)
        zonalstat_target = calculate_zonal_statistics(target, clipped_target)

        shutil.move(rastered_target, target)
        shutil.move(clipped_target, target)
        shutil.move(zonalstat_target, target)


if __name__ == '__main__':
    main()

[MarcNAV]

Danke, dass du dir so Mühe gemacht hast Black Jack. Ich habe dein Programm getestet und es läuft einwandfrei und sieht 10mal aufgeräumter aus;) Noch eine Frage:

Ich möchte noch gerne, dass wenn eine Datei beschädigt ist, er sie unverändert lässt und mit den anderen weiter macht.

Hab das bisher so gemacht:

Code: Alles auswählen

def main():
    try:
        Bearbeitungsschritte

    except:
        pass
        print 'Processing failed: ' + filename

if __name__ == '__main__':
    main()

Das blöde ist nur, dass ermir dann trotzdem lauter Ordner aus der beschädigten Datei macht und auch das Programm abbricht.

[EyDu]

Zeig mal den gesamten Code, dein Beispiel ist zu abstrakt. Wahrscheinlich versetzt du dein Programm irgendwo in einen ungültigen Zustand.

[MarcNAV]

Sorry war zu undeutlich formuliert. Ich meinte den Code von BlackJack. Für diesen hätte ich gerne noch eine Ausnahmebedingung für beschädigte .asc Dateien. Sollte eine .asc Datei beschädigt sein, dann soll er diese nicht weiter beachten(nur printen 'Processing failed', filename) und mit der Nächsten fortfahren. zerbrech mir schon ne Stunde den Kopf und weiß nich weiter.