Zum Schleifen von Ordnern zum Stapeln von Clip-Rastern nach Polygon mit Python und QGIS?

9

Ich benutze Python und QGIS 2.0. Ich versuche, Raster in einem Ordner nach einem Polygon-Feature zu schneiden. Es ist das erste Mal, dass ich "PyQGIS" benutze (sagen wir mal). Wie auch immer, ich bekomme mein einfaches Skript nicht zum Laufen, jeder Vorschlag wäre sehr dankbar!

import qgis.core, qgis,utils
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()

CLIP= "C:/Users/unim/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/study_area_foscagno.shp"
INPUT_FOLDER="C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/LE71930282000259EDC00"
OUTPUT= "C:/Users/unim/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/foscagno_pyqgis/"


for RASTER in INPUT_FOLDER.tif
do
    echo "Processing $RASTER"
    gdalwarp -q -cutline CLIP -crop_to_cutline -of GTiff RASTER OUTPUT+ "clip_"+ RASTER
done

QgsApplication.exitQgis()

Im Folgenden sind die Verbesserungen aufgeführt, die ich seitdem vorgenommen habe, ohne dass das Skript funktioniert, aber ich denke, ich komme vielleicht näher ...

import qgis.core, qgis.utils, os, fnmatch
from osgeo import gdal

CLIP= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/study_area_foscagno.shp"
INPUT_FOLDER= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/LE71930282000259EDC00/DNs2Reflectance_LE71930282000259EDC00"
OUTPUT= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/Cloud_mask_AltaValtellina/clip_2_foscagno"

def findRasters (path, filter):
    for root, dirs, files in os.walk(path):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

for raster in findRasters (INPUT_FOLDER, '*.tif'):
    print (raster)
    outRaster = OUTPUT + '/clip_' + raster
    cmd = 'gdalwarp -dstnodata 0 -q -cutline CLIP -crop_to_cutline %s %s' % (raster, outRaster)
    os.system (cmd)

Ich denke, dass der Befehl "gdal" möglicherweise nicht stimmt, da die Funktion "print" ihre Aufgabe ordnungsgemäß erfüllt, aber keine Datei in die Ausgabe geschrieben wird und ich auch keinen Fehler erhalte. Übrigens war es schwierig, eine einfache Dokumentation zur GDAL-Codierung zu finden ...

umbe1987
quelle
Zunächst einmal mischen Sie Python und Bash mit GDAL-Skripten. Sie können dies nur mit gdal tun oder müssen Sie pyqgis verwenden?
Nathan W
Vielen Dank, ich würde gerne Python verwenden, da dies nur ein Ausgangspunkt für ein größeres Skript wäre. Ist es möglich, es wie bei arcpy mit einer Problemumgehung zu verwenden?
umbe1987
Das CLIPim cmdAusdruck ist das Problem. Wenn Sie eine Variable in eine Zeichenfolge einfügen, wird diese nicht gelesen. Stattdessen würden Sie die Zeichenfolge mit der Variablen verketten.
Antonio Falciano
Ich benutze es jetzt draußen, es gibt keinen Fehler aus und es "druckt" alle ".tif" -Raster richtig. Nachdem ich einige Dinge erledigt habe (z. B. viermal weniger als eine Sekunde pro Fenster geöffnet habe), erhalte ich keine Ausgabe in meinem OUTPUT-Ordner.
umbe1987
Überprüfen Sie die Rasterpfade mit print(cmd)anstelle von os.system(cmd). Ihre outRasterVariable ist nicht korrekt.
Antonio Falciano

Antworten:

9

Ich stimme Nathan zu. Sie müssen Ihr gesamtes Skript pythonisieren. Ersetzen Sie Ihre forSchleife also durch Folgendes:

import os, fnmatch

def findRasters (path, filter):
    for root, dirs, files in os.walk(path):
        for file in fnmatch.filter(files, filter):
            yield file

for raster in findRasters(INPUT_FOLDER, '*.tif'):
    inRaster = INPUT_FOLDER + '/' + raster
    outRaster = OUTPUT_FOLDER + '/clip_' + raster
    cmd = 'gdalwarp -q -cutline %s -crop_to_cutline %s %s' % (CLIP, inRaster, outRaster)
    os.system(cmd)

Hinweis 1: Ich gehe davon aus, dass Ihre Rasterdateien GeoTIFF ( *.tif) sind.
Hinweis 2: -of GTiff Wird in nicht benötigt cmd, da dies das Standardausgabeformat in ist gdalwarp.

Antonio Falciano
quelle
Dankeschön. Es heißt jedoch "os.command (cmd) AttributeError: 'module'-Objekt hat kein Attribut' command '", obwohl das "os" -Modul importiert wurde ...
umbe1987
Es sollte sein os.system, du hast recht.
Antonio Falciano
4

Ich habe es endlich mit diesem sehr einfachen und sauberen Skript geschafft, das GDAL aus Python aufruft, ohne es zu importieren (wie vorgeschlagen, aber mit der Methode "Call ()" anstelle von "os.system ()". Ich hoffe, das könnte helfen!

import os, fnmatch
from subprocess import call
call(["ls", "-l"])

inFolder= 'C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/LE71930282000259EDC00/DNs2Reflectance_LE71930282000259EDC00/'
outFolder= 'C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/Cloud_mask_AltaValtellina/clip_2_foscagno/'

os.chdir (inFolder)

def findRasters (path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

for raster in findRasters (inFolder, '*.tif'):
    (infilepath, infilename)= os.path.split (raster)
    print infilename
    outRaster= outFolder+ 'clip_'+ infilename
    print outRaster
    warp= 'gdalwarp -dstnodata 0 -q -cutline %s -crop_to_cutline -of GTiff %s %s' % ('study_area_foscagno.shp', raster, outRaster)
    call (warp)
umbe1987
quelle
4

Geänderte Version der umbe1987- Lösung für Linux-Benutzer:

import os, fnmatch
from subprocess import call
call(["ls", "-l"])

inFolder= '/run/media/user/SOFT/LANDSAT/Bulk Order 595257/L8 OLI_TIRS/LC81840262015165LGN00/'
outFolder= '/run/media/user/SOFT/LANDSAT/Bulk Order 595257/summer_clipped/'
shp = '/run/media/user/SOFT/LANDSAT/Bulk Order 595257/vector/mask.shp'
os.chdir (inFolder)

def findRasters (path, filter):
    for root, dirs, files in os.walk(path, filter):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

for raster in findRasters (inFolder, '*.TIF'):
    (infilepath, infilename)= os.path.split (raster)
    print infilename
    outRaster= outFolder+ 'clip_'+ infilename
    print outRaster
    warp= 'gdalwarp -cutline \'%s\' -crop_to_cutline -dstalpha \'%s\' \'%s\'' % (shp, raster, outRaster)
    os.system(warp)
RustyDrill
quelle