Добрый день, уважаемые форумчане. Не смог сориентироваться в какой категории создавать тему, надеюсь плюс-минус угадал. Столкнулся с проблемой преобразования координат из плоских топоцентрических МГГТ в геоцентрическую 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 и т.д. Нет проекции Гаусса-Крюгера, а эти преобразования они выполняют по проекции Меркатора. Что тоже меня немного озадачивает. На основании этого, я предполагаю что ключи на геобридж используются вычисленные, но с геодезических координат полученных с проекции Меркатора а не Гаусса-Крюгера, и в этом то и вся моя боль. Если кто то владеет информацией, ключами, советом, теорией, Буду очень благодарен.
Большое спасибо. Скачал на вашем github. В данном документе представлена реализация преобразования из гаусса крюгера для СК, а не МСК. Приращение долготы близкое к нулю я смог получить при значении ординаты:2375591,790874. И на первый взгляд кажется что у вас представлены чуть упрощённые формулы. В принципе используемые мной формулы представлены в первом посте в виде вложений изображений. Возможно, я как то не правильно вас понял, но к сожалению данная реализация мне не помогает :( Большое Вам спасибо, за то, что попытались помочь. Оффтоп (Move your mouse to the spoiler area to reveal the content) Мой расчет: Код: 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 т.е. по широте и долготе за каждый километр невязка увеличивается на пару метров.
Фёдоров Роман, https://geodesist.ru/threads/programmirovanie-preobrazovanija-koordinat.87630/ --- Сообщения объединены, 28 фев 2024, Оригинальное время сообщения: 28 фев 2024 --- https://geodesist.ru/threads/programmirovanie-perevoda-koordinat-iz-odnoj-sistemy-v-druguju.41639/ https://geodesist.ru/threads/programmirovanie-pereschjota-koordinat.22203/
Не заметил как то ваше сообщение ранее... Всё подзабыто конечно. Геодезию люблю, но время всегда было забито тривиальными задачами, не связанные не с математикой, не с геодезией. Прикладная геодезия она такая, заставляет деградировать что ли -) Спасибо большое, тут много полезной информации!
Вообще на сайте сети мггт все разжовано; как уйти от секретности производства полевых работ, и как и как вернуть материалы заказчику в нужной СК.
Плохо. Эллипсоид - 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); } } Единственное там не учтен: Угол поворота осей координат в местной системе координат. Хотя поле для хранения поворота в радианах записал.
