Преобразование координат из МГГТ в WGS-84

Добрый день, уважаемые форумчане.
Не смог сориентироваться в какой категории создавать тему, надеюсь плюс-минус угадал.
Столкнулся с проблемой преобразования координат из плоских топоцентрических МГГТ в геоцентрическую WGS-84. И это для меня превратилось в какую то принципиальную битву.
Согласно Положению о пространственной МСК Москвы:
1. МГГТ - Проецируется на эллипсоид Бесселя по проекции Гаусса-Крюгера.
2. Преобразовывать координаты надо цепочкой: МГГТ-ПМСК-СК95-ПЗ90-WGS-84.
3. Ключи для перехода МГГТ-ПМСК есть, а вот от ПМСК к СК95 под грифом: не хотим делиться. А ключи остальной цепочки-СК-95-ПЗ90-WGS-84 представлены в госте.
Итак, сайт геобридж, достаточно неплохо трансформирует МСК МГГТ в СК WGS84. Я достал ключи перехода оттуда, но вот только мои результаты с сайтом не бьются. В координатах 0,0 расхождение в метр, а вот отдаляясь от точки пересечения, расхождение увеличивается.
Я организовывал перевод через пространственные координаты. т.е.:
x1,y1 -> B1,L1 -> X1,Y1,Z1 -> X2, Y2, Z2 -> B2,L2
Что могло пойти не так?
- Трансформация XYZ1 в XYZ2 по параметрам гельмерта (Хельмерта) достаточно простая, тут расхождений быть не должно, но подумываю написать функцию для перевода BL1 в BL2. Не знаю насколько увеличиться точность, никогда не сравнивал такие результаты, но предполагаю что они будут идентичны, хотя перевод из gпространственных в геодезические для меня не однозначен. Переход от пространственных в геодезические я выполняю и по гостовской формуле и по формуле Огородовой.
- А вот работа с проекцией гаусса меня больше всего напрягает. Во 1.) В Положение о ПМСК написано: Переводим согласно Госту, но в госте указана упрощенная формула с коэффициентами для эллипсоида Красовского, а у нас Бесселя. Пришлось достать свои старые работы с МИИГАиК. Формулы на всякий случай приложу картинками.
- На сколько я понял, в ПО ГИС, mapinfo arcgis и т.д. Нет проекции Гаусса-Крюгера, а эти преобразования они выполняют по проекции Меркатора. Что тоже меня немного озадачивает.

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

 

Таблицы MS Excel для преобразования координат: [ellipsoid]: gauss-kruger.xls

 

Большое спасибо. Скачал на вашем github. В данном документе представлена реализация преобразования из гаусса крюгера для СК, а не МСК. Приращение долготы близкое к нулю я смог получить при значении ординаты:2375591,790874.

И на первый взгляд кажется что у вас представлены чуть упрощённые формулы. В принципе используемые мной формулы представлены в первом посте в виде вложений изображений.
Возможно, я как то не правильно вас понял, но к сожалению данная реализация мне не помогает :(
Большое Вам спасибо, за то, что попытались помочь.

Оффтоп

Мой расчет:

Код:

public class MSK
{
public double B;//Широта точки начала отсчета координат
public double L;//Долгота осевого меридиана
public double m;//Масштабный коэффициент на осевом меридиане
public double dx;//Сдвиг системы координат по оси абсцисс X (координаты начала местной системы в местной системе) 
public double dy;//Сдвиг системы координат по оси ординат Y, (координаты начала местной системы в местной системе) 
public double r;//Угол поворота осей координат местной системы в точке начала местной системы координат
public Ellipsoid el;//Эллипсоид
public MSK(double p_b, double p_l, double p_m, double p_dx, double p_dy, double p_r, Ellipsoid _el) {
B = p_b; 
L = p_l; 
m = p_m; 
dx = p_dx; 
dy = p_dy; 
r = p_r;
el = _el;
}
}

Код:

public class Ellipsoid
{
public double a;//Большая полуось эллипсоида
public double b;//Малая полуось эллипсоида
public double alfa;//Альфа, коэфициент сжатия 1/f
public double c;//Фокальное расстояние
public double e;//Экцентрисет
public double e2;//квадрат Экцентрисет
public double _e;//второй Экцентрисет
public double _e2;//квадрта второго Экцентрисет
public Ellipsoid(double semiAxesA, double semiAxesB) {
a = semiAxesA;
b = semiAxesB;
alfa = (a - b) / a;
c = Math.Sqrt(Math.Pow(a, 2) - Math.Pow(b, 2));
e = c / a;
e2 = Math.Pow(e, 2);
_e = c / b;
_e2 = Math.Pow(_e, 2);
}
public Ellipsoid(double semiAxesA, double _alfa, bool ch) {
a =semiAxesA;
alfa= _alfa;
b = a * (1.0 - alfa);
c = Math.Sqrt(Math.Pow(a, 2) - Math.Pow(b, 2));
e = c/a;
e2 = Math.Pow(e, 2);
_e = c / b;
_e2 = Math.Pow(_e, 2);
}
public double GetNd(double Bdeg) {
double res = a/Math.Sqrt(1-e2*Math.Pow(Math.Sin(SI.Conv_DecToRad(Bdeg)),2));
return res;
}
public double GetNr(double Brad) {
double res = a/Math.Sqrt(1-e2*Math.Pow(Math.Sin(Brad),2));
return res;
}
}

Код:

public class PointXY
{
public double X,Y;
public PointXY(double x, double y) {
X = x;
Y = y;
}
}
public class PointBL
{
public double Bd, Br, Bcos, Bsin;
public double Ld, Lr, Lcos, Lsin;
public PointBL(double b, double l) {
Bd = b; 
Br = SI.Conv_DecToRad(b);
Bcos = Math.Cos(Br); Bsin = Math.Sin(Br);
Ld = l;
Lr = SI.Conv_DecToRad(l);
Lcos = Math.Cos(Lr); Lsin = Math.Sin(Lr);
}
}

Код:

static public class Gaus
{
static public PointBL convertXY(PointXY pt, MSK mSK) {
double x = pt.X - mSK.dx;
double y = pt.Y - mSK.dy;
double e0 = 1.0 - (1.0 / 4.0) * mSK.el.e2;
e0 -= (3.0 / 64.0) * Math.Pow(mSK.el.e, 4.0);
e0 -= (5.0 / 256.0) * Math.Pow(mSK.el.e, 6.0);
e0 -= (175.0 / 16384.0) * Math.Pow(mSK.el.e, 8.0);
e0 -= (441.0 / 65536.0) * Math.Pow(mSK.el.e, 10.0);
//e0 -= (4851 / 1048576.0) * Math.Pow(FirstE, 12);
double c2 = 3.0/8.0 * mSK.el.e2;
c2 += 3.0 / 16.0 * Math.Pow(mSK.el.e, 4.0);
c2 += 213.0 / 2048.0 * Math.Pow(mSK.el.e, 6.0);
c2 += 255.0 / 4096.0 * Math.Pow(mSK.el.e, 8.0);
c2 += 166479.0 / 655360.0 * Math.Pow(mSK.el.e, 10.0);
//c2 +=  * Math.Pow(FirstE, 12);
double c4 = 21.0/256.0*Math.Pow(mSK.el.e, 4.0);
c4 += 21.0 / 256.0 * Math.Pow(mSK.el.e, 6.0);
c4 += 533.0 / 8192.0 * Math.Pow(mSK.el.e, 8.0);
c4 -= 120563.0 / 327680.0 * Math.Pow(mSK.el.e, 10.0);
//c4
double c6 = 151.0/6144.0*Math.Pow(mSK.el.e, 6.0);
c6 += 147.0 / 4096.0 * Math.Pow(mSK.el.e, 8.0);
c6 += 2732071.0 / 9175040.0 * Math.Pow(mSK.el.e, 10.0);
//c6
double c8 = 1097.0/131072.0 * Math.Pow(mSK.el.e, 8.0);
c8 -= 273697.0 / 4587520.0 * Math.Pow(mSK.el.e, 10.0);
//c8
double Bx = x/(mSK.el.a*e0);
Bx += c2 * Math.Sin((2.0 * x) / (mSK.el.a * e0));
Bx += c4 * Math.Sin((4.0 * x) / (mSK.el.a * e0));
Bx += c6 * Math.Sin((6.0 * x) / (mSK.el.a * e0));
Bx += c8 * Math.Sin((8.0 * x) / (mSK.el.a * e0));
double Nx = mSK.el.a / Math.Sqrt(1.0-mSK.el.e2*Math.Pow(Math.Sin(Bx),2.0));
double nu = Math.Sqrt(mSK.el._e2)*Math.Cos(Bx);
double A1 = 1.0 / (Nx * Math.Cos(Bx));
double A2 = 1.0/2.0 * (Math.Tan(Bx)/Math.Pow(Nx,2.0))*(-1.0-Math.Pow(nu,2.0));
double A3 = -1.0/6.0*(A1/Math.Pow(Nx,2.0))*(1.0+2.0*Math.Pow(Math.Tan(Bx),2.0)+Math.Pow(nu,2.0));
double A4 = -1.0/ 12.0 * (A2 / Math.Pow(Nx, 2.0)) * (5.0 + 3.0 * Math.Pow(Math.Tan(Bx), 2.0) + Math.Pow(nu, 2.0) - 9.0 * Math.Pow(nu, 2.0) * Math.Pow(Math.Tan(Bx), 2.0) - 4.0 * Math.Pow(nu, 4.0));
double A5 = 1.0 / 120.0 * (A1 / Math.Pow(Nx, 4.0)) * (5.0 + 28.0 * Math.Pow(Math.Tan(Bx), 2.0) + 24.0 * Math.Pow(Math.Tan(Bx), 4.0) + 6.0 * Math.Pow(nu, 2.0) + 8.0 * Math.Pow(nu, 2.0) * Math.Pow(Math.Tan(Bx), 2.0));
double A6 = 1.0 / 360.0 * (A2 / Math.Pow(Nx, 4.0)) * (61.0 + 90.0 * Math.Pow(Math.Tan(Bx), 2.0) + 45.0 * Math.Pow(Math.Tan(Bx), 4.0) + 46.0 * Math.Pow(nu, 2.0) - 252.0 * Math.Pow(nu, 2.0) * Math.Pow(Math.Tan(Bx), 2.0) - 90.0 * Math.Pow(nu, 2.0) * Math.Pow(Math.Tan(Bx), 4.0));
double A7 = -1.0 / 5040.0 * (A1 / Math.Pow(Nx, 6.0)) * (61.0 + 662.0 * Math.Pow(Math.Tan(Bx), 2.0) + 1320.0 * Math.Pow(Math.Tan(Bx), 4.0) + 720.0 * Math.Pow(Math.Tan(Bx), 6.0));
double A8 = -1.0 / 20160.0 * (A2 / Math.Pow(Nx, 6.0)) * (1385.0 + 3633.0 * Math.Pow(Math.Tan(Bx), 2.0) + 4095.0 * Math.Pow(Math.Tan(Bx), 4.0) + 1575.0 * Math.Pow(Math.Tan(Bx), 6.0));
double B = Bx + A2*Math.Pow(y,2.0)+A4*Math.Pow(y,4.0)+A6*Math.Pow(y,6.0)+A8*Math.Pow(y,8.0);
double L = A1*y+A3*Math.Pow(y,3.0)+A5*Math.Pow(y,5.0)+A7*Math.Pow(y,7.0);
double resB = mSK.B + (B*(180.0/Math.PI));
double resL = mSK.L + (L*(180.0/Math.PI));
PointBL res = new PointBL(resB,resL);
return res;
}
static public PointBL convertXY_2(PointXY pt, MSK mSK) {
double n1_=mSK.el.e2/2;
double n2_=Math.Pow(mSK.el.e,4.0)/8.0;
double n3_=Math.Pow(mSK.el.e,6.0)/16.0;
double n4_=1.0/(6.0*(1.0-mSK.el.e2));
double n5_=(2.0-9.0*mSK.el.e2)/(12.0*(1.0-mSK.el.e2));
double n6_=(4.0*mSK.el.e2-39.0*Math.Pow(mSK.el.e,4.0))/(48.0*(1.0-mSK.el.e2));
double n7_=(3.0*Math.Pow(mSK.el.e,4.0))/(16.0*(1.0-mSK.el.e2));
double n8_=(5.0+6.0*mSK.el.e2+3.0*Math.Pow(mSK.el.e,4.0))/120.0;
double n9_=(192.0-240.0*mSK.el.e2-123.0*Math.Pow(mSK.el.e,4.0))/1280.0;
double n10_=(32.0-1376.0*mSK.el.e2+609.0*Math.Pow(mSK.el.e,4.0))/3840.0;
double n11_=(mSK.el.e2-69.0*Math.Pow(mSK.el.e,4.0))/240.0;
double n12_=(61.0+46.0*mSK.el.e2)/(5040.0*(1.0-mSK.el.e2));
double n13_=(958.0-1361.0*mSK.el.e2)/(10080.0*(1.0-mSK.el.e2));
double n14_=(358.0-4395.0*mSK.el.e2)/(10080.0*(1.0-mSK.el.e2));
double n15_=(815.0*mSK.el.e2-2.0)/(10080.0*(1.0-mSK.el.e2));
double n16_=0.0038;
double n17_=0.0524;
double n18_=(18270.0-113789.0*mSK.el.e2)/362880.0;
double n19_=(1636.0-72123.0*mSK.el.e2)/362880.0;
//
double k1_=1.0/(4.0*(1.0-mSK.el.e2));
double k2_=mSK.el.e2/(2.0*(1.0-mSK.el.e2));
double k3_=Math.Pow(mSK.el.e,4.0)/(4.0*(1.0-mSK.el.e2));
double k4_=(5.0+6.0*mSK.el.e2+3.0*Math.Pow(mSK.el.e,4.0))/48.0;
double k5_=(1.0+14.0*mSK.el.e2+15.0*Math.Pow(mSK.el.e,4.0))/24.0;
double k6_=(8.0*mSK.el.e2+31.0*Math.Pow(mSK.el.e,4.0))/24.0;
double k7_=(2.0*Math.Pow(mSK.el.e,4.0))/3.0;
double k8_=(61.0+107.0*mSK.el.e2)/1440.0;
double k9_=(16.0+333.0*mSK.el.e2)/720.0;
double k10_=(2.0+87.0*mSK.el.e2)/180.0;
double k11_=(17.0*mSK.el.e2)/90.0;
double k12_=(277.0-1108.0*mSK.el.e2)/16128.0;
double k13_=(29.0-116.0*mSK.el.e2)/4480.0;
double k14_=(41.0-164.0*mSK.el.e2)/3360.0;
double k15_=(17.0-68.0*mSK.el.e2)/5040.0;
//
double C0 = 1.0+(3.0*mSK.el.e2)/4.0+(45.0*Math.Pow(mSK.el.e,4.0))/64.0+(176.0*Math.Pow(mSK.el.e,6.0))/256.0+(11025.0*Math.Pow(mSK.el.e,8.0))/16384.0;
double C2 = (3.0*mSK.el.e2)/4.0+(15.0*Math.Pow(mSK.el.e,4.0))/16.0+(525.0*Math.Pow(mSK.el.e,6.0))/512.0+(2205.0*Math.Pow(mSK.el.e,8))/2048.0;
double C4 = (15.0*Math.Pow(mSK.el.e,4.0))/64.0+(105.0*Math.Pow(mSK.el.e,6.0))/256.0+(2205.0*Math.Pow(mSK.el.e,8.0))/4096.0;
double C6 = (35.0*Math.Pow(mSK.el.e,6.0))/512.0+(315.0*Math.Pow(mSK.el.e,8.0))/2048.0;
//
double a0=mSK.el.a*(1.0-mSK.el.e2)*C0;
double beta=pt.X/a0;
double q2=C2/(2.0*a0)+C2*C4/(8.0*Math.Pow(a0,2.0))-Math.Pow(C2,3.0)/(16.0*Math.Pow(a0,3.0));
double q4 = -C4/(4.0*a0)+Math.Pow(C2,2.0)/(4.0*Math.Pow(a0,2.0));
double q6 = C6/(6.0*a0)-(3.0*C2*C4)/(8.0*Math.Pow(a0,2.0))+3.0*Math.Pow(C2,3.0)/(16.0*Math.Pow(a0,3.0));
double p1_=q2+2.0*q4+3.0*q6;
double p2_=4.0*q4+16.0*q6;
double p3_=16.0*q6;
double B0 = beta+Math.Sin(2.0*beta)*(p1_-p2_*Math.Pow(Math.Sin(beta),2.0)+p3_*Math.Pow(Math.Sin(beta),4.0));
//
double _b=B0;
double A2 = -1.0/(Math.Pow(mSK.el.a,2.0)*Math.Pow(Math.Cos(_b),2.0))*Math.Sin(2.0*B0)*(k1_-k2_*Math.Pow(Math.Sin(B0),2.0)+k3_*Math.Pow(Math.Sin(B0),4.0));
double A4 = 1.0/(Math.Pow(mSK.el.a,4.0)*Math.Pow(Math.Cos(_b),4.0))*Math.Sin(2.0*B0)*(k4_-k5_*Math.Pow(Math.Sin(B0),2.0)+k6_*Math.Pow(Math.Sin(B0),4.0)-k7_*Math.Pow(Math.Sin(B0),6.0));
double A6 = -1.0/(Math.Pow(mSK.el.a,6.0)*Math.Pow(Math.Cos(_b),6.0))*Math.Sin(2.0*B0)*(k8_-k9_*Math.Pow(Math.Sin(B0),2.0)+k10_*Math.Pow(Math.Sin(B0),4.0)-k11_*Math.Pow(Math.Sin(B0),6.0));
double A8 = 1.0/(Math.Pow(mSK.el.a,8.0)*Math.Pow(Math.Cos(_b),8.0))*Math.Sin(2.0*B0)*(k12_-k13_*Math.Pow(Math.Sin(B0),2.0)+k14_*Math.Pow(Math.Sin(B0),4.0)-k15_*Math.Pow(Math.Sin(B0),6.0));
double B1 = 1.0/(mSK.el.a*Math.Cos(_b))*(1.0-n1_*Math.Pow(Math.Sin(B0),2.0)-n2_*Math.Pow(Math.Sin(B0),4.0)-n3_*Math.Pow(Math.Sin(B0),6.0));
double B3 = -1.0/(Math.Pow(mSK.el.a,3.0)*Math.Pow(Math.Cos(_b),3.0))*(n4_+n5_*Math.Pow(Math.Sin(B0),2.0)-n6_*Math.Pow(Math.Sin(B0),4.0)-n7_*Math.Pow(Math.Sin(B0),6.0));
double B5 = 1.0/(Math.Pow(mSK.el.a,5.0)*Math.Pow(Math.Cos(_b),5.0))*(n8_+n9_*Math.Pow(Math.Sin(B0),2.0)+n10_*Math.Pow(Math.Sin(B0),4.0)-n11_*Math.Pow(Math.Sin(B0),6.0));
double B7 = -1.0/(Math.Pow(mSK.el.a,7.0)*Math.Pow(Math.Cos(_b),7.0))*(n12_+n13_*Math.Pow(Math.Sin(B0),2.0)+n14_*Math.Pow(Math.Sin(B0),4.0)-n15_*Math.Pow(Math.Sin(B0),6.0));
double B9 = -1.0/(Math.Pow(mSK.el.a,9.0)*Math.Pow(Math.Cos(_b),9.0))*(n16_+n17_*Math.Pow(Math.Sin(B0),2.0)+n18_*Math.Pow(Math.Sin(B0),4.0)+n19_*Math.Pow(Math.Sin(B0),6.0));
//
double B = B0+A2*Math.Pow(pt.Y,2.0)+A4*Math.Pow(pt.Y,4.0)+A6*Math.Pow(pt.Y,6.0)+A8*Math.Pow(pt.Y,8.0);
double L = B1*pt.Y+B3*Math.Pow(pt.Y,3.0)+B5*Math.Pow(pt.Y, 5.0)+B7*Math.Pow(pt.Y,7.0)+B9*Math.Pow(pt.Y,9.0);
double resB = mSK.B+(B*(180/Math.PI));
double resL = mSK.L+(L*(180/Math.PI));
return new PointBL(resB, resL);
}
}

--- Сообщения объединены, 27 фев 2024, Оригинальное время сообщения: 27 фев 2024 ---

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

--- Сообщения объединены, 27 фев 2024 ---

Я сверил. Результат разложения у нас с вами одинаковый. Но вот только geobridge как я понимаю считает как то по другому.

 

Выходит, что расчёты верны. "Неверным" остаётся человекозависимое: сдвиги "не те", ключи "не те". Всё это скрыто от любопытных глаз. Надо "колупать замазку". ;)

 

Да, на сайте в консоли упоминается проекция Меркатора, а не Гаусса-Крюгера. У меня почему то такое ощущение, что они развёртывают плоские координаты по другому, но очень хорошо подобрали(вычислили) ключи под себя.

 

Фёдоров Роман, много раз поднимали тему за ПМСК. И параметры здесь приводили. Это зафиксировано на конкретную эпоху. Ну ты знаешь, и знаешь какой СК.
Ты мне объясни, зачем тебе понадобилась геоцентрическая??

--- Сообщения объединены, 27 фев 2024, Оригинальное время сообщения: 27 фев 2024 ---

5

 

Ну изначально для выгрузки точек с мггт в kml, для геологов. Ну они мне обычно скидывают файл кадовский с блоками скважин enggeo, я выгружал точки, кидал в геобридж, и так этот вопрос решался.
Ну решил я внезапно автоматизировать этот процесс. Написал основную структуру, функции перевода, подглядел ключи с геобридж (не ПМСК, а МГГТ-WGS84). Тестирую, а оно не бьется. А у них бьётся.
И вот и любопытно дико, и охота всё доделать.
Ну как так, одни и те же ключи, а результат разный с ними. А расчет у меня, такой как учили в университете. (Хотя допускаю что может где то я накосячил, но пока не нашёл где).
Какой то специфической важной задачи конечно же нет.

 

Чем дальше я углубляюсь, тем больше чувствую себя дураком.
Я уже предполагаю, что я просто что то не то делаю в целом. Я сделал тест:
Преобразовал точки (0,0),(0,1000),(1000,0) в геобридж и у себя.
Точки 1,2,3 - Геобридж.
Точки 4,5,6 - Мои.

1. Странное отклонение точки 1 и 4 в 21 метр под углом в 45 градусов, говорит, что в этом виновата не проекция (так как точки 0,0 строго определены в геодезических координатах B L).
2. Оси повернуты верно.
3. Что то не так с определением долготы
Очень интересно познать источник своей глупости

--- Сообщения объединены, 28 фев 2024, Оригинальное время сообщения: 28 фев 2024 ---

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

Но смещение в 23 метра все равно остается. Хотя ключи одинаковые.

--- Сообщения объединены, 28 фев 2024 ---

Если взять данные с геобридж

Код:

{
    "srsCodeInput": "EPSG:6335000",
    "queue": [],
    "srsCode": "EPSG:6335000",
    "srsAuth": "epsg",
    "srsProjNumber": "6335000",
    "defData": "+title= Московская СК (МГГТ)    +proj=tmerc +lat_0=55.66666666667 +lon_0=37.5 +k=1 +x_0=16.098 +y_0=14.512 +ellps=bessel +towgs84=316.151,78.924,589.650,-1.57273,2.69209,2.34693,8.4507 +units=m +no_defs",
    "title": " Московская СК (МГГТ)    ",
    "projName": "tmerc",
    "lat0": 0.9715666169435683,
    "long0": 0.6544984694978736,
    "k0": 1,
    "x0": 16.098,
    "y0": 14.512,
    "ellps": "bessel",
    "datum_params": [
        316.151,
        78.924,
        589.65,
        -0.0000076248102069140055,
        0.000013051620627781706,
        0.000011378237726064032,
        1.0000084507
    ],
    "units": "m",
    "a": 6377397.155,
    "rf": 299.1528128,
    "ellipseName": "Bessel 1841",
    "b": 6356078.962818189,
    "a2": 40671194472602.09,
    "b2": 40399739781579.94,
    "es": 0.006674372231802045,
    "e": 0.0816968312225269,
    "ep2": 0.006719218799174659,
    "axis": "enu"
}

Можно наблюдать строчки:

Код:

"x0": 16.098, 
"y0": 14.512,

Я не понимаю в связи с чем это добавлено, но если я вношу в прямоугольные исходные x y эти поправки, то мои вычисления определяются как:
1.108,1.552
1.220,1003.835
1002.958,1.552
т.е. по широте и долготе за каждый километр невязка увеличивается на пару метров.

 

Не заметил как то ваше сообщение ранее... Всё подзабыто конечно. Геодезию люблю, но время всегда было забито тривиальными задачами, не связанные не с математикой, не с геодезией. Прикладная геодезия она такая, заставляет деградировать что ли -)

Спасибо большое, тут много полезной информации!

 

Вообще на сайте сети мггт все разжовано; как уйти от секретности производства полевых работ, и как и как вернуть материалы заказчику в нужной СК.

 

Роман, а у Вас есть статистические данные (RINEX) с баз МГГТ?

 

Нет.

 

Плохо.

Эллипсоид - WGS-84 не просто так в ПМСК г. Москвы использован.

ПМСК г. Москвы и WGS-84 жёстко связаны между собой параметрами ИГД. И найти их проблем больших нет.

 

С ключами ПМСК всё понятно в принципе, официальный нет, но они мне не очень и нужны оказались. Мне достаточно было взять параметры на geobridge от мггт к wgs84.
Я застрял на проекции гаусса. Хотя вовсе не ожидал такого исхода.
По ссылкам на темы, которые добрые форумчане приложили выше, я наткнулся на очень полезное сообщение от stout: https://geodesist.ru/threads/preobr...-ploskoj-sistemy-v-druguju.42027/#post-479425
И коэффициенты посчитал и приведенные методы реализовал. кроме вот алгоритма Кленшоу который упоминается в сообщении. и вижу расхождение по x в случае если просто перехожу в геодезические, и перевожу повторно в проекцию.

Код:

x:0;y:0
B:54.9999704712036;L:37.5
x:-3.22434993926436;y:0

Если я до этого писал перевод плюс минус как учили в универе, то сейчас я взял готовый алгоритм, который используется в ПО форумчан, но у меня что то не так. Я взял центр проекции. и ошибка получается явно есть в переводе из геодезических в прямоугольные. Я сначала согрешил на малое количество доступных мне цифр в мантиссе и перешел к decimal на C#. Но это не помогло :-(. В свободное время постоянно возвращаюсь к этой задачке, но понять свою ошибку так и не получается.

 

См. Таблицы MS Excel для преобразования координат: [ellipsoid]: gauss-kruger.xls
Там же и сверяйся.

PS: Понимаю, что повтор. Но что же ещё остаётся.

 

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

Но почему у меня не получилось по источникам от stout... Я не понимаю. Я опечатки поправил конечно свои. и на 1 км перевода туда и обратно невязка меньше см у меня сейчас. И несколько сантиметров на 10км. Коэффициенты разложения у нас с вами разные. Я думаю вы ознакомлены с A HIGHLY ACCURATE WORLD WIDE ALGORITHM FOR THE TRANSVERSE MERCATOR MAPPING (ALMOST) и The Universal Grids and the Transverse Mercator and Polar Stereographic Map Projections, на которые ссылался stout. У меня и с английским все печально, но проблема же явно не у них, а что то сделал или не так понял я. Приложу вырезки формул по которым переводил я в геодезические.

 

Да, @stout не раз упоминал, что я пользую достаточно грубое решение, и что есть более элегантные и сходящиеся разложения. Но мне хватает моих таблиц. C Proj у них расхождений не наблюдалось.

 

В скриншоте из ворда, у меня опечатка с бетами. Написал их из другой колонки. В целом у меня все сходиться с геобридж по вашему @zvezdochiot коду. Большое вам спасибо. Надеюсь потом все таки попробую осилить более точные методы для себя. Вдруг кому пригодиться. Ну и просто как результат темы, код на C#

Код:

namespace Coord
{
public class Projection
{
public string Name { get; }//Имя проекции
public string Description { get; }//Описание проекции
public string Type { get; }//Тип проекции
public decimal Lat { get; }//Начальная широта в радианах
public decimal Lon {  get; }//Начальная долгота в радианах
public decimal Scale { get; }//Масштабный коэффициент на осевом меридиане
public decimal Dn { get; }//Сдвиг координаты север
public decimal De {  get; }//Сдвиг координаты восток
public decimal r { get; }//Угол поворота осей координат в местной системе координат в радианах
public Ellipsoid Ellips { get; }//Эллипсоид
public Projection (string name, string description,Ellipsoid ellipsoid, string type, decimal latDec, decimal lonDec, decimal scale, decimal rotation, decimal dx, decimal dy) {
Name = name;
Description = description;
Ellips = ellipsoid;
Type = type;
Lon = DecMath.Rad(lonDec);
Lat = DecMath.Rad(latDec);
Scale = scale;
r = DecMath.Rad(rotation);
Dn = dx;
De = dy;
}
public PtBL XY2BL(PtXY point) {
PtBL res;
if (Type == "tMerc") {
res = XY2BL_tMerc(point);
} else {
res = null;
}
return res;
}
public PtXY BL2XY(PtBL point) {
PtXY res;
if(Type == "tMerc") {
res = BL2XY_tMerc(point);
} else {
res = null;
}
return res;
}
private PtBL XY2BL_tMerc(PtXY p) {
decimal bt1 = 0.5m * Ellips.n1 - 2.0m/3.0m*Ellips.n2 + 37.0m/96.0m*Ellips.n3-1.0m/360.0m*Ellips.n4;
decimal bt2 = 1.0m/48.0m*Ellips.n2+1.0m/15.0m*Ellips.n3-437.0m/1440.0m*Ellips.n4;
decimal bt3 = 17.0m/480.0m*Ellips.n3-37.0m/840.0m*Ellips.n4;
decimal bt4 = 4397.0m/161280.0m*Ellips.n4;
 
decimal kA = Ellips.e2 + Ellips.e4 + Ellips.e6 + Ellips.e8;
decimal kB = -(7.0m * Ellips.e4 + 17.0m * Ellips.e6 + 30.0m * Ellips.e8)/6.0m;
decimal kC = (224.0m * Ellips.e6 + 889.0m * Ellips.e8)/120.0m;
decimal kD = -(4279.0m * Ellips.e8) / 1260.0m;
 
 
decimal x1 = (Lat2X(Lat) + (p.x - Dn)) / (Scale * Ellips.R);//Lat2X(Lat)
decimal y1 = (p.y - De) / (Scale*Ellips.R);
 
decimal xp = x1 - bt1 * DecMath.Sin(2.0m * x1) * DecMath.Cosh(2.0m * y1)  - bt2 * DecMath.Sin(4.0m * x1) * DecMath.Cosh(4.0m * y1) - bt3 * DecMath.Sin(6.0m * x1) * DecMath.Cosh(6.0m * y1)  - bt4 * DecMath.Sin(8.0m * x1) * DecMath.Cosh(8.0m * y1);
decimal yp = y1 - bt1 * DecMath.Cos(2.0m * x1) * DecMath.Sinh(2.0m * y1) - bt2 * DecMath.Cos(4.0m * x1) * DecMath.Sinh(4.0m * y1) - bt3 * DecMath.Cos(6.0m * x1) * DecMath.Sinh(6.0m * y1) - bt4 * DecMath.Cos(8.0m * x1) * DecMath.Sinh(8.0m * y1);
decimal bp = DecMath.Arcsin(DecMath.Sin(xp)/DecMath.Cosh(yp));
decimal dl = DecMath.Arctan(DecMath.Sinh(yp)/DecMath.Cos(xp));
decimal L = DecMath.Dec(Lon + dl);
decimal B = bp + DecMath.Sin(bp)*DecMath.Cos(bp)*(kA + kB * DecMath.DecimalPow(DecMath.Sin(bp),2) + kC * DecMath.DecimalPow(DecMath.Sin(bp),4) + kD * DecMath.DecimalPow(DecMath.Sin(bp),6));
PtBL res = new PtBL(DecMath.Dec(B), L);
return res;
}
private PtXY BL2XY_tMerc(PtBL p) {
decimal bt1 = 0.5m * Ellips.n1 - 2.0m/3.0m * Ellips.n2 + 5.0m/16.0m * Ellips.n3 + 41.0m / 180.0m * Ellips.n4;
decimal bt2 = 13.0m / 48.0m * Ellips.n2 -3.0m/5.0m*Ellips.n3 + 557.0m/1440.0m*Ellips.n4;
decimal bt3 = 61.0m/240.0m*Ellips.n3 - 103.0m/140.0m*Ellips.n4;
decimal bt4 = 49561.0m/161280.0m*Ellips.n4;
decimal kA = Ellips.e2;
decimal kB = (5.0m * Ellips.e4 - Ellips.e6)/6.0m;
decimal kC = (104.0m * Ellips.e6-45.0m*Ellips.e8)/120.0m;
decimal kD = (1237.0m * Ellips.e8)/1260.0m;
 
decimal dL = p.L - Lon;
decimal bp = p.B - DecMath.Sin(p.B) * DecMath.Cos(p.B) * (kA + kB * DecMath.DecimalPow(DecMath.Sin(p.B),2) + kC * DecMath.DecimalPow(DecMath.Sin(p.B),4) + kD * DecMath.DecimalPow(DecMath.Sin(p.B),6));
decimal xp = DecMath.Arctan(DecMath.Tan(bp)/DecMath.Cos(dL));
decimal yp = DecMath.Arctanh(DecMath.Cos(bp)*DecMath.Sin(dL));
decimal x = Scale * Ellips.R * (xp + bt1 * DecMath.Sin(2.0m * xp) * DecMath.Cosh(2.0m * yp) + bt2 * DecMath.Sin(4.0m * xp) * DecMath.Cosh(4.0m * yp) + bt3 * DecMath.Sin(6.0m * xp) * DecMath.Cosh(6.0m * yp) + bt4 * DecMath.Sin(8.0m * xp) * DecMath.Cosh(8.0m * yp)) - Lat2X(Lat);
decimal y = Scale * Ellips.R * (yp + bt1 * DecMath.Cos(2.0m * xp) * DecMath.Sinh(2.0m * yp) + bt2 * DecMath.Cos(4.0m * xp) * DecMath.Sinh(4.0m * yp) + bt3 * DecMath.Cos(6.0m * xp) * DecMath.Sinh(6.0m * yp) + bt4 * DecMath.Cos(8.0m * xp) * DecMath.Sinh(8.0m * yp));
x += Dn; y+=De;
PtXY res = new PtXY(DecMath.Round(x, 4), DecMath.Round(y, 4));
return res;
}
public decimal Lat2X(decimal B) {
decimal bt1 = 0.5m * Ellips.n1 - 2.0m/3.0m * Ellips.n2 + 5.0m/16.0m * Ellips.n3 + 41.0m / 180.0m * Ellips.n4;
decimal bt2 = 13.0m / 48.0m * Ellips.n2 -3.0m/5.0m*Ellips.n3 + 557.0m/1440.0m*Ellips.n4;
decimal bt3 = 61.0m/240.0m*Ellips.n3 - 103.0m/140.0m*Ellips.n4;
decimal bt4 = 49561.0m/161280.0m*Ellips.n4;
decimal kA = Ellips.e2;
decimal kB = (5.0m * Ellips.e4 - Ellips.e6)/6.0m;
decimal kC = (104.0m * Ellips.e6-45.0m*Ellips.e8)/120.0m;
decimal kD = (1237.0m * Ellips.e8)/1260.0m;
decimal dL = 0;
decimal bp = B - DecMath.Sin(B) * DecMath.Cos(B) * (kA + kB * DecMath.DecimalPow(DecMath.Sin(B),2) + kC * DecMath.DecimalPow(DecMath.Sin(B),4) + kD * DecMath.DecimalPow(DecMath.Sin(B),6));
decimal xp = DecMath.Arctan(DecMath.Tan(bp)/DecMath.Cos(dL));
decimal yp = DecMath.Arctanh(DecMath.Cos(bp)*DecMath.Sin(dL));
decimal x = Scale * Ellips.R * (xp + bt1 * DecMath.Sin(2.0m * xp) * DecMath.Cosh(2.0m * yp) + bt2 * DecMath.Sin(4.0m * xp) * DecMath.Cosh(4.0m * yp) + bt3 * DecMath.Sin(6.0m * xp) * DecMath.Cosh(6.0m * yp) + bt4 * DecMath.Sin(8.0m * xp) * DecMath.Cosh(8.0m * yp));
return DecMath.Round(x,20);
}
}
}

Код:

namespace Coord
{
public class Ellipsoid
{
public decimal a;//Большая полуось
public decimal b;//малая полуось
public decimal f;//коэфициент сжатия
public decimal n1;//Третье сжатие (n)
public decimal n2;//Третье сжатие (n) вторая степень
public decimal n3;//Третье сжатие (n) третья степень
public decimal n4;//Третье сжатие (n) четвертая степень
public decimal n5;//Третье сжатие (n) пятая степень
public decimal n6;//Третье сжатие (n) шестая степень
public decimal n7;//Третье сжатие (n) седьмая степень
public decimal e1;//первый экцентриситет
public decimal e2;//первый экцентриситет вторая степень
public decimal e3;//первый экцентриситет третья степень
public decimal e4;//первый экцентриситет четвертая степень
public decimal e5;//первый экцентриситет пятая степень
public decimal e6;//первый экцентриситет шестая степень
public decimal e7;//первый экцентриситет седьмая степень
public decimal e8;//первый экцентриситет восьмая степень
public decimal R;//меридиональный изопериметрический радиус
public Ellipsoid(decimal a, decimal f) {
this.a = a; 
this.f = 1.0m / f;
this.b = a * (1.0m - this.f);
//this.n1 = (this.a -this.b)/(this.a+this.b);
this.n1 = this.f / (2.0m - this.f);
this.n2 = DecMath.DecimalPow(this.n1, 2);
this.n3 = DecMath.DecimalPow(this.n1, 3);
this.n4 = DecMath.DecimalPow(this.n1, 4);
this.n5 = DecMath.DecimalPow(this.n1, 5);
this.n6 = DecMath.DecimalPow(this.n1, 6);
this.n7 = DecMath.DecimalPow(this.n1, 7);
this.e2 = 1.0m - (DecMath.DecimalPow(this.b,2)/ DecMath.DecimalPow(this.a, 2));
this.e1 = DecMath.Sqrt(this.e2);
this.e3 = DecMath.DecimalPow(this.e1, 3);
this.e4 = DecMath.DecimalPow(this.e1, 4);
this.e5 = DecMath.DecimalPow(this.e1, 5);
this.e6 = DecMath.DecimalPow(this.e1, 6);
this.e7 = DecMath.DecimalPow(this.e1, 7); 
this.e8 = DecMath.DecimalPow(this.e1, 8); 
this.R = this.a / (1.0m + this.n1) * (1.0m + 0.25m * this.n2 + 0.05625m * this.n4);
}
public PtXYZ BL2XYZ(PtBL point) {
decimal N = a / DecMath.Sqrt(1.0m - e2 * DecMath.DecimalPow(DecMath.Sin(point.B),2));
decimal H = 0.0m;
decimal x = (N+H)*DecMath.Cos(point.B)*DecMath.Cos(point.L);
decimal y = (N+H)*DecMath.Cos(point.B)*DecMath.Sin(point.L);
decimal z = ((1.0m-e2)*N+H)*DecMath.Sin(point.B);
PtXYZ res = new PtXYZ(x,y,z);
return res;
}
public PtBL XYZ2BL(PtXYZ point) {
// Начальное приближение широты и долготы
decimal longitude = DecMath.Arctan2(point.Y, point.X);
decimal p = DecMath.Sqrt(point.X * point.X + point.Y * point.Y);
decimal e2Prime = e2 / (1 - e2);
decimal theta = DecMath.Arctan2(point.Z, p * (1 - e2));
decimal latitude = DecMath.Arctan((point.Z + e2Prime * a * DecMath.Sin(theta) * DecMath.Sin(theta) * DecMath.Sin(theta)) / (p - e2 * a * DecMath.Cos(theta) * DecMath.Cos(theta) * DecMath.Cos(theta)));
// Итерационное уточнение широты
decimal prevLatitude = latitude;
for(int i = 0; i < 7; ++i) { // Ограничим количество итераций для избежания бесконечного цикла
decimal sinLatitude = DecMath.Sin(latitude);
decimal N = a / DecMath.Sqrt(1 - e2 * sinLatitude * sinLatitude);
decimal h = p / DecMath.Cos(latitude) - N; // Высота над эллипсоидом, предполагается равной 0 для упрощения
// Уточнение широты
decimal nextLatitude = DecMath.Arctan((point.Z / p) / (1 - e2 * N / (N + h)));
if(DecMath.Abs(nextLatitude - prevLatitude) < 1e-10m) {
break;
}
prevLatitude = latitude;
latitude = nextLatitude;
}
decimal B = DecMath.Dec(latitude);
decimal L = DecMath.Dec(longitude);
return new PtBL(B, L);
}
}
}

Код:

public class DecMath
{
public const decimal PI = 3.141592653589793238462643383m;
public static decimal Sign(decimal x) {
if(x > 0.0m) {
return 1;
} else if(x == 0) {
return 0;
} else {
return -1;
}
}
public static decimal Round(decimal value, int digits) {
decimal multiplier = DecimalPow(10, digits);
return Math.Round(value * multiplier) / multiplier;
}
public static decimal Abs(decimal dec) {
if (dec < 0.0m) {
return dec * -1.0m;
} else {
return dec;
}
}
public static decimal Rad(decimal dec) {
decimal res = dec * (PI/180);
return res;
}
public static decimal Rad(int gr, int min, decimal sec) {
decimal dec = Dec(gr, min, sec);
decimal res = dec * (PI/180);
return res;
}
public static decimal Dec(decimal rad) {
decimal res = rad * (180/PI);
return res;
}
public static decimal Dec(int gr, int min, decimal sec) {
decimal res = gr + (min/60.0m)+(sec/3600.0m);
return res;
}
public static decimal DecimalPow(decimal baseValue, int exponent) {
 
if(exponent < 0) {
if(baseValue == 0) {
throw new DivideByZeroException("Деление на ноль.");
}
baseValue = 1 / baseValue;
exponent = -exponent;
}
decimal result = 1;
for(int i = 0; i < exponent; i++) {
result *= baseValue;
}
return result;
}
public static decimal Factorial(int n) {
decimal result = 1;
for(int i = 1; i <= n; i++) {
result *= i;
}
return result;
}
 
public static decimal Exp(decimal x) {
const int maxIterations = 100; 
const decimal tolerance = 1e-28m; 
decimal exp = 1m; 
decimal term = 1m; 
 
for(int i = 1; i < maxIterations; i++) {
term *= x / i;
if(Abs(term) < tolerance) break; 
exp += term;
}
 
return exp;
}
public static decimal Sqrt(decimal S, decimal tolerance = 1e-28m) {
if(S < 0) throw new ArgumentOutOfRangeException(nameof(S), "S не может быть отрицательным.");
if(S == 0) return 0;
 
decimal x = S / 2; 
decimal prevX = 0;
 
while(true) {
prevX = x;
x = (x + S / x) / 2;
if(Math.Abs(x - prevX) < tolerance) break;
}
return x;
}
public static decimal Ln(decimal x) {
if(x <= 0) {
throw new ArgumentOutOfRangeException(nameof(x), "Аргумент должен быть больше 0.");
}
 
// Если x = 1, ln(1) = 0 по определению.
if(x == 1) return 0;
 
if(x > 2) {
int powerOfTwo = 0;
while(x > 2) {
x /= 2;
powerOfTwo++;
}
// ln(x*2^n) = ln(x) + n*ln(2)
return Ln(x) + powerOfTwo * Ln(2);
}
 
// Для 0 < x < 2 используем разложение в ряд Тейлора для ln((1+x)/(1-x))
decimal y = (x - 1) / (x + 1);
decimal ySq = y * y;
decimal sum = 0;
decimal term = y;
int n = 1;
 
while(term != 0) {
sum += term / n;
term *= ySq;
n += 2;
}
 
return 2 * sum;
}
public static decimal NormalizeAngle(decimal angle) {
decimal period = 2 * PI;
angle %= period;
if(angle < 0) {
angle += period;
}
if(angle > PI) {
angle -= period;
}
return angle;
}
public static decimal Cos(decimal x) {
x = NormalizeAngle(x); 
const int maxIterations = 100; 
const decimal tolerance = 1e-28m; 
decimal cos = 1m; 
decimal term = 1m; 
decimal xSquared = x * x; 
 
for(int i = 1; i <= maxIterations; i += 2) {
term *= -xSquared / (i * (i + 1)); 
cos += term;
if(Abs(term) < tolerance) break; 
}
 
return cos;
}
public static decimal Sin(decimal x) {
x = NormalizeAngle(x); 
const int maxIterations = 100; 
const decimal tolerance = 1e-28m; 
decimal sin = 0m;
decimal term = x; 
decimal xSquared = x * x; 
 
for(int i = 1; i <= maxIterations; i += 2) {
sin += term;
term *= -xSquared / ((i + 1) * (i + 2));
if(Abs(term) < tolerance) break; 
}
 
return sin;
}
public static decimal Tan(decimal x) {
decimal cosX = Cos(x);
if(cosX == 0) {
throw new DivideByZeroException("Косинус угла равен 0, тангенс не определен.");
}
return Sin(x) / cosX;
}
public static decimal Cosh(decimal x) {
const int maxIterations = 100;
const decimal tolerance = 1e-28m;
decimal cosh = 1m; // cosh(0) = 1
decimal term = 1m; 
decimal xSquared = x * x;
 
for(int i = 2; i <= maxIterations; i += 2) {
term *= xSquared / ((i - 1) * i);
cosh += term;
if(Abs(term) < tolerance) break;
}
 
return cosh;
}
public static decimal Sinh(decimal x) {
const int maxIterations = 100;
const decimal tolerance = 1e-28m;
decimal sinh = x; 
decimal term = x;
decimal xSquared = x * x;
 
for(int i = 3; i <= maxIterations; i += 2) {
term *= xSquared / ((i - 2) * (i - 1));
sinh += term;
if(Abs(term) < tolerance) break;
}
 
return sinh;
}
public static decimal Arctan(decimal x) {
if(x > 1) {
return PI / 2 - Arctan(1 / x);
} else if(x < -1) {
return -PI / 2 - Arctan(1 / x);
}
const int maxIterations = 100; 
const decimal tolerance = 1e-28m; 
decimal arctan = x; 
decimal term = x;
decimal xSquared = x * x;
 
for(int i = 1; i <= maxIterations; i++) {
term *= -xSquared; 
decimal delta = term / (2 * i + 1);
if(Abs(delta) < tolerance) break; 
arctan += delta;
}
 
return arctan;
}
public static decimal Arctan2(decimal y, decimal x) {
if(x > 0) {
return Arctan(y / x);
}
if(x < 0 && y >= 0) {
return Arctan(y / x) + PI;
}
if(x < 0 && y < 0) {
return Arctan(y / x) - PI;
}
if(x == 0 && y > 0) {
return PI / 2;
}
if(x == 0 && y < 0) {
return -PI / 2;
}
if(x == 0 && y == 0) {
 
return 0;
}
 
return 0;
}
public static decimal Arctanh(decimal x) {
if(x <= -1 || x >= 1) {
throw new ArgumentOutOfRangeException(nameof(x), "Аргумент должен быть в интервале (-1, 1)");
}
 
const int maxIterations = 100;
const decimal tolerance = 1e-28m;
decimal arctanh = 0m;
decimal term = x;
decimal xSquared = x * x;
 
for(int i = 1; i <= maxIterations; i += 2) {
if(i > 1) term *= xSquared * (i - 2) / i;
arctanh += term / i;
if(Abs(term / i) < tolerance) break;
}
 
return arctanh;
}
public static decimal Arcsin(decimal x) {
if(x < -1 || x > 1) {
throw new ArgumentOutOfRangeException(nameof(x), "Аргумент должен быть в интервале [-1, 1].");
}
 
if(x == -1) return -PI / 2;
if(x == 1) return PI / 2;
if(x == 0) return 0;
 
decimal denominator = Sqrt(1.0m - x * x);
if(denominator == 0) {
throw new DivideByZeroException("Деление на ноль в вычислении Arcsin.");
}
decimal res = Arctan(x / denominator);
return Arctan(x / denominator);
}
}

Единственное там не учтен: Угол поворота осей координат в местной системе координат. Хотя поле для хранения поворота в радианах записал.