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=0). В этом случае, я делаю еще «пару танцев с бубном»:

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"