Do zapętlania folderów do grupowania klipów rastrowych według wielokąta za pomocą Pythona i QGIS?

9

Używam Pythona i QGIS 2.0. Próbuję przyciąć rastry do folderu według jednej funkcji wielokąta. To pierwszy raz, gdy używam (powiedzmy) „PyQGIS”, wcześniej byłem przyzwyczajony do gry. Tak czy inaczej, mój prosty skrypt nie działa, każda sugestia byłaby bardzo mile widziana!

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()

Poniżej znajdują się ulepszenia, które wprowadziłem od tej pory, ale nie uruchamiają skryptu, ale myślę, że mogę być coraz bliżej ...

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)

Myślę, że może być coś złego w poleceniu „gdal”, ponieważ funkcja „print” wykonuje swoje zadanie prawidłowo, ale żaden plik nie jest zapisywany na wyjściu, ani nie otrzymuję żadnego błędu. Nawiasem mówiąc, trudno było znaleźć łatwą dokumentację na temat kodowania gdal ...

umbe1987
źródło
Na początek miksujesz Python i Bash ze skryptami gdal. Możesz to zrobić po prostu za pomocą gdal czy potrzebujesz pyqgis?
Nathan W
dziękuję, chciałbym użyć Pythona, ponieważ byłby to tylko punkt wyjścia do większego skryptu. Czy można go używać tak jak w przypadku Arcpy z pewnym obejściem?
umbe1987,
CLIPW cmdekspresji jest problem. Jeśli wstawisz zmienną do łańcucha, nie zostanie ona odczytana. Zamiast tego połączysz łańcuch ze zmienną.
Antonio Falciano
Używam go teraz na zewnątrz, nie wyświetla żadnych błędów i odpowiednio „drukuje” wszystkie rastry „.tif”. Jednak po zrobieniu kilku rzeczy (takich jak otwieranie 4 razy na mniej niż sekundę okna), nie otrzymuję żadnych danych wyjściowych w moim folderze OUTPUT.
umbe1987,
Sprawdź ścieżki rastrowe z print(cmd)zamiast os.system(cmd). Twoja outRasterzmienna jest nieprawidłowa.
Antonio Falciano

Odpowiedzi:

9

Zgadzam się z Nathanem. Musisz spytonować cały skrypt. Więc zamień swoją forpętlę na coś takiego:

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)

Uwaga 1: Przypuszczam, że twoje pliki rastrowe to GeoTIFF ( *.tif).
Uwaga 2: -of GTiff nie jest wymagana w cmd, ponieważ jest to domyślny format wyjściowy w gdalwarp.

Antonio Falciano
źródło
Dziękuję Ci. Mówi jednak: „os.command (cmd) AttributeError: obiekt„ module ”nie ma atrybutu„ command ””, chociaż moduł „os” został zaimportowany ...
umbe1987,
Powinno być os.system, masz rację.
Antonio Falciano
4

W końcu poradziłem sobie z tym bardzo prostym i czystym skryptem, który wywołuje GDAL z Pythona bez importowania go (zgodnie z sugestią, ale przy użyciu metody „Call ()” zamiast „os.system ()”. Mam nadzieję, że to może pomóc!

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
źródło
4

Zmodyfikowana wersja rozwiązania umbe1987 dla użytkowników systemu Linux:

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
źródło