Нарезка тайлов из растров MrSID c *.tab привязкой

Эта страница является черновиком статьи.

Итак была задача сформировать тайловый массив в файловой системе в формате .../Z/X/Y.png (не TMS!) Исходными данными был набор растров в формате *.sid (MrSID) с фалами привязки в формате *.tab (MapInfo). При этом система координат отсутствует, вернее - CoordSys Nonearth Units "m" , сами координаты в МКС-64

Сразу ринулся естественно к GDAL. Не имея достаточного опыта работы с GDAL, в голове изначально был следующий алгоритм (не совсем неправильный): 1. Преобразовать все *.sid в GeoTIFF 2. Объединить все полученные на п1 в единый 4. Для единого растра назначить СК, а затем перепроецировать в EPSG:3857 5. Нарезать тайлы с использованием gdal2tiles

Не получилось, а именно 1. Файлы привязки *.tab не "схватывались" 2. Когда разобрался с п1, объединение шло очень долго 3. После перпроецирования получился артифакт, небольшой поворот и черная неравномерная окантовка вокруг растра 4. gdal2tiles "генерит" согласно спецификации TMS. То есть направление оси Y с юга на север

Не буду долго описывать ход экспериментов и т.д., сразу напишу как я сделал (но тоже не факт что оптимально) Примеры кода будут под винду, bat-команды (я образно)

1. GDAL при конвертации и т.д. использует геопривязку. При этом если геопривязка есть в самом файле, то внешние файлы привязки, например *.tab игнорируются Решение: Убираем геопривязку из файлов *.sid

Код

   rem Обязательно с параметром -co "PROFILE=BASELINE",
   rem Иначе геопривзяка из файла не будет удалена и дальнейшие команды не будут использовать фалы привязки
   for %%f in (*.sid) do (
      gdal_translate -co "PROFILE=BASELINE" -q %%f %%~nf.tif
   )

2. GDAL, как и, например, QGIS, ArcGIS, прекрасно работаю с так называемыми виртуальными наборами растров (VRT), которые представляют из себя просто XML. То есть нет необходимости мержить (объединять) все растры в один, чтобы потом использовать для других операций.

Код

   rem Построение виртуального набора растров для всех полученных *.tif
   rem Назначение системы координат с геопривязкой "внтури" файла
   rem В данном случае уже будут использоваться файлы геопривязки
   gdalbuildvrt -overwrite -q -a_srs "+proj=tmerc +lat_0=0 +lon_0=30 +k=1 +x_0=... +y_0=... +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs" -vrtnodata "255 255 255" -hidenodata -resolution highest MergedMsk.vrt *.tif

При этом устанавливаем параметры для сокрытия пустых областей

3. Дабы избежать артефакта черной каймы в результирующем растре, оказывается можно не делать перепроицирование в в EPSG:3857, а сразу отдать виртуальный набор растров скрипту gdal2tiles То есть вот так делать не обязательно, и не надо

Код

   rem Перепроецирование виртуального набора растров
   rem gdalwarp -of VRT -overwrite -t_srs EPSG:3857 MergedMsk.vrt Merged3857.vrt 

Сразу "генерим" растры модифицированным питон-скриптом:

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

   rem Создание тайлового кэша. Используется модифицированный код, так как стандартный в схеме TMS
   python C:\_Work\_Projects\nRsn\CoreLibs\gdal-1.11.2\bin\gdal\python\scripts\gdal2tilesXYZ.py -z 11-18 MergedMsk.vrt pyramid

Что я дополнительно сделал:

Код:

   rem Далее можно прогнать OptiPNG и PNGOUT
   rem Можно из проводника http://blog.sijey.ru/2011/04/04/optipng_pngout_integration/

Также можно удалить все пустые тайлы (очень очень большое количество). Определить их легко, их размер 4 096 байт.

Как результат, всего за 1 сутки я построил тайловый кэш ортофотопланов 11-18 уровня на территорию СПб.

UPD: Относительно долго не мог понять почему GDAL не принимает файл привязки после "очистки" файлов от привязки. Оказалось поставщик растров сделал *.tab в таком формате

Код:

   !table
   !version 300
   !charset WindowsLatin1
   Definition table
   File   "2222-04.sid"
   Type"Raster"
   (87000,88000)(0,0) Label "Point 1",
   (88000,88000)(1666,0) Label "Point 2",
   (87000,87000)(0,1666) Label "Point 3",
   (88000,87000)(1666,1666) Label "Point 4"
   CoordSys Nonearth Units "m"

Для MapInfo не критично, а вот GDAL-у не понятно. Оказалось банально (но я долго не мог понять): В Type"Raster" нужен пробел... Type "Raster"