Преобразование координат с помощью полиномиальных преобразов

[ALEUT]

Вот возникла ещё одна проблемка. Пытаюсь преобразовать координаты с помощью полиномиальных преобразований с WGS84 в локальные, но выходит какая-то лажа. Есть несколько точек привязки (информация взята из map файла карты Ozi Explorer):
Point01,xy,18933,03532,in, deg, 53, 56.4338,N, 27, 40.3640,E, grid,,,,N
Point02,xy,16166,01398,in, deg, 53, 57.5766,N, 27, 37.6107,E, grid,,,,N
Point03,xy,13875,05416,in, deg, 53, 55.1843,N, 27, 35.5749,E, grid,,,,N
Point04,xy,05823,12840,ex, deg, 53, 50.5447,N, 27, 28.1406,E, grid,,,,N

MMPXY,1,0,0
MMPXY,2,26733,0
MMPXY,3,26733,18360
MMPXY,4,0,18360
MMPLL,1,27.366520,53.964417
MMPLL,2,27.795383,53.978703
MMPLL,3,27.808283,53.802109
MMPLL,4,27.379421,53.787823

0;0;27.366520;53.964417
26733;0;27.795383;53.978703
26733;18360;27.808283;53.802109
0;18360;27.379421;53.787823

Я их и использовал. Тестовую точку взял как 53.90291,27.48884, в локальных должно получится 7290, 6770 с незначительными отклонениями, а получаются вообще отрицательные координаты, которых не может быть.

[Максим Дубинин]

Вот что получается у меня, беру 3 из 4 точек, делаю аффинное преобразование, привожу код на R:

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

x = c(0,26733,26733)
y = c(0,0,18360)
x2 = c(27.36652,27.794383,27.808283)
y2 = c(53.964417,53.978703,53.802109)
mat <- matrix( c(1, 1, 1,x2, y2), nrow = 3, ncol = 3)
aninv = solve(mat, x)
bninv = solve(mat, y)

получаем коэффициенты

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

> aninv
[1] -1970083.127    62316.505     4905.033
> bninv
[1] 5501077.137    3462.285 -103694.770

тестируем

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

testpoint = c(27.48884,53.90291)
xpred = aninv[1] + aninv[2]* testpoint[1] + aninv[3]* testpoint[2]
ypred = bninv[1] + bninv[2]* testpoint[1] + bninv[3]* testpoint[2]

получается: xpred = 7320.861, ypred = 6801.461
ошибка:

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

sqrt((7320.861-7290)^2+(6801.461-6770)^2)
[1] 44.07035

[ALEUT]

А если полиномом? Вот 6 точек:
5310;11210;27.45996;53.85775
10658;12435;27.54727;53.84891
17980;11170;27.66289;53.86579
16106;4890;27.62995;53.92577
15094;2335;27.61077;53.95008
11324;4883;27.55240;53.92338
4988;6280;27.45102;53.90684

x = c(5310, 10658,17980,16106,15094,11324,4988)
y = c(11210,12435,11170,4890,2335,4883,6280)
x2 = c(27.45996,27.54727,27.66289,27.62995,27.61077,27.55240,27.45102)
y2 = c(53.85775,53.84891,53.86579,53.92577,53.95008,53.92338,53.90684)

[Максим Дубинин]

считаем:

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

x = c(5310, 10658,17980,16106,15094,11324)
y = c(11210,12435,11170,4890,2335,4883)
x2 = c(27.45996,27.54727,27.66289,27.62995,27.61077,27.55240)
y2 = c(53.85775,53.84891,53.86579,53.92577,53.95008,53.92338)
mat <- matrix( c(1, 1, 1, 1, 1, 1, x2, y2, x2^2, x2*y2, y2^2), nrow = 6, ncol = 6)
aninv = solve(mat, x)
bninv = solve(mat, y)

коэффициенты

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

aninv
[1] 79372078.038   777007.838 -3378967.426     8021.566   -21478.112 36886.886
bninv
[1] -15391607.893    505362.039    413051.406      1018.563    -10353.857 -2129.111

тестируем

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

testpoint = c(27.48884,53.90291)
xpred = aninv[1] + aninv[2]* testpoint[1] + aninv[3]* testpoint[2] + aninv[4]* testpoint[1]^2 + aninv[5]* testpoint[1]*testpoint[2] + aninv[6]*testpoint[2]^2
ypred = bninv[1] + bninv[2]* testpoint[1] + bninv[3]* testpoint[2] + bninv[4]* testpoint[1]^2 + bninv[5]* testpoint[1]*testpoint[2] + bninv[6]*testpoint[2]^2

Результат:

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

xpred
[1] 7324.237
ypred
[1] 6755.628

Ошибка:

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

sqrt((7324.237-7290)^2+(6755.628-6770)^2)
[1]37.13121

[ALEUT]

В очередной раз меня выручили)