GIS-Lab - Вклад [ru]

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

2015-03-20T12:04:04Z

SergeyRyzhkov:

Итак была задача сформировать тайловый массив в файловой системе в формате .../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"

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

2015-03-20T12:00:42Z

SergeyRyzhkov: Новая страница: «Итак была задача сформировать тайловый массив в файловой системе в формате .../Z/X/Y.png (не TMS…»

Итак была задача сформировать тайловый массив в файловой системе в формате .../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"