Пакетная обработка gdal (в bash)

[mak-vas]

Это конечно больше по линуксу вопрос, чем по gdal, но всё же. Есть набор файлов raster01.tif - raster10.tif. Есть соответствующие им файлы-маски (mask[01-10]). Выходные файлы - raster_warped[01-10]. Вот команда:

Код: Выделить всё

gdalwarp -cutline mask??.shp raster??.tif raster_warped??.tif

Но, естественно, это не работает =), т.к. аутпут-файл не понимает "*" и "?" маски в баше.

[Александр Мурый]

А если сделать цикл типа такого?

Код: Выделить всё

for i in $(seq -w 01 10); do gdalwarp -cutline mask$i.shp raster$i.tif raster_warped$i.tif; done

[mak-vas]

Отлично, работает! Но позвольте задать ещё пару вопросов.

1. Если имена растров представлены вот так вот: P-43-1,2, P-43-3,4 (сдвоенные листы топокарт), то чтобы обойти каждый файл в цикле (писать скрипты на баше не умею, уж простите -я про логику тут думаю, синтаксис команды-то конечно неправильный) надо сделать что-то вроде:

Код: Выделить всё

for i in $(seq -w 0 1); do gdalwarp P-43-$(i*2+1),$(i*2+2).tif P43-$(i*2+1),$(i*2+2)_warped.jpg; done

Собственно вопрос в том, как это корректно сделать.

2. Вот в этой теме (последнее собщение) вы говорили про создание геотифа из гиф+мап.
Меня интересует копирование разрешения файла по x и y с помощью gdalinfo.
Есть ли альтернатива этому методу, которую можно было бы автоматизировать (как в вашем цикле выше)?

[Voltron]

Меня интересует копирование разрешения файла по x и y с помощью gdalinfo.
Есть ли альтернатива этому методу, которую можно было бы автоматизировать (как в вашем цикле выше)?

Посмотрите в сторону связки GDAL+Python (вводная статья). В принципе, можно размеры растра и на bash получить, только муторно.

[Александр Мурый]

* Voltron
Каждый питон.. тьфу, каждый кулик своё болото хвалит.

Код: Выделить всё

:~$ gdalinfo raster.jpg | grep "Size is" | cut -d" " -f3,4 | sed -e 's/,//'

1000 1500

Честно, не вижу, что здесь муторного. Фильтруем вывод gdalinfo (grep), разделяем на части и выводим только 3-ю и 4-ю (cut), убираем запятую (sed). Всё в одну строку.
Можно и так:

Код: Выделить всё

:~$ gdalinfo raster.jpg | grep "Size is" | tr -d 'Size is' | tr ',' ' '

1000 1500

В shell-скрипте помещаем в переменную, например:

Код: Выделить всё

size=$(gdalinfo raster.jpg | grep "Size is" | cut -d" " -f3,4 | sed -e 's/,//')

Можно увидеть такой же вывод размера на питоне? Статью читал, но прошу понимающих помочь (надо учиться).

* mak-vas
Такая логика, как у вас написано, не пойдёт (seq тут не нужен, в одну строку тоже не сделать). Лучше исп-ть утилиту find для рекурсивного обхода каталога и поправки с помощью gdalwarp всех растров.
Прикладываю подобный скрипт, написанный когда-то "на коленке". Скрипт надо поместить в каталог с растрами, сделать исполняемым (chmod +x gdal_proj_check.sh), запустить (./gdal_proj_check.sh). Русские имена файлов и пробелы и прочие нехорошие знаки перевариваются. На выходе -- поправленные с учётом их размера растры по шаблону "имя_растра__warped.tif".

gdal_proj_check.gz

(1.21 КБ) 1122 скачивания

[Voltron]

Каждый питон.. тьфу, каждый кулик своё болото хвалит.
<------------------------- skip ---------------------------------------->
Честно, не вижу, что здесь муторного. Фильтруем вывод gdalinfo (grep), разделяем на части и выводим только 3-ю и 4-ю (cut), убираем запятую (sed). Всё в одну строку.

Да я собственно и не спорю, что на bash тоже можно и ничего сложного там нет (не windows же). Просто явно размеры нужны не просто так, возможно, и другие параметры понадобятся и хитрая обработка. Потому и предложил

Можно увидеть такой же вывод размера на питоне? Статью читал, но прошу понимающих помочь (надо учиться).

А как же bash и «свое болото»?
Собственно вот

Код: Выделить всё

from osgeo import gdal

ds = gdal.Open( "/path/to/raster.tiff" )

print "Raster size", ds.RasterXSize, "x", ds.RasterYSize

Результат

Код: Выделить всё

Raster size 8251 x 7501

[Александр Мурый]

Просто явно размеры нужны не просто так, возможно, и другие параметры понадобятся и хитрая обработка.

"Подрихтовать" у топокарт информацию о СК и обрезать по вектору -- это явно не тот случай. Вот когда несколько каналов и т.п. -- то да, питон + GDAL.

А как же bash и «свое болото»?

Надо учиться питону, он сейчас везде + всё-таки полноценный ЯП. А от баша никуда не деться (ну разве что только пересесть на zsh).

Результат

Код: Выделить всё

Raster size 8251 x 7501

Отлично, а как получить чистые размеры "8251 7501" ?

[Александр Мурый]

..а как получить чистые размеры "8251 7501" ?

Туплю чего-то.. Конечно же:

Код: Выделить всё

# -*- coding: utf-8 -*-              

#!/usr/bin/env python

                                   

from osgeo import gdal

ds = gdal.Open( "/path/to/raster.tif" )

print ds.RasterXSize, ds.RasterYSize

[mak-vas]

* amuriy
Спасибо большое за gdal_proj_check!

[Александр Мурый]

* mak-vas
На здоровье
Главное, что без всяких ухищрений с получением размеров растра можно с помощью gdalwarp -r near получить то же самое --- правильно привязанный растр.

[mak-vas]

Решил свою проблему следующим циклом. Оно мне и нужно было. Потом подразберусь с башем и вмонтирую это в сценарий amuriy (надеюсь он по открытой лицензии распространяется? ):

Код: Выделить всё

X="100k--p42-107_108
100k--p42-119_120"
for x in $X
do
gdalwarp -s_srs EPSG:2502 -t_srs EPSG:28413 -dstalpha -cutline $x'.shp' -of GTiff -co "COMPRESS=LZW" $x'.tif' $x'_tr.tif'
done

У меня стояло задание — скачанные с "поехали" карты трансформировать в геотиф, обрезать по векторным маскам и преобразовывать в нужную проекцию. Естественно с минимальным участием человека.
Вроде бы разобрался со всем, остался последний штрих. Как вытащить из геотифа информацию о коде EPSG? Какая строка какой команды содержит его в себе? Я хочу, чтобы в вышеприведённом цикле проекция источника назначалась сама, без моего участия. gdalinfo не даёт этой информации в явном виде (или я не прав?).

[Александр Мурый]

Потом подразберусь с башем и вмонтирую это в сценарий amuriy (надеюсь он по открытой лицензии распространяется? ):

Ну да, типа GPL и всё такое.

Как вытащить из геотифа информацию о коде EPSG?[/b] Какая строка какой команды содержит его в себе? Я хочу, чтобы в вышеприведённом цикле проекция источника назначалась сама, без моего участия. gdalinfo не даёт этой информации в явном виде (или я не прав?).

А зачем точно указывать исходную проекцию? gdalwarp читает входную проекцию из файла (если она не пропущена там), можно указывать только выходную.

[mak-vas]

Действительно так и есть, благодарю!