Аффинное преобразование координат MapInfo

[inouterer]

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

[Mitrich]

может поможет, только вчера искал, правда немного другое:
подробнее ссылка в этом посте
viewtopic.php?t=1494

помогите разобраться?

* Three shift parameters specifying the distance, in meters, to shift the ellipsoid along each of its axes. These parameters are usually denoted by dX, dY, and dZ. They may also be denoted by X, Y, and Z, or by u, v, and w.

* Three rotation parameters specifying the angle, in arc-seconds, to rotate the ellipsoid around each of its axes. These parameters are usually denoted by EX, EY, and EZ. They may also be denoted by X, Y, and Z, or by , , and .
* A scale correction factor specifying the amount, in parts per million, to adjust the size of the ellipsoid. This parameter is denoted by the letter m, or sometimes .

* The longitude of the prime meridian, in degrees east of Greenwich. The prime meridian specifies which location on earth is assigned longitude 0°. Most datums use Greenwich as the prime meridian, so this parameter is usually zero. However, some datums use a different location as the prime meridian. For example, the NTF datum uses Paris as its prime meridian, which is 2.33722917 degrees east of Greenwich. If the NTF datum is used in a coordinate system, all longitudes in that coordinate system are relative to Paris instead of Greenwich.


A custom datum can be defined in any coordinate system definition. The Appendix describes how to define a coordinate system. To define a custom datum in a coordinate system, use datum number 9999 followed by the datum parameters, in this order:

9999, EllipsoidNumber, dX, dY, dZ, EX, EY, EZ, m, PrimeMeridian

Some datums specify only an ellipsoid and shift parameters (dX, dY, dZ), with no rotation parameters, scale correction, or prime meridian. In those cases, use datum number 999 instead of 9999, to simplify the definition: (See Example LD-3 below):

999, EllipsoidNumber, dX, dY, dZ

[Донецков]

Есть ли простой способ нахождения параметров:

dX, dY, dZ, EX, EY, EZ, m, PrimeMeridian

по минимум 3 парам координат в известной и искомой системе?

[Игорь Белов]

Есть ли простой способ нахождения параметров:

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

dX, dY, dZ, EX, EY, EZ, m, PrimeMeridian

по минимум 3 парам координат в известной и искомой системе?

Утилита от Ракурса?

Теперь предположим, что «простой способ» — программирование, хотя бы в электронных таблицах. Вам нужно 8 параметров? Не так всё страшно.

Наборы координат должны быть в трёхмерной системе X, Y, Z ("Earth centered, Earth fixed").

PrimeMeridian должен быть известен заранее, поскольку он связан с вращением функционально. 8 − 1 = 7

Находим «центры масс» точек в известной и искомой системах координат и переносим начала координатных систем в эти центры. Параметры переноса нейтрализованы: 7 − 3 = 4

Находим суммы квадратов векторов точек в новых системах. Масштаб равен квадратному корню отношения сумм. Перемасштабируем одну из систем. Масштаб исключён: 4 − 1 = 3

Теперь самое трудное: для нахождения трёх неизвестных параметров вращения EX, EY, EZ задействуем метод наименьших квадратов.

Вернёмся к исходным системам. Используя «центры масс», вычислим параметры переноса dX, dY, dZ как разницу между имеющимися координатами и трансформированными из другой системы с нулевым переносом.

Пересчитаем углы из радиан в секунды, а масштаб в поправку к единице в миллионных долях. Всего и делов-то!

[Донецков]

!!! одной писанины почти на страницу...

[Игорь Белов]

Ерунда, каких-то несколько тысяч знаков в кодах.

[Донецков]

Если серьезно, то готового решния по этой теме похоже нема..., а если и есть, то платно или под сукном...

[Boris]

все куда-то ушли в сторону подбора нового эллипсоида, а вопрос был про афинное преобразование, задаваемое в операторе Coordsys. Если я угадал, то вот текст, объясняющий порядок:

The following example defines a coordinate system called DCS that is derived from UTM Zone 10 coordinate system using the affine transformation.

x1 = 1.57x - 0.21y + 84120.5
y1 = 0.19x + 2.81y - 20318.0
In this transformation, (x1, y1) represents the DCS derived coordinates, and (x, y) represents the UTM Zone 10 base coordinates. If the DCS coordinates are measured in feet, the CoordSys clause for DCS would be as follows:

CoordSys Earth
Projection 8, 74, "m", -123, 0, 0.9996, 500000, 0
Affine Units "ft", 1.57, -0.21, 84120.5, 0.19, 2.81, -20318.0

для теоретической основы: http://oldskola1.narod.ru/Jakovlev/Jakovlev013.htm

[Игорь Белов]

Если серьезно, то готового решния по этой теме похоже нема..., а если и есть, то платно или под сукном...

Собственно, программу я написал несколько лет назад. Для себя, без графического интерфейса. Не формула Бурса-Вольфа, а преобразование Гельмерта для сравнительно больших углов вращения последовательно вокруг осей X, Y и Z. Для малых углов, естественно, тоже годится. Могу поделиться.

все куда-то ушли в сторону подбора нового эллипсоида, а вопрос был про афинное преобразование, задаваемое в операторе Coordsys.

Boris, как всегда, прав. Однако изложенный мной алгоритм является решением и этой задачи. Пространство двумерное вместо трёхмерного. Вместо масштаба и матрицы вращения 3×3 — матрица 2×2 из параметров A, B, D, E. Параметры переноса C и F.
X₂ = A X₁ + B Y₁ + C
Y₂ = D X₁ + E Y₁ + F
Бонус общего афинного преобразования: пары параметров A, B и D, E вычисляются независимо. Вместо МНК для четырёх неизвестных A, B, D, E — два МНК для двух неизвестных.
Ну а настоящие геодезисты используют только конформное преобразование: D = −B, E = A.
X₂ = A X₁ + B Y₁ + C
Y₂ =−B X₁ + A Y₁ + F

[Донецков]

Boris, спасибо, что поставили тему на "путь истинный"...
Интересует именно программа вычисления этих параметров, чтобы записать в файл PRJ, а не преобразования одной таблицы в другую (такие есть в поставке MapInfo).
И, естественно, ErnieBoyd прав, что в решении этой конкретной задачи будет применяться конформное преобразование: D = −B, E = A:
X₂ = A X₁ + B Y₁ + C
Y₂ =−B X₁ + A Y₁ + F
во всяком случае в имеющихся определениях МСК другого не встречал...
P.s. все новое это хорошо забытое старое, откапал в своих закромах программу Recalc (которую когда-то бесплатно скачал с сайта GLAB@2007), использовал ее исключительно для быстрого пересчета таблиц из/в план-схема - проекция..., так вот она вычисляет эти 6 параметров, правильность мсам не проверял, т.к. пока сравнить не счем...
Когда-то видел аналогичную программу, только она сразу генерировала текстовую строку с параметрами преобразования для добавления в PRJ...

[Игорь Белов]

Интересует именно программа вычисления этих параметров, чтобы записать в файл PRJ...

Хм, мысль интересная. Правда, программа не должна писать в системный файл MAPINFOW.PRJ. Прилично оставить это пользователю.
Итак, прожект.
Чтобы не зависеть от ошибок пересчёта между системами координат, характерных для MapInfo, лучше написать программу на MapBasic'е или, скажем, C++ с использованием MapInfo API.
Перед пуском программы в MAPINFOW.PRJ должна присутствовать выходная проекция, только без афинного преобразования.
Входные данные — открытый слой в исходной системе координат (не план-схема, а на Земле), содержащий точечные объекты. В таблице присутствуют колонки с координатами в целевой системе координат.
Запускаем программу. В окне диалога выбираем:

• название таблицы;
• колонки координат в конечной СК X и Y;
• колонки для остаточных невязок δX и δY;
• координаты центра преобразования во входной СК X₀₀ и Y₀₀;
• координаты центра преобразования в выходной СК X₀ и Y₀;
• конечную СК (выходная проекция без афинного преобразования);
• метод преобразования — конформное или афинное;
• писать или не писать остаточные невязки.

В текстовый файл программа пишет вычисленные параметры (скажем, в виде строки "Affine Units "m", ..."). В соответствующие колонки таблицы пишутся невязки.
Ну как?

[Игорь Белов]

Простите, inouterer. Если Вы ещё не нашли ответа, вот он.
Итак, после параметров проекции (и перед границами проекции, если они даются) записываются шесть параметров афинного преобразования A, B, C, D, E, F.

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

A = m cos α
B =−m sin α
C = X₀ − (A X₀₀ + B Y₀₀)
D = m sin α
E = m cos α
F = Y₀ − (D X₀₀ + E Y₀₀)

m — масштаб (у Вас единица);
α — разворот (знак соответствует вращению координатных осей по часовой стрелке);
X₀₀, Y₀₀— координаты центра преобразования в исходной системе координат;
X₀, Y₀— координаты центра преобразования в трансформированной системе координат.

А, вот, есть в MapInfo Professional Help:
"Working with Coordinate Systems" → "Working with Coordinate Systems and Projections" → "Understanding Affine Transformations" → "Description of an Affine Transformation"

[Донецков]

Правда, программа не должна писать в системный файл MAPINFOW.PRJ. Прилично оставить это пользователю

- естественно, это должно остаться на "совести" пользователя...

[Игорь Белов]

mifreg — программа, осуществляющая афинное преобразование.
История такова. Очень давно в конторе, куда я пришёл, в Easy Trace под DOS отвекторизовали кучу нетрансформированных сканов в СК на основе растров. Пришлось написать программу для трансформирования векторов в СК топокарты.

Последовательность действий.
Зарегистрировать скан карты в MapInfo.
Экспортировать из Easy Trace векторы в файлы MIF/MID.
Запустить mifreg. Указать пути:

• Raster registration TAB file: файл регистрации растра;
• Raw MIF file: входной MIF в координатах растра;
• Output MIF file: трансформированный MIF.

Радиокнопки определяют начало исходной системы координат на растре:

• MapInfo coordinate — левый верхний пиксел;
• EasyTrace coordinate — левом нижний пиксел.

Жмём [Proceed].
Помимо выходного MIF/MID к файлу регистрации растра создаётся LOG.
LOG содержит такие данные:

• параметры прямого афинного преобразования;
• параметры обратного афинного преобразования;
• альтернативные параметры;
• невязки по точкам привязки.

Альтернативные параметры я вычислял для трансформирования растра в Фотошопе. Нужно изображение повернуть на первый угол, изменить размеры в соответствии с масштабами по главным осям, повернуть на второй угол.

Теперь главное.
Программа может использоваться не только для пар "пиксельные координаты" — "мировые координаты", но для любых пар координат. Возьмём за основу любой файл привязки растра TAB и с помощью хорошего редактора и электронных таблиц создадим новый, в котором для каждой точки первая пара чисел будет координатами в выходной СК, а вторая — во входной. В качестве входного MIF/MID укажем любой в системе CoordSys NonEarth (план-схема). Проверим, что активна радиокнопка [MapInfo coordinate]. Результаты найдём в логе.