GDAL+MapInfo: разрезаем изображение прямоугольной сеткой

Постановка задачи

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

Для решения задачи будет использоваться следующее ПО: утилиты GDAL 1.9.1 версии и MapInfo актуальной (9.5) версии. Ниже в командных инструкциях используется синтаксис командной строки Windows XP.

Полигональная сетка

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

Небольшое отступление

Самый известный вариант сетки — это сетка планшетов СК-42 (и аналогичных). Она стандартна, удобно создана для разных масштабов, и имеется у всех, кому нужна. С ней есть только один небольшой минус — она не прямоугольна. А смысл в том, чтобы резать изображение, не нарушая его пиксельной структуры — без пересчетов и искажений. Утилита gdalwarp позволяет обрезать изображение векторной сеткой, записанной в файл. Это не наш случай, т.к. он производит пересчет пикселей под векторный контур. Идея нашей обработки - извлечение прямоугольных участков растра, заданных двумя вершинами. Для этого у gdal_translate есть два механизма - задать координаты растра и задать координаты в проекции растра, а программа сама пересчитает эти координаты в координаты растра (и выведет их нам при работе).

Создание контура изображения

Для начала создадим внешний контур изображения. Тут есть два способа — простой и сложный, выбор зависит от свойств конкретных данных и частностей задачи.

Простой способ

Простая рамка изображения создаётся утилитой gdaltindex.

Например, для всех растров в папках, начинающихся с «MRGN_»:

for /d %b in (MRGN_*) do @for %a in (%b\*.tif) do if exist "%~fa" gdaltindex -skip_different_projection MRGN_SP5_2011_index.shp "%~fa"

для одной папки с растрами:

gdaltindex -skip_different_projection index.shp "*.tif"

Сложный способ

Нужно учитывать, что значимая часть растра может занимать только часть матрицы, оставшаяся часть растра нас не интересует, т.к. представляет собой «черное» поле (NODATA).

1. Cоздадим новый растр из предположения, что есть NODATA=0 и все остальное > 0. Используем утилиту gdal_calc, чтобы получить ненулевые значения (в примере обработка целой папки с растрами в фрмате TIF):

for /r %a in (*.tif) do gdal_calc -A "%~fa" --outfile="%~dpna_c.tif" --calc="(A>0)*127+(A==0)*0"  --NoDataValue=0 --type=Byte –overwrite

2. Сжимаем полученные файлы с помощью опции COMPRESS=CCITTFAX4 с очевидным уменьшением разрядности до 1 бита через gdal_translate:

for /d %b in (MRGN_*) do @for %a in (%b\%~b_c.tif) do gdal_translate -of "GTiff" -a_nodata 127 -co "NBITS=1" -co "TFW=YES" -co "COMPRESS=CCITTFAX4" "%~fa" "%~dpna_msk.tif"

Это делаем для того, чтобы однозначно избавиться от всех значений больше 0, и радикально уменьшить размер растра, если мы задумаем его хранить как маску.

Создаем контура значимых (непустых) частей снимков через утилиту векторизации gdal_polygonize:

for /d %b in (MRGN_*) do @for %a in (%b\*_c_msk.tif) do if exist "%~fa" gdal_polygonize "%~fa" -b 1 -f "ESRI Shapefile" %b\%~na_poly.shp

Теперь надо открыть полученный вектор в какой-то ГИС. MapInfo часто не может открыть векторные данные, получившиеся в результате работы GDAL-утилит, т.к. в них есть тестовые поля длиной 255 символов (у MapInfo ограничение длины строки в 254), так что контур можно перевести программно через ogr2ogr:

ogr2ogr -f "MapInfo File" MRGN_SP5_2011_poly_mi.tab MRGN_SP5_2011_poly.shp

Темы объединения векторов касаться не будем, это тоже можно сделать с помощью ogr2ogr.

Откроем контур (-а) растра (-ов) в MapInfo. С вероятностью 25%, в случае когда мы векторизовали значимый контур изображения через gdal_polygonize, в векторе будут присутствовать мелкие дырки от пикселей, которые имели значение NODATA внутри значимой части изображения.

• Для борьбы с ними в MapInfo разгруппируем все составные объекты, превращая «дырки» в отдельные объекты.
• После этого удалим все объекты, имеющие площадь меньше некоторой критичной (например, 1 кв.км). Должны остаться только контура изображений.

Теперь откроем контуры планшетов в ГИС. Проделаем с ними следующие превращения (большинство избыточны):

• Увеличим размер каждого планшета на некоторую величину (не обязательно, но если потом придется выводить растр с перекрытием, это даст запас пикселей по краям). Я увеличиваю размера планшета приблизительно на 5%, для простоты создания буфера я перед созданием буфера пересчитываю их в метры.
• Переводим полученный вектор в проекцию (или систему координат) растра.
• Новый растр открываем и находим для каждого контура MBR (minimum bounding rectangle). В MapInfo это стандартная функция, которую можно применить к объектам слоя, затирая их (это команда на MapBasic, которая выполняется в окне MapBasic внутри MapInfo):

update cutSHP_m50_000_UTM37 set obj=MBR(buffer(obj,20,2000,"m"))

что буквально означает следующее: "Для всех объектов слоя cutSHP_m50_000_UTM37 создать MBR для буфера вокруг каждого объекта с приращением 2000 метров и числом вершин при моделировании дуги равным 20".

• Надо расширить атрибутику векторных данных, т.к. в дальнейшем нам потребуются числовые значения координат вершин планшетов (под "планшетами" здесь понимаются «геометрические объекты сетки»). Этот пункт обязателен для исполнения (ниже станет ясно, почему):

Alter Table "cutSHP_m50_000_UTM37" ( add ulX float, ulY float, lrX float, lrY float)

• Теперь надо физически удалить планшеты, которые не пересекают значимый контур изображения. В MapInfo отбор таких объектов делается с помощь пространственного запроса. В этом же пункте можно удалить и избыточные объекты сетки, т.к. мы их создали с перекрытием.
• Теперь существенный момент: надо обрезать планшеты внешней границей растра (той, что получена в результате работы gdal_polygonize). Это обязательно надо сделать, т.к. планшеты, которые будут выступать за физические границы растра, не будут обработаны gdal_translate. В MapInfo это делается стандартной командой «Erase outside».
• Последний пункт – вычисление координат углов каждого планшета и сохранение этих координат в ранее созданные поля ulX float, ulY float, lrX float, lrY float. В MapInfo для этого существуют стандартные функции и команда Update:

update cutSHP_m50_000_UTM37 set ulx=ObjectGeography( object, 1 ),uly=ObjectGeography( object, 4 ),lrx=ObjectGeography( object, 3 ),lry=ObjectGeography( object, 2 )

• Внимание! Параметр -projwin у gdal_translate имеет вид

-projwin ulx uly lrx lry

Такая запись отличается от обычной, которая начинается с ЛЕВОГО НИЖНЕГО угла. Здесь, как в растре – начало в ЛЕВОМ ВЕРХНЕМ углу, а вторая пара координат соответствует ПРАВОМУ НИЖНЕМУ углу.

• Сохраняем векторные данные в формат MIF/MID так, чтобы они имели следующий вид:

Id, ulx, uly, lrx, lry

Id – уникальный номер (не обязательно, но полезно, чтобы задать упорядоченные имена для создаваемых файлов-фрагментов), X -координата верхнего левого угла, Y - координата верхнего левого угла, X координата нижнего правого угла, Y координата нижнего правого угла.

Это может быть не только формат MIF/MID, но и любой другой текстовый формат. В этом случае будем использовать только MID – текстовый файл атрибутивных данных. Все, что нам нужно – это координаты углов планшетов, чтобы извлекать подмножества из исходного растра.

Все «танцы с бубном» завершены.

Повторим, что нам необходимо было получить:

- 4 координаты 2-х углов сетки, которые лежат только в пределах обрезаемого (скорее разбираемого) растра и имеют ту же СК, что и обрезаемый растр. Если это у нас есть, то вперед - интерационо разберем растр на части.

Обрезка растра

Само использование GDAL, а именно, gdal_translate, гораздо короче:

1. Создаем набор команд для вычисления контуров в виде BAT-файла.
2. Команда выполняется в каталоге, где лежит файл Q1.MID (так я назвал экспорт координат углов), файлы выполнения сохраняются в каталог "root\1\", исходный файл находится по адресу "f:\test\test.tif":

for /f "tokens=1-5 delims=," %a in (Q1.MID) do @echo @if not exist "root\1\test_%~a.tif" gdal_translate -of GTIFF -projwin %b %c %d %e -a_nodata 0 -co "TFW=YES" -co "compress=LZW" "f:\test\test.tif" "root\1\test_%~a.tif">>root\1\do_%~a.bat

1. Можно сразу же запускать команды на выполнение в цикле for командной строки DOS, но создание набора BAT-файлов позволяет воспользоваться параллельным вычислением фрагментов или вернуться к прерванному вычислению, если понадобиться.

Разберём командные инструкции, написаные выше. Команда

for /f "tokens=1-5 delims=," %a in (Q1.MID) do

означает "Построчно считать файл Q1.MID и выполнить над каждой строчкой операцию, записанную после do":

• Использовать в качестве разделителя параметров строки символ «,»
• Используя разделитель, получить 5 параметров – переменные от а до e
• %а – уникальный номер, %b %c %d %e - координаты углов планшета из файла Q1.MID

Если не существует файл с именем "root\1\test_%~a.tif" , то:

• Выполнить команду gdal_translate
• над исходным файлом "f:\test\test.tif"
• с результирующим файлом "root\1\test_%~a.tif"
• -a_nodata 0 - условием "нет данных" (NODATA) для пикселей со значение 0 во всех каналах
• -of GTIFF - выходным форматом GeoTIFF
• -co "TFW=YES" -co "compress=LZW" с дополнительными параметрами «Создавать world file» и «Сжимать растр алгоритмом LZW»
• -projwin %b %c %d %e – вырезать фрагмент, заданный указанными координатами, в системе координат растра (gdal_translate переведет их в координаты растра перед тем, как вырезать фрагмент, о чем программа информирует нас по ходу выполнения). Если границы выйдут за пределы растра, программа завершится сбоем и сообщением об ошибке.

И, наконец, запускаем на выполнение все созданные BAT-файлы:

for %a in (do_*.bat) do "%~fa"