Тема не вопрос, а просто делюсь опытом.
Естественно вопрос избитый, тысячу раз обмусоленный, но по опыту знаю, что иногда надо 100 прочитать об одном и том же из разных источников, чтобы осмыслить все до конца.
Итак была задача сформировать тайловый массив в файловой системе в формате .../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
)
Код: Выделить всё
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"
В Type"Raster" нужен пробел... Type "Raster"