В качестве примера я возьму следующую плитку http://a.tile.openstreetmap.org/3/4/2.png и сохраню ее как «4_2.png».
Координаты WGS84 этой плитки можно вычислить или прочитать там , щелкнув соответствующую плитку:
0 66.51326044311185 45 40.97989806962013 (West North East South)
Как правильно привязать плитку (используя gdal для создания географического или другого формата с географической привязкой), чтобы:
- Растровое изображение не нужно растягивать (= пиксели в геотифе точно такие же, как в исходном растровом изображении)
- Полученное изображение будет открыто в нужном месте в средстве просмотра / редакторе ГИС (как, например, в TatukGIS Free Viewer )?
(Отредактировано 19 сентября 2011 г. для уточнения моего вопроса и включения моих выводов)
Мой вывод:
Сначала я понял, что третья идея (см. Ниже) была правильной. Я открыл геотиф в средстве просмотра ГИС и сравнил отображаемые координаты с ожидаемыми. Геотиф из второй идеи, кажется, сдвинут на 2 пикселя к северу. Вот почему я считал идею 3 (или 4) правильной.
Но если вы попытаетесь использовать плитку с гораздо более высоким уровнем масштабирования, геотиф из идеи 3 окончательно сместится на юг. Глупо было сравнивать координаты на плитке уровня масштабирования 3. Границы страны на таком уровне масштабирования были упрощены, так что сравнение не дает хороших результатов.
Дэн С. был прав, изображение плитки уже в EPSG: 3857. Тогда вторая идея является правильной (и дает хороший результат и при высоких уровнях масштабирования)
Первая идея: EPSG: 4326
Код EPSG для координат WGS84 - EPSG: 4326 . Поэтому я просто использую координаты WGS84, чтобы геозонить плитку как геотиф с помощью gdal_translate :
gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 66.51326044311185 45 40.97989806962013 -a_srs EPSG:4326 4_2.png t4326.tif
Получившаяся карта отображается в нужном месте, но я боюсь, что проекция неверна и что в середине плитки может произойти смещение. После долгих попыток проверить это путем перепроецирования карты с помощью gdalwarp, я скачал демо-версию Global Mapper, и это похоже на случай (те же границы, что и в идее 3, но смещение внутри фрагмента). Изображение должно быть растянуто, чтобы можно было использовать EPSG: 4326 координат.
Вторая идея: EPSG: 3857
Эта плитка использует проекцию «web mercator» (псевдоним Google Map проекция), которая теперь имеет код EPSG: EPSG: 3857 (псевдоним EPSG: 900913). Я просто конвертирую координаты, используя gdaltransform :
gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.51326044311185
0 10018754.1713946 0
45 40.97989806962013
5009377.08569731 5009377.08569731 0
Мои координаты в метрах:
0 10018754.1713946 5009377.08569731 5009377.08569731 (West North East South)
Теперь я могу использовать gdal_translate для генерации геотиффа:
gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 10018754.1713946 5009377.08569731 5009377.08569731 -a_srs EPSG:3857 4_2.png t3857.tif
У меня сложилось впечатление, что это не правильно, потому что границы карт смещены к северу. Кажется, это правильная идея.
Третья идея: EPSG: 3857 - EPSG: 4055
Я читал, что «веб-меркатор» использует координаты WGS84, но рассматривают их так, как если бы они были сферическими координатами. Из-за разницы между геодезической и геоцентрической широтой (см. Википедию о широте ), значения широты не будут одинаковыми на эллипсоиде или на сфере. Я обнаружил, что EPSG: 4055 - это код сферических координат на сфере, основанный на WGS84.
Преобразование координат в EPSG: 4055:
gdaltransform -s_srs EPSG:4326 -t_srs EPSG:4055
0 66.51326044311185
0 66.3722684317026 -17964.0621483233
45 40.97989806962013
45 40.7894557844857 -9152.84527519904
Соответствующие сферические координаты тогда:
0 66.3722684317026 45 40.7894557844857 (West North East South)
Затем я делаю, как будто те координаты, которые все еще на эллипсоиде (EPSG: 4326) и преобразовываю их в веб-меркаторе:
gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857
0 66.3722684317026
0 9979483.26733298 0
45 40.7894557844857
5009377.08569731 4981335.86590183 0
Результирующие координаты отличаются от координат по идее 2:
0 9979483.26733298 5009377.08569731 4981335.86590183 (West North East South)
Теперь мне просто нужно записать координаты на карту:
gdal_translate -of Gtiff -co tfw=yes -a_ullr 0 9979483.26733298 5009377.08569731 4981335.86590183 -a_srs EPSG:3857 4_2.png t3857_through_4055.tif
Эта третья идея, кажется, дает лучшие результаты. Но я не уверен, правильно ли это. Если идея 3 верна, существует ли код EPSG для выполнения этой операции за один шаг?
Четвертая идея: EPSG: 3857 через towgs84 = 0,0,0,0,0,0,0
gdal (и, видимо, тоже epsg) определяет EPSG: 3857 следующим образом:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
тогда как пространственные ссылки.org вот так:
+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Если я использую определение изatialreference.org, я получу правильные координаты за один шаг (ну, я все еще не знаю, являются ли они «правильными» координатами, но, по крайней мере, они такие же, как в идее 3):
gdaltransform -s_srs EPSG:4326 -t_srs "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
0 66.51326044311185
0 9979483.26733298 -17964.0621483233
45 40.97989806962013
5009377.08569731 4981335.86590183 -9152.84527519904
Почему такая разница в определениях EPSG: 3857?
источник
Функции, необходимые для генерации глобальных плиток, используемых в сети. Он содержит классы, реализующие преобразование координат для:
GlobalMercator (на основе EPSG: 900913 = EPSG: 3785 ) для Google Maps, Yahoo Maps, Microsoft Bing Maps совместимых плиток
GlobalGeodetic (на основе EPSG: 4326) для базовой карты OpenLayers и плиток, совместимых с Google Планета Земля
http://svn.osgeo.org/gdal/sandbox/klokan/gdal2tiles.py
источник
О моем дополнительном вопросе в четвертых идеях: почему существует такая разница в определениях EPSG: 3857 между определением в gdal и на сайтеrereference.org :
Основное отличие состоит в том, что gdal использует «+ nadgrids = @ null» и пространственную привязку «+ towgs84 = 0,0,0,0,0,0,0». Согласно документации или формату PROJ.4 оба параметра используются для базовых преобразований. Я нашел интересный комментарий от Микаэля Риттри на сервере рассылки Proj4 :
Использование «+ towgs84 = 0,0,0,0,0,0,0» не представляется правильным или, по крайней мере, только при определенных условиях.
источник