Я пишу простую утилиту для обрезки пакетов многоканальных растровых файлов геотифов в одну и ту же (меньшую) область. Используя gdalwarp, я легко могу обрезать файл, используя шейп-файл с одним полигоном:
gdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif
Однако фактическая область, в которую я хочу обрезать, всегда будет первоначально определяться другим растровым файлом геотифов, а не шейп-файлом. Было бы хорошо, если бы я мог использовать экстент этого растра в качестве файла отсечения, но я не уверен, как это сделать. Неудивительно, что следующее не работает (оно не вызывает ошибку, оно просто ничего не производит):
gdalwarp -cutline clipper.tif-crop_to_cutline input.tif output.tif
Итак, мой вопрос, есть ли способ снабдить растром gdalwarp -cutline
? Альтернативно, есть ли другая функция gdal, которая может обрезать растр, используя другой растр? Если ни то, ни другое невозможно, существует ли очень простой способ создания шейп-файла с одним полигоном, определяемым экстентом растра?
Этот код будет заключен в более обширный скрипт на python, поэтому я могу использовать утилиты gdal для командной строки или любые привязки python для gdal.
В качестве примечания, я знаю, что я мог бы просто сделать отсеченный шейп-файл, который охватывает размер моего растра в QGIS. Я могу сделать это, если не найду простого решения, но в конечном итоге я буду использовать эту утилиту на десятках, если не на сотнях областей, как часть большого автоматизированного анализа, поэтому я предпочел бы не делать утомительного ручной шаг, даже если это очень легко.
gdaltindex clipper.shp clipper.tif
, а затемgdalwarp -cutline clipper.shp -crop_to_cutline input.tif output.tif
Для неправильных полигонов и при условии, что ваш растровый файл geotiff является двоичным, вы можете использовать GDAL_Calc :
Этот запрос будет заполнять 0, где Mask.tif <= 0 и BigImage, где Mask> 0. Для этого оба растра должны иметь одинаковый размер ячейки, строки и столбцы. Чтобы извлечь те же экстенты, используйте GDAL_Translate с
-projwin ulx uly lrx lry
параметром (прямоугольник находится в проецируемых координатах), но убедитесь, что прямоугольник projwin не распространяется по краям любого растра.Подставьте значения для поля projwin, полученного из Mask.
источник
Решение в Python напрямую, без придания формы:
источник