Я использую Python и QGIS 2.0. Я пытаюсь обрезать растры в папке одним полигоном. Это первый раз, когда я использую (скажем, «PyQGIS»), я раньше привык к Arcpy. В любом случае, мой простой сценарий не работает, любое предложение будет высоко оценено!
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()
Ниже приведены улучшения, которые я сделал с тех пор, пока сценарий не работал, но я думаю, что могу быть ближе ...
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)
Я думаю, что может быть что-то не так в команде "gdal", так как функция "print" выполняет свою работу в основном, но файл не записывается в вывод, и я не получаю никакой ошибки. Кстати, трудно найти простую документацию по кодированию gdal ...
CLIP
Вcmd
выражении является проблемой. Если вы поместите переменную в строку, она не будет прочитана. Вместо этого вы должны объединить строку с переменной.print(cmd)
вместоos.system(cmd)
. ВашаoutRaster
переменная неверна.Ответы:
Я согласен с Натаном. Вы должны питонизировать весь ваш скрипт. Поэтому замените ваш
for
цикл на что-то вроде следующего:Примечание 1: я предполагаю, что ваши растровые файлы - GeoTIFF (
*.tif
).Примечание 2:
-of GTiff
не требуется вcmd
, потому что это формат вывода по умолчанию вgdalwarp
.источник
os.system
, ты прав.Я наконец справился с этим очень простым и чистым скриптом, который вызывает GDAL из Python без его импорта (как предложено, но с использованием метода «Call ()» вместо «os.system ()». Надеюсь, это могло бы помочь!
источник
Модифицированная версия решения umbe1987 для пользователей Linux:
источник