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 ...
CLIP
nacmd
expressã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ávelprint(cmd)
no lugar deos.system(cmd)
. SuaoutRaster
variável não está correta.Respostas:
Estou de acordo com Nathan. Você precisa pythonizar todo o seu script. Portanto, substitua seu
for
loop por algo como o seguinte:Nota 1: Suponho que seus arquivos de varredura sejam GeoTIFF (
*.tif
).Nota 2:
-of GTiff
não é necessáriocmd
, porque é o formato de saída padrão emgdalwarp
.fonte
os.system
, você está certo.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!
fonte
Versão modificada da solução de umbe1987 `para usuários de Linux:
fonte