Привязка растровой карты Арктики

[MYS]

Здесь вам вот что пытаются объяснить.

Нельзя привязать карту, не используя информацию о ее проекции. То, что вам удалось, по вашим словам, это сделать раньше - удачное стечение обстоятельств. Например, карта была в географических координатах, которые вы ей и предлагали, давая точки привязки.

Ваша карта Арктики использует метры для измерения координат, привязывать ее надо, давая точкам координаты в метрах (см. пост Александр Мурый выше), и эти координаты должны соответствовать конкретной проекции. То, что на карте нарисована градусная сетка, ничего не значит, и использовать ее для привязки, давая координаты в градусах, нельзя.

Ну, и конечно очень интригует ваше высказывание в первом посте: "Программный проект не использует никаких библиотек или компонент, специфических именно для ГИС" - это, чтобы легких путей не искать?

[Ariki]

Вы хотите сказать, что надо сначала пиксели переводить в координаты полярной стереографической проекции, а из нее - в градусы широты-долготы ?

Да, с единственной поправкой, что в данном случае проекция всё-таки полярная равнопромежуточная. В этом легко убедиться, измерив линейкой расстояния между параллелями на карте. Формулы проекции можно найти в Википедии.

Картографическая проекция - это функция, отображающая географические координаты на плоскость. У разных карт она разная, поэтому, если вы хотите поддерживать разные типы проекций, рекомендую всё-таки использовать библиотеку PROJ.4, поддерживающую десятки разных проекций. Далеко не все из них можно аппроксимировать полиномами, как вам уже намекнули, в том числе и из-за топологических разрывов. Вообще аппроксимацию допустимо использовать только в случаях, когда нет другого выхода, при пересчёте между схожими проекциями, и то на ограниченных территориях, а никак не для целого полушария.

А вот когда вы уже пересчитали координаты из градусов в единицы проекции, вы можете использовать полиномы, чтобы пересчитать эти координаты в пиксели, заодно убрав искажения отсканированного материала. Хотя обычно хватает аффинного преобразования или преобразования подобия. Коэффициенты такого преобразования можно хранить в world-файле или метаданных самого растра - почитайте статьи на ГИС-Лабе.

Для работы с геопривязанными растрами есть замечательная библиотека GDAL.

[KotAlex]

Большое спасибо участникам MYS и Ariki, от которых я наконец-то получил вменяемые ответы. Неужели все остальные не могут ответить так же ?

Нельзя привязать карту, не используя информацию о ее проекции. То, что вам удалось, по вашим словам, это сделать раньше - удачное стечение обстоятельств. Например, карта была в географических координатах, которые вы ей и предлагали, давая точки привязки.

Значит, если сама карта - в географических координатах и точкам привязки я тоже задаю географические координаты, то вопрос о проекциях карты отпадает ? Использовать удачное стечение обстоятельств - в этом, по-моему, ничего незаконного нет. Или, по-вашему, такой способ все-таки некорректен ?

Ну, и конечно очень интригует ваше высказывание в первом посте: "Программный проект не использует никаких библиотек или компонент, специфических именно для ГИС" - это, чтобы легких путей не искать?

Реализованное мною приложение вводит файл изображения карты на стандартную канву, на которой делаются стандартные графические операции - в том числе и с точками привязки. Вводится-записывается стандартный текстовый файл привязки точек. Модуль полиномиального преобразования я целиком взял из источника, приведенного в головном сообщении. Что тут требуется такого специфического из области ГИС ? А во-вторых, мне нужно отрабатывать новые технологии, для которых заведомо нет реализации в стандартных библиотеках ГИС.

[KotAlex]

Картографическая проекция - это функция, отображающая географические координаты на плоскость. У разных карт она разная,

Хорошо, если я правильно понимаю, что такое проекция, то для приведенной мною в головном сообщении карты Арктики проекция выглядит так : сама проекция - в плоскости изображения карты, начало координат проекции - левый нижний угол карты, оси направлены : ось Y - строго вверх ( совпадает с левой границей карты ), ось X - строго горизонтально ( совпадает с нижней границей карты ). И единицы измерений по осям - метры ( или что там еще ... ) ?

...десятки разных проекций. Далеко не все из них можно аппроксимировать полиномами, как вам уже намекнули, в том числе и из-за топологических разрывов.

То, что на карте Арктики будет разрыв и использовать тот же метод привязки не удастся - это мне было понятно с самого начала. Потому и задал вопрос. Понятно было, что здесь надо будет изгиляться, но хотелось, все же, найти более простое решение в рамках прежнего.

[KotAlex]

https://ru.wikipedia.org/wiki/%D0%A2%D0 ... 0%B8%D1%8F
разрежь картинку пополам

Не пойдет - картинка нужна целая.

[trir]

нет, центр СК - на полюсе
[offtop]почему люди любят изобретать велосепеды с квадратными колёсами?[/offtop]

[KotAlex]

нет, центр СК - на полюсе
[offtop]почему люди любят изобретать велосепеды с квадратными колёсами?[/offtop]

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

По поводу СК - позвольте уточнить : метрика по обеим осям - строго равномерная вдоль оси и одинаковая для обеих осей ?

[Ariki]

Единицы проекции - это просто декартовы координаты на плоскости, на которую спроецирована поверхность эллипсоида. Они могут соответствовать к метрам или другим единицам на местности в определённых точках карты, но в целом масштаб любой карты - величина непостоянная в пределах этой карты. Тот масштаб, что подписан на карте, является усреднённым либо относится к некоторой точке или линии на карте, имеющей нулевые искажения. Невозможно отобразить сферическую поверхность на плоскость, не исказив расстояния.

[KotAlex]

... в целом масштаб любой карты - величина непостоянная в пределах этой карты. Тот масштаб, что подписан на карте, является усреднённым либо относится к некоторой точке или линии на карте, имеющей нулевые искажения. Невозможно отобразить сферическую поверхность на плоскость, не исказив расстояния.

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

[gamm]

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

географические координаты - это долгота и широта. В этих координатах карты не делают, только глобусы (карты всегда в какой-то проекции), и тем более масштабов в долготе/широте не считают ... т.е. в реальности все обстоит с точностью до наоборот. Но читать вас интересно (в пятницу вечером - особенно)

[Denis Rykov]

Согласен. С небольшим уточнением : масштаб, подписанный на карте, относится к географическим координатам, в которых сделана карта.

У вас на карте даже масштабная линейка есть и там указаны единицы измерения - "км". Откуда вы взяли географические координаты?

[KotAlex]

У вас на карте даже масштабная линейка есть и там указаны единицы измерения - "км". Откуда вы взяли географические координаты?

Согласен. Я действительно лопухнулся. Непростительно. Сегодня у меня был очень тяжелый день.

[Denis Rykov]

Попробовал привязать ваш растр, вот что получилось:

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

$ gdal_translate -of GTiff -gcp 393.911 364.442 -2.15198e+06 -2.56463e+06 -gcp 394.431 1063.67 2.15198e+06 -2.56463e+06 -gcp 1228.15 1065.75 2.15198e+06 2.56463e+06 -gcp 1229.97 366.002 -2.15198e+06 2.56463e+06 "/home/rykov/tmp/arctic.jpg" "/tmp/arctic.jpg"
$ gdalwarp -t_srs "+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r cubicspline -order 1 -co COMPRESS=NONE "/tmp/arctic.jpg" "/home/rykov/tmp/arctic.tif"

Координаты в полярной равнопромежуточной проекции можно рассчитать используя PROJ.4 (утилита cs2cs, есть баиндинги в Python):

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

$ cs2cs +init=epsg:4326 +to +proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
40 60
2151984.08	-2564634.76 0.00
-40 60
-2151984.08	-2564634.76 0.00
140 60
2151984.08	2564634.76 0.00
-140 60
-2151984.08	2564634.76 0.00

[KotAlex]

...

Большое спасибо за проделанную работу, однако, из приведенной вами информации мало что понял. Пожалуйста, поясните :

1.

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

$ gdal_translate -of GTiff -gcp 393.911 364.442 -2.15198e+06 -2.56463e+06 -gcp 394.431 1063.67 2.15198e+06 -2.56463e+06 -gcp 1228.15 1065.75 2.15198e+06 2.56463e+06 -gcp 1229.97 366.002 -2.15198e+06 2.56463e+06 "/home/rykov/tmp/arctic.jpg" "/tmp/arctic.jpg"
$ gdalwarp -t_srs "+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -r cubicspline -order 1 -co COMPRESS=NONE "/tmp/arctic.jpg" "/home/rykov/tmp/arctic.tif"

Что здесь делается, для чего и что получаем ?

2.

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

$ cs2cs +init=epsg:4326 +to +proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
40 60
2151984.08	-2564634.76 0.00
-40 60
-2151984.08	-2564634.76 0.00
140 60
2151984.08	2564634.76 0.00
-140 60
-2151984.08	2564634.76 0.00

Здесь я понял так, что вы для пары координат каждой точки, заданных в долготе-широте, получаете пары координат в полярной равнопромежуточной проекции ? Верно или нет ?

3.
В приложенном файле записано такое же изображение карты, как у меня, повернутое на 90 градусов. Что это означает ? Для чего это изображение получено ?

[Denis Rykov]

1. С помощью утилиты gdal_translate мы задаём контрольные точки (сопоставляем пиксельные координаты выбранных точек и их спроецированные координаты), затем с помощью gdalwarp трансформируем исходный растр и назначаем ему выбранную систему координат. Подробнее можно прочитать тут. Ещё остались вопросы?

2. Всё правильно, есть ещё веб морда к cs2cs вот, пример расчета во вложении.

3. Во вложении привязанный GeoTIFF, то есть если вы загрузите его в какую-нибудь ГИС (например, QGIS), то перемещая указатель мыши будите видеть в статустной строке не пиксельные координаты, а реальные координаты в проекции растра (метры). Кроме того, в том же QGIS можно включить нужные canvas units (например, градусы) и он будет автоматом переводить координаты из системы координат растра в градусы и отображать и в статусной строке.