преобразование координат и нахождение точки на карте.

[edoroe]

Доброго времени суток.

Что-то я делаю не так для отображения точки с координатами lat0 long0 на карте.

Опишу как делаю. Просьба к многоуважаемой публике указать на мои ошибки.

1) координаты точки получаю в wgs84
2) карта привязана для Ozi. из девяти точек привязки получаю ближайшие 6 к входной точке (lat0 long0)
в результате анализа точек привязки получаю шесть пар x/y координаты в графическом файле и lat/long этих точек
тоесть
x1 y1 lat1 long1
....
x6 y6 lat6 long6

3) преобразовываю все семь (водная точка и 6 точек привязки) lat/long координаты в прямоугольную систему UTM
тоесть получаю
x1 y1 ux1 uy1
....
x6 y6 ux6 uy6

В результате этих преобразований у меня есть 6 точек привязки с координатами в графическом файле (пиксели) и соответствия в прямоугольной системе UTM

и есть входная точка с UTM координатами (я их преобразовал WGS84->UTM)

Тоесть необходимо найти координаты пикселей для входной точки. Тоесть решить полином по шести точкам.
Используя статью
http://gis-lab.info/qa/polynom-calc.html

a) я строю матрицу
1 ux1 uy1 ux1*ux1 ux1*uy1 uy1*uy1
.....
1 ux6 uy6 ux6*ux6 ux6*uy6 uy6*uy6

б) нахожу обратную матрицу для вышеприведенной матрицы
в) перемножаю одностолбцовую матрицу (сорри, не знаю как обозвать ее) составленную из иксовых координат в графическом файле
x1
...
x6
и умножаю ее на обратную, получаю одностолбцовую матрицу шести коэфициентов (a0 - a5) для нахождения Х координаты входной точки
тоже самое проделываю для Y координат
тоесть перемножаю
y1
...
y6, на обратную матрицу полученую в пункте Б), получая коэфициенты (b0-b5) для нахождения Y координаты входной точки

г) вычисляю Х координату в графическом файле входной точки по формуле
Х = $a0 + $a1*$ux0 + $a2*$uy0 + $a3*$ux0*$ux0 + $a4*$ux0*$uy0 + $a5*$uy0*$uy0

в результате получаю бред, координаты не соответствуют действительности.

Где я ошибаюсь ?

Заранее благодарствую за помощь.

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

у вас приложение свое, иначе зачем это вручную считаете?

если свое приложение, проверьтесь сначала по примеру с известными координатами
http://gis-lab.info/qa/polynom-calc-examples.html

[edoroe]

у вас приложение свое, иначе зачем это вручную считаете?

приложение свое. предложеный тест не прошел.
спасибо за помощь. буду разбираться что с матрицами напутал.

[edoroe]

у вас приложение свое, иначе зачем это вручную считаете?

приложение свое. предложеный тест не прошел.
спасибо за помощь. буду разбираться что с матрицами напутал.

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

Когда я на вход тестового решения подал большие координаты, те что в тесте есть как результат, то на выходе получил 499.999 -299.999 - тоесть то что в тесте на входе.

Тоесть нахождение коефициентов работает правильно.
И функция расчета пиксельных координат тоже работает правильно, но для вот такого набора данных

шесть базовых точек:
x y ux uy
-----------------------------
0 1100 654243 5508937
0 2200 654437 5502146
1099 1100 661043 5509135
0 0 654050 5515719
1099 2200 661245 5502344
1099 0 660841 5515917
------------------------------
для входной точки по UTM ux0=656183 uy0=5508415
я получаю -173.30807153438218 1179.2341383858584

причем "Y" координата, примерно та что надо.

Есть какие-либо соображения ?

[edoroe]

Нет я ошибся. Тест прошел ...

Кстати, побаловался с выбором шести точек из девяти, которые откалиброваны в *map файле.

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

Сейчас, при неправильном результате, выбираю шесть ближайших до искомой точек.

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

матрица может выродится если точки в линию выстроились.

[edoroe]

матрица может выродится если точки в линию выстроились.

ну, этого к счасью нет:

0 1100 654243 5508937
0 2200 654437 5502146
1099 1100 661043 5509135
0 0 654050 5515719
1099 2200 661245 5502344
1099 0 660841 5515917

первые две колонки это x и y координаты точек в графическом файле.

[edoroe]

шесть базовых точек:
x y ux uy
-----------------------------
0 1100 654243 5508937
0 2200 654437 5502146
1099 1100 661043 5509135
0 0 654050 5515719
1099 2200 661245 5502344
1099 0 660841 5515917
------------------------------
для входной точки по UTM ux0=656183 uy0=5508415
я получаю -173.30807153438218 1179.2341383858584

причем "Y" координата, примерно та что надо.

Есть какие-либо соображения ?

подставил эти значения в эксельный пример, взятый из
http://gis-lab.info/qa/polynom-calc-examples.html
получил тот-же результат что и выдает моя программа

А должна быть X координата кокло 300 а не -173

В чем причина ?

[gamm]

простая линейная модель на этих данных выдает 310.792 1193.199 (для квадратичной данных не хватает)

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

Спойлер

lm(formula = Xold ~ Xnew + Ynew, data = p)
Residuals: 0.03914 0.66591 -0.01354 -0.70677 -0.67868 0.69394
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.316e+05 3.298e+02 -398.93 3.47e-08 ***
Xnew 1.615e-01 9.516e-05 1696.94 4.51e-10 ***
Ynew 4.706e-03 5.845e-05 80.51 4.22e-06 ***
Multiple R-squared: 1, Adjusted R-squared: 1

lm(formula = Yold ~ Xnew + Ynew, data = p)
Residuals: 0.4874 -0.2235 0.4859 -0.2639 -0.2628 -0.2230
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.902e+05 2.026e+02 4393.58 2.60e-11 ***
Xnew 4.716e-03 5.846e-05 80.67 4.20e-06 ***
Ynew -1.619e-01 3.590e-05 -4510.66 2.40e-11 ***
Multiple R-squared: 1, Adjusted R-squared: 1

[edoroe]

простая линейная модель на этих данных выдает 310.792 1193.199 (для квадратичной данных не хватает)
у вас проблемы со счетом ... кстати, модель у вас в точности линейная.

Огромное спасибо ! Вы прям как Моисей - приоткрыли мне глаза на сущность явлений !

Я понял. Буду пересчитывать по трем точкам.
И понял почему на картах большого масштаба, когда точки привязки находятся в разных UTM зонах, вообще координаты на графическом файле ложны. Это из-за того, что я получаю UTM координаты в пределах шестиградусной зоны. Отсюда следующий вопрос.

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

Что почитать на предмет нахождения формул для пересчета точек, находящихся в разных UTM зонах.
Или как в таком случае поступают нормальные люди ?

[gamm]

1) Я понял. Буду пересчитывать по трем точкам.

2) Что почитать на предмет нахождения формул для пересчета точек, находящихся в разных UTM зонах.
Или как в таком случае поступают нормальные люди ?

1) лучше по 5, конвертом, и контроль еще по 4

2) тут на соседней ветке выложены формулы для UTM , переводите "соседнюю зону" в градусы, а потом в единую зону. ссылка

[edoroe]

1) лучше по 5, конвертом, и контроль еще по 4

Не понял, как по пяти.
Использую формулы по 3 http://gis-lab.info/qa/affine-math.html
и по шести: http://gis-lab.info/qa/polynom-calc.html

а как по пяти ?

[gamm]

1) лучше по 5, конвертом, и контроль еще по 4

Не понял, как по пяти.
Использую формулы по 3 http://gis-lab.info/qa/affine-math.html
и по шести: http://gis-lab.info/qa/polynom-calc.html

а как по пяти ?

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

X*alpha=Y
X'*X*alpha=X'*Y
alpha_bar=(X'*X)^-1*X*Y, X' - транспонирование.

А еще лучше учите R, ниже програмулька, которая решает вашу задачку

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

p<-read.table("foo.dat",header=TRUE)
names(p)
str(p)
#[1] "Xold" "Yold" "Xnew" "Ynew"

lm.x<-lm(Xold~Xnew+Ynew,data=p)
summary(lm.x)

lm.y<-lm(Yold~Xnew+Ynew,data=p)
summary(lm.y)

new.data<-data.frame(Xnew=656183,Ynew=5508415)

x.prd<-predict(lm.x,newdata=new.data)
y.prd<-predict(lm.y,newdata=new.data)
print(c(x.prd,y.prd))

lm.x2<-lm(Xold~Xnew*Ynew+I(Xnew^2)+I(Ynew^2),data=p)
summary(lm.x2)

lm.y2<-lm(Yold~Xnew*Ynew+I(Xnew^2)+I(Ynew^2),data=p)
summary(lm.y2)

new.data<-data.frame(Xnew=656183,Ynew=5508415)

x.prd2<-predict(lm.x2,newdata=new.data)
y.prd2<-predict(lm.y2,newdata=new.data)
print(c(x.prd2,y.prd2))

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

Xold Yold Xnew Ynew
0 1100 654243 5508937
0 2200 654437 5502146
1099 1100 661043 5509135
0 0 654050 5515719
1099 2200 661245 5502344
1099 0 660841 5515917

[edoroe]

преобразование координат и нахождение точки на карте.

Все. Сделал.

Отдельная благодарность
Максиму Дубинину и gamm

Поступил так:
формулу для нахождения искомой точки по шести известным http://gis-lab.info/qa/polynom-calc.html
преобразовал до вида для девяти точек. И использую все девять полученных коэффициентов.

Проверил на картах гугли полученых посредством sasplanet на уровнях z15-z8
Все отлично, кроме уровня z8 - но это вся Украина в размере 1800x1800 пикселей.
Загрузил для успокоения эту карту в OziExplorer - он тоже на этой карте дает ошибку пикселей 7-9.

Думая, из-за крупности карты с этим можно смириться.