gdal_translate: привязанный растр из тайлового источника

[Denis Rykov]

Потребовалось на основе тайловых данных сервиса openstreetmap.org подготовить набор геопривязанных растров, представляющих собой территории определённых городов, отрендеренные на нужном масштабном уровне. Решил попробовать выполнить данную задачу при помощи консольных утилит GDAL. В списке форматов, поддерживаемых GDAL-ом, находим OGC Web Map Service и чуть ниже на этой же страницы видим описание мини-драйвера TMS, отвечающего за поддержку тайловых источников данных. Как им пользоваться смотрим в разделе Examples. Нас интересует пример OpenStreetMap TMS Service Example:

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

gdal_translate -of PNG -outsize 512 512 frmt_wms_openstreetmap_tms.xml openstreetmap.png

Как видно из представленного примера, вначале мы должны представить описание нашего тайлового сервиса в виде XML-файла и затем скормить его утилите gdal_translate. Создаём файл osm.xml и помещаем в него следующий код:

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

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:900913</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache />
</GDAL_WMS>

Теперь попробуем получить растр на территорию села Фершампенуаз, имеющего охват 6656381.223504174, 7076012.868322925, 6662582.175551272, 7082935.574421264 в проекции тайлов 900913:

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

gdal_translate -of PNG -projwin 6656381.223504174 7082935.574421264 6662582.175551272 7076012.868322925 -outsize 512 512 osm.xml fershampenuaz512.png

Получили растр.

fershampenuaz512.png (111.89 КБ) 12054 просмотра

Но не понятно на основе тайлов какого уровня он построен?? Предполагаю, что GDAL должен сначала по указанному bbox-у и размеру выходного изображения рассчитать разрешение входного растра и определить к разрешению какого масштабного уровня оно более близко, то есть вычислить координату тайлов z, после этого запросить нужные тайлы, вырезать кусок и вместить в квадрат указанного размера (512px), но не уверен. В любом случае при таком способе (явном указании размера выходного изображения в пикселах) на выходе мы будем почти всегда иметь деформированный по одной из осей растр. Попробуем задать размер выходного растра в процентах:

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

gdal_translate -of PNG -projwin 6656381.223504174 7082935.574421264 6662582.175551272 7076012.868322925 -outsize 10% 10% osm.xml fershampenuaz10pct.png

fershampenuaz10pct.png (333.34 КБ) 12054 просмотра

Получившийся растр имеет нужные пропорции. Но возникает 2 вопроса. Первый - что за проценты? Второй - какой масштабный уровень при этом используется?

P.S. Что касается привязки, то для этого достаточно указать -of GTiff вместо -of PNG.
P.P.S Файл osm.xml можно открыть в QGIS как растровый слой.

[Дмитрий Барышников]

Что бы получить четко без искажений нужный уровень, необходимо задать размеры выходного растра, таким образом, что бы не происходило преобразований, т.е. пикселы копировались 1 : 1.
При указании размера или процентов - это все от исходного размера изображения (в случае OSM оно виртуальное, но легко рассчитывается через bbox = 67108864 х 67108864 пикселов). Т.е выходная картинка соответствует некоему промежуточному уровню.
Программно для этого имеются методы GDALBandGetBestOverviewLevel и GDALDatasetGetBestOverviewLevel (метод внутренний для нужд GDAL). Причем, после вызова такого метода, вы получаете не только необходимый уровень пирамид, но и размеры выходного растра и смещение внутри пирамиды.

Например:
1. На входе растр или его участок размером 25130 х 15438.
2. На ходе хотим 866 х 532
3. После вызова метода для первого слоя получим: уровень пирамид 4, выходной размер 785 х 482.

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

Резюмирую, по-моему для решения вышеозначенной задачи только кодить.
Можно попробовать завернуть все в VRT со сглаживанием при помощи AveragedSource.
Может быть еще есть идеи как без программирования решить?

[Denis Rykov]

Спасибо, немного прояснилось, но всё равно остались вопросы. Итак рассмотрим на конкретном примере, что мы имеем.

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

gdal_translate -of PNG -projwin 6656381.223504174 7082935.574421264 6662582.175551272 7076012.868322925 -outsize 10% 10% osm.xml fershampenuaz10pct.png

1. Вычисляем размеры запрашиваемого изображения в проекции тайлов 900913:

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

   x = 6662582-6656381=6201(м);
   y = 7082935-7076012=6923(м);

2. Вычислим сколько пикселов займёт данное изображение на максимальном 18 масштабном уровне (0.5972 - разрешение на данном уровне):

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

   X = x/0.5972 = 6201/0.5972 = 10383;
   Y = y/0.5972 = 6923/0.5972 =  11592;

3. Так как мы указали размер выходного растра 10% на 10%, то получаем размер в пикселах: 1038x1159, что соответствует действительности (можно посмотреть размер растра из первого сообщения в теме).
4. Определим разрешение выходного растра:

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

   Rx = x/X = 6201/1038 = 6;
   Ry = y/Y = 6923/1159 = 6

   Получается, что мы попали между двумя уровнями - 14(9.5546) и 15(4.7773). Какой-то из них (какой?) берет GDAL и формирует на их основе выходное изображение (растягивая или сжимая до вычисленных размеров).

Таким образом, чтобы получить выходное изображение, построенное на базе определённого масштабного уровня, нужно варьировать размером выходного растра. Допустим, мы хотим получит растр на базе тайлов 10 уровня (разрешение 152.87), тогда размер растра можно вычислить так:

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

X = x/R = 6201/152.87 = 40;
Y = y/R = 6923/152.87 = 45;

И команда будет выглядеть следующим образом:

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

gdal_translate -of PNG -projwin 6656381.223504174 7082935.574421264 6662582.175551272 7076012.868322925 -outsize 40 45 osm.xml fershampenuaz40x45.png

Если я правильно понял, то ничего программировать не нужно, нужно всего лишь посчитать размеры выходного растра, а для этого достаточно разделить X и Y на разрешение нужного Z-уровня или нет?

[Дмитрий Барышников]

Какой-то из них (какой?) берет GDAL и формирует на их основе выходное изображение (растягивая или сжимая до вычисленных размеров).

1. Уровень выбирается ближайший к необходимому размеру выходного растра. В данном случае должен был быть 15.

Если я правильно понял, то ничего программировать не нужно, нужно всего лишь посчитать размеры выходного растра, а для этого достаточно разделить X и Y на разрешение нужного Z-уровня или нет?

2. В принципе да, просто уже не одной командой. Да еще нужно знать размер пиксела для растра. Кстати у меня gdalinfo дает 0.597164 х 0.597164.

P.S. Интересно, как сочетается минидрайвер для Google с командой gdal_translate - это же нарушение лицензии? Или нет? Так и Sas.Planet не нужен. Сразу выходная картинка с привязкой в нужном формате.

[Kosoku]

Добрый день, коллеги!

Изначально писал в ЛС Денису свои вопросы, но решил вынести на общее обсуждение свой провал.

Задача стоит следующая: через возможности GDAL (gdal_translate и, возможно, gdaltransform) выдергивать с доступного web ресурса (arcgis, google, openstreet... ) растр, центрованный на вводимых пользователем координатах.
Находясь в Петербурге, воюю с попаданием в эти координаты: 59.95, 30.316667

Знакомство с GDAL и web картографией началось с http://gdal.org/frmt_wms.html
Отсюда был взят текст .xml файла и примерно формировался запрос, при этом разброс между координатами верхнего левого и правого нижнего углов необходимого растра давались как +/- 5 градусов от координат Петербурга.
Попасть в указанные координаты не получалось, пока не начали подгонять диапазон Y в исполняемом .xml файле. Когда задали диапазон от -143 до 143 градусов, и в запросе gdal_translate... указали необходимые координаты в градусах, попали в Питер, но изображение было искажено по горизонтали (растянуто).
Долго пытался понять в чем причина, искал формулу изменений диапазонов и расстояний при удалении от экватора и так далее.
Сейчас уже подкипает мозг после этого оказался чудом на вашем форуме, где написал сообщение Денису, который (за что ему огромное спасибо), дал ссылку на эту тему.

Скопировал его шаги (взял текст .xml документа и добавил свои координаты в его запрос командной строки).
Координаты для тестов получал след. образом:
1. Установил ArcGIS Explorer
2. Выделил в нем желаемую область вокруг Петербурга
3. Открыл инфо о выделенной области и взял указанные там координаты (они находились вокруг указанных выше координат Петербурга, потому в них нет сомнения)
4. Через оболочку, поставляющуюся с GDAL запустил gdaltransform где указал трансформацию из 4326 в 900913
5. Полученные после трансформации данные добавил в запрос Дениса gdal_translate и отправил на выполнение.

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

В чем моя ошибка? Куда копать, где смотреть теорию?

[Kosoku]

В добавление к предыдущему посту.

Мой запрос выглядит так:
gdal_translate -of PNG -projwin 3412110.7941 6714546.7817 3650397.9450 6625457.7932 -outsize 10% 10% gis.xml saintP.png

[Kosoku]

Нашел свой же косяк в отсылаемых координатах, переделал запрос и получил центр СПб.
Запрос для Петербурга для .xml Дениса:
gdal_translate -of PNG -projwin 3322980.57 8438614.34 3422980.57 8338614.34 -outsize 10% 10% gis.xml saintP.png

Следующий вопрос:
Получаю адское разрешение, которое не обязательно. Как его уменьшить? Сократить количество тайлов в .xml документе в строке <TileLevel>18</TileLevel>?

[Denis Rykov]

В зависимости от задаваемого вами размера выходного растра используется тот или иной масштабный уровень, вот проверил на ваших данных, вроде работает (выходной PNG получился 1.6 Mb с сохраненными пропорциями):

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

rykov@extensa:~/temp/gdal$ python
Python 2.7.3 (default, Apr 10 2013, 05:46:21) 
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> zoom_level_11 = 76.437
>>> dx = 3422980.57 - 3322980.57
>>> dy = 8438614.34 - 8338614.34
>>> (x,y) = dx/zoom_level_11, dy/zoom_level_11
>>> x,y
(1308.26693878619, 1308.26693878619)
>>> 
rykov@extensa:~/temp/gdal$ gdal_translate -of PNG -projwin 3322980.57 8438614.34 3422980.57 8338614.34 -outsize 1308 1308 osm.xml saintP.png
Input file size is 67108864, 67108864
Computed -srcwin 39119032 19423288 167458 167458 from projected window.
0...10...20...30...40...50...60...70...80...90...100 - done.

[Kosoku]

Денис, спасибо!

Просто я указывал, как у вас в примере, -outsize 10% 10%, и потому получил картинку на 111 Мб с уровнем детализации соответствующим.
Завтра буду штурмовать варианты с растром соответствующего размера. Уровень детализации задается в соответствии с таблицей размеров тайлов?

[Denis Rykov]

Уровень детализации выбирается вами самостоятельно, для tms-профиля global-mercator это стандартные значения, можно посмотреть например тут.