Para pasta de loop para rasters de clipe em lote por polígono usando python e QGIS?

9

Estou usando python e QGIS 2.0. Estou tentando cortar rasters em uma pasta por um recurso de polígono. É a primeira vez que eu uso (digamos) "PyQGIS", eu já estava acostumado ao arcpy antes. De qualquer forma, não consigo que meu script simples funcione, qualquer sugestão seria muito apreciada!

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

Abaixo estão as melhorias que fiz desde agora, mas não consegui que o script funcionasse, mas acho que posso estar chegando mais perto ...

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)

Eu acho que pode haver algo errado no comando "gdal", pois a função "print" faz seu trabalho corretamente, mas nenhum arquivo é gravado na saída, nem recebo nenhum erro. By the way, tem sido difícil gostar de uma documentação fácil sobre codificação gdal ...

umbe1987
fonte
Bem, para começar, você está misturando Python e bash com scripts gdal. Você pode fazer isso usando gdal ou precisa usar o pyqgis?
Nathan W
obrigado, eu gostaria de usar o Python, pois isso seria apenas um ponto de partida para um script maior. É possível usá-lo como fiz com o arcpy com alguma solução alternativa?
umbe1987
O CLIPna cmdexpressão é o problema. Se você colocar uma variável em uma string, ela não será lida. Em vez disso, você concatenaria a string com a variável
Antonio Falciano 22/10
Estou usando-o fora agora, ele não gera nenhum erro e "imprime" todas as rasters ".tif" corretamente. No entanto, depois de fazer algumas coisas (como abrir 4 vezes por menos de um segundo por janela), não recebo nenhuma saída na minha pasta OUTPUT.
Umb1987
Verifique os caminhos de varredura com print(cmd)no lugar de os.system(cmd). Sua outRastervariável não está correta.
Antonio Falciano 23/10

Respostas:

9

Estou de acordo com Nathan. Você precisa pythonizar todo o seu script. Portanto, substitua seu forloop por algo como o seguinte:

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)

Nota 1: Suponho que seus arquivos de varredura sejam GeoTIFF ( *.tif).
Nota 2: -of GTiff não é necessário cmd, porque é o formato de saída padrão em gdalwarp.

Antonio Falciano
fonte
obrigado. No entanto, ele diz "os.command (cmd) AttributeError: o objeto 'module' não tem nenhum atributo 'command'", embora o módulo "os" tenha sido transportado ...
umbe1987 22/13
Deveria estar os.system, você está certo.
Antonio Falciano 22/10
4

Finalmente consegui gerenciar esse script muito simples e limpo, que chama GDAL do Python sem importá-lo (como sugerido, mas usando o método "Call ()" em vez de "os.system ()". Espero que isso ajude!

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
fonte
4

Versão modificada da solução de umbe1987 `para usuários de 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
fonte