Программирование перевода координат из одной системы в другую.

Вопрос снят. Забыл прописать вычисление арктангенса для дальнейшего расчета, исправил и ВСЕ ЗАРАБОТАЛО!!! stout спасибо за поддержку!!!
Что можно сделать для повышения точности алгоритма, помимо уточнения параметров перехода, или параметры Решают Все?

 

Всем доброго времени суток. Стыдно... спрашивать... но факт есть факт... от прикладных программ мозг отупел и одурнел.
Задача стоит в следующем:
Из WGS84 перейти в МСК-ХХ
с формулами WGS84 в СК42 и др. проблем нет. СК42 в плоские тоже нет. С нахождение расстояний по алгоритму Винсенте, нет.
МСК-ХХ в МСК-Х1Х1 проблем нет.
Нужна подсказка и намек.
есть ключ перехода 5 параметров, начальный осевого меридиана, масштаб, смешение север, смешение запад, и др...
к примеру: 109.03333333333, 0, 1, 1250000, -5111057.628 для Чукотки. Все это я так понимаю на геоиде Красовского
То есть по формулам перехожу из WGS84 пересчитываю в геодезические координаты референц-эллипсоид...
А вот далее что? Перейти в плоские СК42? или в 63 и их трансформировать по этим параметрам, Или на основе параметров перехода рассчитать плоские примерно как это происходит по госту для плоских ск42? Если не трудно приведите формулы, или ссылки. Но только не советуйте МапИнфо транскор и др. продукты. Мне не нужно переводит координаты мне нужен алгоритм. Сильно не стыдите, понимаю что это базовые знания в геодезии, но за 12 лет забыл.

 

???

Нет.

Да.
То, что вы привели:

это и есть параметры проекции Гаусса-Крюгера для

Формулы можно и из ГОСТа, можно из Coordinate Conversions and Transformations including Formulas Revised - April 2015.
Самые простые (компактные) и достаточно точные формулы даны в Transverse Mercator: Bowring series
Посмотрите ссылки в сообщении
http://geodesist.ru/forum/threads/Преобразование-координат-из-одной-плоской-системы-в-другую.42027/#post-479425

А это-то вам зачем?
Остаётся только «Главный вопрос жизни, вселенной и всего такого».
Зачем вам это надо? Для какой такой задачи понадобился этот переход?
А также, откуда у вас точные параметры связи из

? Как и кто это определял?

Удивительное рядом.

 

Это было банальное утверждение, что это все параметры перехода проекции Гаусса-Крюгера на геоиде с параметрами Красовского.

Изначально думал обойтись этим алгоритмом для вычисления плоских координат.
Далее приведу не верные рассуждения свои.
Х= это ни что иное как расстояние от экватора до точки по эллипсоиду. По этому рассчитываем длину дуги по одноименной долготе до нужной широты.
Y= это расстояние от начальной широты (для зоны) до точки на заданной широте.
сейчас использую ГОСТ 32453-2013 п.5.4.

Я не говорил, что у меня точные параметры. У меня параметры из ГОСТ 32453-2013. Переход через WGS84 в ПЗ-90.02 далее в СК42

да сам сижу и думаю, что имел в виду. Главное что проблем у меня с этим нет

Полевой контроллер TOPCON FC-250 стоит 130т.р. или сотовый на Android 7т.р. Да и работает быстрее, разница только в IP... Исходя из этих моментов захотелось поупражняться и захотелось написать на Java, что то похожее на Magnet Fild.
В общем кто после работы бухает, кто то играет в танчики-самолетики, стрелялки. А я маюсь ерундой.

 

Оставим в стороне рениксы типа "геоид Красовского" и

Ближе к теме.

Во втором предложении непонятно зачем говорить о долготе, от которой длина дуги ну никак не зависит. А первое предложение справедливо лишь частично. Оно справедливо только для точек на осевом меридиане. Длину дуги меридиана S от экватора до параллели с геодезической широтой B можно записать как

Интересно, что строгая формула для проекции Гаусса-Крюгера выглядит аналогично

Вот только φ здесь не геодезическая широта, а комплекснозначная величина, одна из схем вычисления которой выглядит так:
вычисляем изометрическую широту q

формируем комплексное число q + iΔλ, где Δλ — удаление от осевого меридиана
и решаем относительно sinφ нелинейное уравнение

В качестве начального приближения для sinφ можно принять

Вообще-то, умножение матриц обладает свойством ассоциативности, поэтому перемножив один раз матрицы (и запомнив/сохранив результат), можно напрямую переходить от WGS 84 к СК-42.

Тогда посмотрите C++ and Java Code for Geodesic and Meridian Arc Computations -- by Klaus Hehl.

 

Оффтоп

Соседняя тема по непонятным причинам называцца "Помогите перевести координаты участка из МСК в WGS-84 с точностью 10-15м. "
Есть предположение, что от головы до хвоста (МСК - WGS) и от хвоста до головы (WGS - МСК) примерно одинаково. 10-15 метров.
Этот форум, вроде как, геодезический. ГОСТы, Формулы, однако, хорошо. Но что на выходе? 10-15 метров?
Вместе с тем, есть NTv2. По непонятным причинам невидимый для банальных обывателей на сайте gisa (http://www.gisa.ru/93691.html)
Есть и обкатка NTv2 (http://gis-lab.info/forum/viewtopic.php?f=2&t=13194&start=45).
Интересно, 40 сантиметров лучше чем 10-15 метров?

 

Оффтоп

Родичкин, та тема специально так названа как раз для "банальных обывателей", но уже этого сайта.
Распространение GPS модулей в телефоны и планшеты привело к тому, что получив перевод из МСК люди сами на местности пытаются свои участки выносить.
А что, айфон за 60 т.р. круче геодезиста за 10-15 т.р., да и зачем тратиться.
С последствиями работы этих "умников" приходится сталкиваться уже нормальным геодезиста.

 

Наличие развитой индустрии моделирования нисколько не мешает любителям клепать модельки чего-либо своими руками. Им это просто интересно.

Для тех, у кого вместо карты пачка Беломора, этот вопрос не имеет смысла.

 

Оффтоп

С последствиями работы этих "ГОСТов" приходится сталкиваться относительно нормальным геодезистам ...

Плюс к этому у нас позволительно:
- Вводить в обиход ГСК-2011, которая по факту является (будет являться) очередной за СК-95 реализацией СК-42;
- выпускать в свет новый ГОСТ по пересчету координат, который застрял на WGS_се дрейвейшей реализации, в котором не наблюдается ГСК-2011 и ПЗ-90.11, и в котором отсутствуют поля поправок СК-42 --> СК-95 --> ГСК-2011 и поля поправок к моделям геоидов EGM2008 и EIGEN-4 ....

Всё это , блин, называцца - ФЦП ГЛОНАСС ... в натуре.

 

Оффтоп

Но ведь для этого нужно проделать много работы, в т.ч. и полевых измерений.
И, что важнее, не засекретить её результаты.

Давал бы ГОСТ сантиметровую точность, не было бы этих проблем и этих тем.
99 % геодезистов его тоже хватало бы "за глаза".

 

А на кой .... тогда существуют РОСРЕЕСТР и РОСКАРТОГРАФИЯ в виде единственного исполнителя?
Или они привыкли распихивать государственное бабло по собственным карманам без "проделывания много работы"?
Типа, натянем СК-42 на ФАГС и СГС, на выходе получим ГСК-2011 ??? А по полям побегать?

Отмывка бабла на фоне геобреда первична, получаемый SXFовский бредопродукт - вторичен ...
Правильной дорогой идете господа квазигеодезисты ...

 

Оффтоп

Демаго́г (др.-греч. δημαγωγός) — демократ и диктатор в Древней Греции; также популист, «народный» политик. Первоначально слово не имело негативного оттенка и обозначало то, что впоследствии Аристотель (быть может, из-за дискредитации термина «демагог») передавал через выражение «простат (защитник, представитель интересов) народа». «Простатами народа», то есть демократическими лидерами, на протяжении большей части V в. до н. э. были выходцы из знатных родов, вроде Фемистокла или Перикла. Положение меняется к концу столетия, когда на авансцену политической жизни выходят незнатные «выскочки», вроде владельца кожевенной мастерской Клеона или владельца мастерской ламп Гипербола, с радикальными политическими устремлениями (в современной историографии их зовут «вождями радикальной демократии»). Противники обвиняли их в популизме, политической безответственности, коррупции и игре на самых низких и темных инстинктах толпы. Благодаря им, понятие «демагог» начинает обозначать политика-популиста и приближается к современному значению. В смысле современного термина «популист», оно использовалось ещё в XIX веке, например, для обозначения революционных лидеров. Более того, «Словарь иностранных слов» 1954 г. определяет термин «демагог» как
«политикан, лицо, старающееся создать себе популярность среди народных масс недостойными средствами (извращением фактов, лестью и т. д.)», то есть вполне в древнегреческом смысле этого слова.
Разительная и не утратившая актуальности сатира на политическую демагогию дана в комедии Аристофана «Всадники», направленной лично против Клеона. Клеон выведен в ней под именем Пафлагонца, раба выжившего из ума старика Демоса («Народа»), обманывающего и обирающего своего господина.
В. И. Даль уже так определяет демагога — «крайний демократ, добивающийся власти во имя народа; тайный возмутитель; поборник безначалия, желающий ниспровергнуть порядок управления», а демагогию как «господство власти народа, черни в управлении»

По теме:
Сеточные поправки, несомненно, позволяют получить более точный результат. Описание форматов можно найти в Инете.
Причина, по которой нет РукДоков и числовых значений поправок на просторы РФ - мне достоверно не известна.
Однако, единообразие пересчета с их использованием, предполагает их официальное опубликование РОСРЕЕСТРом.

 

stout, огромное тебе СПАСИБО. Что не сообщение, то просто кладезь информации. Очень понравился документ
IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015 To facilitate improvement, this document is subject to revision. The current version is available at www.epsg.org. Как там все грамотно расписано и понятно объяснено! Изучаем дальше =)

кусок рабочего кода метода на Java (раскрыть) кусок рабочего кода метода на Java (свернуть)

// http://www.iogp.org/pubs/373-07-2.pdf
// проекция трансмеркатора
public static double[] TMERC (double B, double L, double H,
double EGD[], double Key[]) {
double π=Math.PI;
double a=EGD[0];// большая полуось
double f=EGD[1];// полярное сжатие
double e2 = f* 2 - f * f;// Квадрат первого эксцентриситета элипсоида
double b=a-f*a;// получаем маоую полуось из большой полуоси и полярного сжатия
double eh2=(a*a-b*b)/(b*b);// Квадрат второго эксцентриситета элипсоида
//λo ϕo - требуеться перевести в радианы
// a=6377563.396; //для проверки потом удалить
// f=1/299.32496;//для проверки потом удалить
/*λo,ϕo - изначально передаеться углы в формате ГГММСС,СССС
* GeodeZ.trGGMMSStoGG(L) переводит в вид ---> ГГ.ГГГГГГГГ
*/
double λo=GeodeZ.trGGMMSStoGG(Key[0]); //Осевой меридиан (Longitude of natural origin) Central Meridian
double ϕo=GeodeZ.trGGMMSStoGG(Key[2]);// Начальная широта (Latitude of natural origin) Lat0
double ko=Key[1]; // масштаб (Scale factor at natural origin) Scale
double FE =Key[3];//Смешение Восток (False easting) East0
double FN =Key[4];// Смешение север (False northing) North0
// _________________________Данные для контроля и проверки результата
λo=Math.toRadians(λo); //Осевой меридиан (Longitude of natural origin) Central Meridian
ϕo=Math.toRadians(ϕo);// Начальная широта (Latitude of natural origin) Lat0
// ko=0.9996012717; // масштаб (Scale factor at natural origin) Scale
// FE =400000.00;//Смешение Восток (False easting) East0
//FN =-100000.00;// Смешение север (False northing) North0
//_________________________________

double ϕ= Math.toRadians(B);
double λ=Math.toRadians(L);
// System.out.println ("Квадрат первого эксцентриситета элипсоида="+e2);
//System.out.println ("Квадрат второго эксцентриситета элипсоида="+eh2);
// Вычисление констант
double n=f/(2-f);
b=(a/(1+n))*(1+n*n/4.0+Math.pow(n, 4)/64.00);//заменяем переменную b на другое значение константы для расчета
double n2=Math.pow(n, 2);
double n3=Math.pow(n, 3);
double n4=Math.pow(n, 4);
double h1 = n / 2.0 - (2.0 / 3.0) * n2 + (5.0 / 16.0) * n3 + (41.0 / 180.0) * n4;
double h2 = (13.0 / 48.0) * n2 - (3.0 / 5.0) * n3 + (557.0 / 1440.0) * n4;
double h3 = (61.0 / 240.0) * n3 - (103.0 / 140.0) * n4;
double h4 = (49561.0/161280.0) * n4;
System.out.println ("константы"+" n="+n+" B="+b+" h1="+h1+" h2="+h2+" h3="+h3+" h4="+h4);
double Mo=0;
double Qo=0;
double βo=0;
double ξO0=0;
double ξO1=0;
double ξO2=0;
double ξO3=0;
double ξO4=0;
double ξO=0;

if (ϕo == 0){Mo=0;} else{if (ϕo==π/2){Mo=b*(π/2);} else {if (ϕo==(-π/2)){Mo=b*(-π/2);}else {
double tanϕo=Math.tan(ϕo);
double sinϕo=Math.sin(ϕo);
double e=Math.sqrt(e2);
double e_sinϕo=e*sinϕo;
Qo=Math.log(tanϕo+Math.sqrt(tanϕo*tanϕo+1))-(e*Math.log(Math.sqrt(1-e_sinϕo*e_sinϕo)/(1-e_sinϕo)));
βo=Math.atan(Math.sinh(Qo));
ξO0=βo;
ξO1 = h1* Math.sin(2.0*ξO0);
ξO2 = h2* Math.sin(4.0*ξO0);
ξO3 = h3*Math.sin(6.0*ξO0);
ξO4 = h4*Math.sin(8.0*ξO0);
ξO = ξO0+ ξO1+ ξO2+ ξO3+ ξO4;
Mo = b *ξO;

// System.out.println ("Qo="+Qo+" βo="+βo +" Mo="+Mo);
// System.out.println("ξO1="+ξO1);
// System.out.println("ξO2="+ξO2);
// System.out.println("ξO3="+ξO3);
// System.out.println("ξO4="+ξO4);
} } }

/*
Mo = a[(1 – e2/4 – 3e4 /64 – 5e6 /256 –....)ϕO – (3e2 /8 + 3e4 /32 +
45e6 /1024+....)sin2ϕO + (15e4 /256 + 45e6 /1024 +.....)sin4ϕO –
(35e6 /3072 + ....)sin6ϕO + .....]
*/

double tanϕ=Math.tan(ϕ);
double sinϕ=Math.sin(ϕ);
double e=Math.pow(e2,0.5);
double e_sinϕ=e*sinϕ;
double Q=Math.log(tanϕ+Math.sqrt(tanϕ*tanϕ+1))-(e*Math.log(Math.sqrt(1-e_sinϕ*e_sinϕ)/(1-e_sinϕ)));
double β=Math.atan(Math.sinh(Q));
double cosβsinDλ=(Math.cos(β)*Math.sin(λ-λo));
double η0=Math.log(Math.sqrt(1-cosβsinDλ*cosβsinDλ)/(1-cosβsinDλ));
double ξ0 = Math.asin (Math.sin (β)* Math.cosh (η0));
double ξ1 = h1 * Math.sin(2 * ξ0) * Math.cosh(2 * η0);
double η1 = h1 * Math.cos(2 * ξ0) * Math.sinh(2 * η0);
double ξ2 = h2 * Math.sin(4 * ξ0) * Math.cosh(4 * η0);
double η2 = h2 * Math.cos(4 * ξ0) * Math.sinh(4 * η0);
double ξ3 = h3 * Math.sin(6 * ξ0) * Math.cosh(6 * η0);
double η3 = h3 * Math.cos(6 * ξ0) * Math.sinh(6 * η0);
double ξ4 = h4 * Math.sin(8 * ξ0) * Math.cosh(8 * η0);
double η4 = h4 * Math.cos(8 * ξ0) * Math.sinh(8 * η0);
double ξ = ξ0 + ξ1 + ξ2 + ξ3 + ξ4;
double η = η0 + η1 + η2 + η3 + η4;
double XY [] = {(FN + ko* (b* ξ-Mo)),(FE + ko* b *η)};
return XY ;}

 

Так и делаем.

 

Родичкин, съездите в отпуск, например в Крым. Тут Роскартография длинные щупальца распустила. Полным ходом идёт обследование пунктов, спутниковые наблюдения на них. Программа предусм. перенаблюдать все пункты 1-2 классов. И это будет сделано. Какая для вас к чёрту разница сколько будет стоить?! Что вы заладили!
-=13=-, всё будет хо-ро-шо! Работа идёт!

 

Коли тема так называться. Привожу пример рабочих, проверенных методов преобразования координат на java. Может пригодиться кому. Код не причесан, отладочные ремарки не убраны. Очень много ошибок по Русскому языку в ремарках. Прошу терпимо отнестись к этому.
Метод trnsfrm и transform - ГОСТ 32453-2013 5.3 Преобразование геодезических координат. Первый дает Дельты, второй готовый результат
Метод SK42_WGS84 выдает параметры перехода от СК-42 в WGS84 (один из вариантов) для тестов, а так это все храниться в XML
Методы EGD......... выдает параметры эллипсоидов. (опять жен для отладки) в сборке работает по другому =)
Методы MSK......... выдает параметры ключей перехода проекции
Метод TMERC - переводит из геодезических координат в координаты на плоскости проекции Тарнс Меркатора. по ключу перехода переводит как туда так и обратно.
Методы arcth, asinh - гиперболические функции недостающие в библиотеке Math.
Метод trGGMMSStoGG трансформирует данные из вида ГГММСС.ССС в ГГ.ГГГГ (данные XML представлены в таком виде)
Метод PloscToCk42ug - трансформирует ГОСТ 32453-2013 5.4 преобразует плоские прямоугольные в геодезические СК42
Метод ck42uglToPlosc -трансформирует ГОСТ 32453-2013 5.4 преобразует геодезические плоские прямоугольные в СК42
Метод returnGMS вернет строку вида GG°MM′SS″
_________________
Что бы запустить нужно написать отдельно метод main, но это уже основы Java.
Методы все проверенны и обкатаны. Но еще не оптимизированы. Ошибок нет. Но есть оговорки. О них можно подробно посмотреть в ссылках в ремарках внутри тела программы.

Java код

Код:

public class GeodeZ {
// трансформация геодезических координат из системы А в B
// Yes = true из A B, Yes = false из B A
public static double[] transform(boolean Yes, double B, double L, double H,
double p[]) {
double[] G = GeodeZ.trnsfrm(Yes, B, L, H, p[0], p[1], p[2], p[3],
p[4], p[5], p[6], p[7], p[8], p[9], p[10], p[11]);
// B = (G[0] / 3600 + B * 2) / 2;
// L = (G[1] / 3600 + L * 2) / 2;
// H = (G[2] + H * 2) / 2;
 
// G = GeodeZ.trnsfrm(Yes, B, L, H, p[0], p[1], p[2], p[3], p[4], p[5],
// p[6], p[7], p[8], p[9], p[10], p[11]);
G[0] = B + G[0] / 3600;
G[1] = L + G[1] / 3600;
G[2] = H + G[2];
return G;
 
public static double[] trnsfrm(boolean Yes, // Yes = true из A B, Yes = false из B A
double B,// Широта в системе А
double L,// Долгота в системе А
double H,// Высота в системе А
double Aa,// Большая полуось А
double Ea,// Сжатие элепсоида А
double Ab,// Большая полуось B
double Eb,// Сжатие элепсоида B
double Dx,// параметры геоцентрических смещений по оси X
double Dy,// параметры геоцентрических смещений по оси Y
double Dz,// параметры геоцентрических смещений по оси Z
double Wx,// параметры геоцентрических поворота по оси X
double Wy,// параметры геоцентрических поворота по оси Y
double Wz,// параметры геоцентрических поворота по оси Z
double Mm,// масштабны коофициент
double P// = 206264.8062; ro- Число угловых секунд в радиане
) {
// из A в B если требуеться в другую сторону, то Yes=false
double[] dBdLdH = { 0, 0, 0 };
 
 
 
double BRadi = Math.toRadians(B);
double LRadi = Math.toRadians(L);
double SinB = Math.sin(BRadi);
double CosB = Math.cos(BRadi);
double SinL = Math.sin(LRadi);
double CosL = Math.cos(LRadi);
double Cos2B = Math.cos(BRadi * 2);
double TanB = Math.tan(BRadi);
double E2a = Ea * 2 - Ea * Ea;// Квадрат эксцентриситета элипсоида
// системы А
double E2b = Eb * 2 - Eb * Eb;// Квадрат эксцентриситета элипсоида
// системы В
 
double A = (Ab + Aa) / 2;// Средная главная ось
double E2 = (E2b + E2a) / 2;
double De2 = E2b - E2a;
double DA = Ab - Aa;
double M = A * (1 - E2) * Math.pow((1 - E2 * SinB * SinB), (-1.5));
double N = A * Math.pow((1 - E2 * SinB * SinB), (-0.5));
double dB = P
/ (M + H)
* (N / A * E2 * SinB * CosB * DA + ((N * N) / (A * A) + 1) * N
* SinB * CosB * De2 / 2 - (Dx * CosL + Dy * SinL)
* SinB + Dz * CosB) - Wx * SinL * (1 + E2 * Cos2B) + Wy
* CosL * (1 + E2 * Cos2B) - P * Mm * E2 * SinB * CosB;
 
double dL = P / ((N + H) * CosB) * (-Dx * SinL + Dy * CosL) + TanB
* (1 - E2) * (Wx * CosL + Wy * SinL) - Wz;
 
double dH = -A / N * DA + N * SinB * SinB * De2 / 2
+ (Dx * CosL + Dy * SinL) * CosB + Dz * SinB - N * E2 * SinB
* CosB * (Wx / P * SinL - Wy / P * CosL) + ((A * A) / N + H)
* Mm;
double Ink;
if (Yes) {
Ink = 1;
} else {
Ink = -1;
}
dBdLdH[0] = dB * Ink;
dBdLdH[1] = dL * Ink;
dBdLdH[2] = dH * Ink;
return dBdLdH;
}
 
public static double[] SK42_WGS84() {
final double Aa = 6378245;
final double Ea = 1 / 298.3;
final double Ab = 6378136;
final double Eb = 1 / 298.257223563;
final double Dx = 23.9; //23.92;
final double Dy = -141.3;//-141.27
final double Dz = -80.9;
final double Wx = 0;
final double Wy = -0.35;
final double Wz = -0.82;
final double Mm = -0.12*Math.pow(10,-6);
final double P = 206264.8062;
double[] R = { Aa, Ea, Ab, Eb, Dx, Dy, Dz, Wx, Wy, Wz, Mm, P };
 
public static double[] EGDgostPZ90() {
final double A = 6378136;
final double E = 1 / 298.25784;
double P[]={0,0};
P[0]=A;
P[1]=E;
return P;
 
}
public static double[] EGDgostKrosovskogo() {
final double A = 6378245;
final double E = 1 / 298.3;
double P[]={0,0};
P[0]=A;
P[1]=E;
return P;
}
public static double[] EGDgostWGS84() {
final double A = 6378139;
final double E = 1 / 298.257223563;
double P[]={0,0};
P[0]=A;
P[1]=E;
return P;
}
 
 
 
 
public static double[]MSK_23Zone1 ()
{
double []B={0,0,0,0,0};
B[0]=375900.00000; //<P0 name="Central Meridian"
B[1]=1.0000000000;//<P1 name="Scale">1.0000000000</P1>
B[2]=00000.00000; //<P2 name="Lat0">00000.00000</P2>
B[3]=1300000;//<P3 name="East0">1300000</P3>
B[4]=-4511057.628; //<P4 name="North0">-4511057.628</P4>
return B;}
 
 
//ключ перехода
public static double[]MSKtest ()
{
double []B={0,0,0,0,0};
B[0]=-020000.00000; //<P0 name="Central Meridian"
B[1]=0.9996012717;;//<P1 name="Scale">1.0000000000</P1> // ko=0.9996012717; // масштаб (Scale factor at natural origin) Scale
B[2]=490000.00000; //<P2 name="Lat0">00000.00000</P2>
B[3]=400000.00;//<P3 name="East0">1300000</P3> // FE =400000.00;//Смешение Восток (False easting) East0
B[4]=-100000.00; //<P4 name="North0">-4511057.628</P4>//FN =-100000.00;// Смешение север (False northing) North0
return B;}
//параметры элепсоида
public static double[] EGDtest() {
final double A = 6377563.396;
final double E = 1 / 299.32496;
double P[]={0,0};
P[0]=A;
P[1]=E;
return P;
}
 
 
// http://www.iogp.org/pubs/373-07-2.pdf
// проекция трансмеркатора
/*
* пердаем в метод B и L широта долгота в формате ГГ.ГГГГГ либо в B и L в
* координатах Х и Y в ММ.МММММ boolean toPl = true если из геогр.
* координаты в проекцию boolean toPl false если из проекции в гегр.
* координаты EGD[] масив исходных параметров геоида EGD[0]- большая полуось
* EGD[1]- полярное сжатие Key[] массив ключа перехода B[0]=-ГГMMCC.СССС;
* //<P0 name="Central Meridian" B[1]=X.XXXXXX;;//<P1
* name="Scale">1.0000000000</P1> // ko=0.9996012717; // масштаб (Scale
* factor at natural origin) Scale B[2]=ГГММСС.СССС; //<P2
* name="Lat0">00000.00000</P2> B[3]=ММММММ.ММММ;//<P3
* name="East0">1300000</P3> // FE =400000.00;//Смешение Восток (False
* easting) East0 B[4]=ММММММ.ММММ; //<P4
* name="North0">-4511057.628</P4>//FN =-100000.00;// Смешение север (False
* northing) North0
*/
public static double[] TMERC(double B, double L, boolean toPl,
double EGD[], double Key[]) {
double[] XY = { 0, 0 };
double π = Math.PI;
double a = EGD[0];// большая полуось
double f = EGD[1];// полярное сжатие
double e2 = f * 2 - f * f;// Квадрат первого эксцентриситета элипсоида
double b = a - f * a;// получаем маоую полуось из большой полуоси и
// полярного сжатия
double eh2 = (a * a - b * b) / (b * b);// Квадрат второго
// эксцентриситета элипсоида (он
// в принципе здесь не нужен)
double e = Math.sqrt(e2); // первый ксцентриситета элипсоида
// λo ϕo - требуеться перевести в радианы
/*
* λo,ϕo - изначально передаеться углы в формате ГГММСС,СССС
* GeodeZ.trGGMMSStoGG(L) переводит в вид ---> ГГ.ГГГГГГГГ
*/
double λo = GeodeZ.trGGMMSStoGG(Key[0]); // Осевой меридиан (Longitude
// of natural origin)
// Central Meridian
double ϕo = GeodeZ.trGGMMSStoGG(Key[2]);// Начальная широта (Latitude of
// natural origin) Lat0
double ko = Key[1]; // масштаб (Scale factor at natural origin) Scale
double FE = Key[3];// Смешение Восток (False easting) East0
double FN = Key[4];// Смешение север (False northing) North0
// _________________________Данные для контроля и проверки результата
λo = Math.toRadians(λo); // Осевой меридиан (Longitude of natural
// origin) Central Meridian
ϕo = Math.toRadians(ϕo);// Начальная широта (Latitude of natural origin)
// Lat0
// _________________________________
double ϕ = Math.toRadians(B);
double λ = Math.toRadians(L);
// System.out.println ("Квадрат первого эксцентриситета элипсоида="+e2);
// System.out.println
// ("Квадрат второго эксцентриситета элипсоида="+eh2);
// Вычисление констант
double n = f / (2 - f);
/* заменяем переменную b на другое значение константы для расчета */
b = (a / (1 + n)) * (1 + n * n / 4.0 + Math.pow(n, 4) / 64.00);
double n2 = Math.pow(n, 2);
double n3 = Math.pow(n, 3);
double n4 = Math.pow(n, 4);
double h1 = n / 2.0 - (2.0 / 3.0) * n2 + (5.0 / 16.0) * n3
+ (41.0 / 180.0) * n4;
double h2 = (13.0 / 48.0) * n2 - (3.0 / 5.0) * n3 + (557.0 / 1440.0)
* n4;
double h3 = (61.0 / 240.0) * n3 - (103.0 / 140.0) * n4;
double h4 = (49561.0 / 161280.0) * n4;
/*System.out.println("константы" + " n=" + n + " B=" + b + " h1=" + h1
+ " h2=" + h2 + " h3=" + h3 + " h4=" + h4);*/
// в этом месте решил обявить все переменные
double Mo = 0;
double Qo = 0;
double βo = 0;
double ξO0 = 0;
double ξO1 = 0;
double ξO2 = 0;
double ξO3 = 0;
double ξO4 = 0;
double ξO = 0;
// System.out.println("ϕo=" + ϕo);
if (ϕo == 0) {
Mo = 0;
} else {
if (ϕo == π / 2) {
Mo = b * (π / 2);
} else {
if (ϕo == (-π / 2)) {
Mo = b * (-π / 2);
} else {
double tanϕo = Math.tan(ϕo);
double sinϕo = Math.sin(ϕo);
double e_sinϕo = e * sinϕo;
Qo = Math.log(tanϕo + Math.sqrt(tanϕo * tanϕo + 1))
- (e * Math.log(Math.sqrt(1 - e_sinϕo * e_sinϕo)
/ (1 - e_sinϕo)));
βo = Math.atan(Math.sinh(Qo));
ξO0 = βo;
ξO1 = h1 * Math.sin(2.0 * ξO0);
ξO2 = h2 * Math.sin(4.0 * ξO0);
ξO3 = h3 * Math.sin(6.0 * ξO0);
ξO4 = h4 * Math.sin(8.0 * ξO0);
ξO = ξO0 + ξO1 + ξO2 + ξO3 + ξO4;
Mo = b * ξO;
// System.out.println ("Qo="+Qo+" βo="+βo +" Mo="+Mo);
// System.out.println("ξO1="+ξO1);
// System.out.println("ξO2="+ξO2);
// System.out.println("ξO3="+ξO3);
// System.out.println("ξO4="+ξO4);
}
}
}
/*
* Mo = a[(1 – e2/4 – 3e4 /64 – 5e6 /256 –....)ϕO – (3e2 /8 + 3e4 /32 +
* 45e6 /1024+....)sin2ϕO + (15e4 /256 + 45e6 /1024 +.....)sin4ϕO –
* (35e6 /3072 + ....)sin6ϕO + .....]
*/
if (toPl) {
double tanϕ = Math.tan(ϕ);
double sinϕ = Math.sin(ϕ);
double e_sinϕ = e * sinϕ;
double Q = Math.log(tanϕ + Math.sqrt(tanϕ * tanϕ + 1))
- (e * Math.log(Math.sqrt(1 - e_sinϕ * e_sinϕ)
/ (1 - e_sinϕ)));
double β = Math.atan(Math.sinh(Q));
double cosβsinDλ = (Math.cos(β) * Math.sin(λ - λo));
double η0 = Math.log(Math.sqrt(1 - cosβsinDλ * cosβsinDλ)
/ (1 - cosβsinDλ));
double ξ0 = Math.asin(Math.sin(β) * Math.cosh(η0));
double ξ1 = h1 * Math.sin(2 * ξ0) * Math.cosh(2 * η0);
double η1 = h1 * Math.cos(2 * ξ0) * Math.sinh(2 * η0);
double ξ2 = h2 * Math.sin(4 * ξ0) * Math.cosh(4 * η0);
double η2 = h2 * Math.cos(4 * ξ0) * Math.sinh(4 * η0);
double ξ3 = h3 * Math.sin(6 * ξ0) * Math.cosh(6 * η0);
double η3 = h3 * Math.cos(6 * ξ0) * Math.sinh(6 * η0);
double ξ4 = h4 * Math.sin(8 * ξ0) * Math.cosh(8 * η0);
double η4 = h4 * Math.cos(8 * ξ0) * Math.sinh(8 * η0);
double ξ = ξ0 + ξ1 + ξ2 + ξ3 + ξ4;
double η = η0 + η1 + η2 + η3 + η4;
XY[0] = (FN + ko * (b * ξ - Mo));
XY[1] = (FE + ko * b * η);
}
else {
// _________________________________________________________________________________
h1 = n / 2.0 - (2.0 / 3.0) * n2 + (37.0 / 96.0) * n3
- (1.0 / 360.0) * n4;
h2 = (1.0 / 48.0) * n2 + (1.0 / 15.0) * n3 - (437.0 / 1440.0) * n4;
h3 = (17.0 / 480.0) * n3 - (37.0 / 840.0) * n4;
h4 = (4397.0 / 161280.0) * n4;
double η = (L - FE) / (b * ko);
double ξ = ((B - FN) + ko * Mo) / (b * ko);
// System.out.println("____________________________________________");
// System.out.println("L="+L+" B="+B);
// System.out.println("____________________________________________");
// System.out.println("ξ="+ξ+" η="+η);
double ξ1 = h1 * Math.sin(2 * ξ) * Math.cosh(2 * η);
double η1 = h1 * Math.cos(2 * ξ) * Math.sinh(2 * η);
double ξ2 = h2 * Math.sin(4 * ξ) * Math.cosh(4 * η);
double η2 = h2 * Math.cos(4 * ξ) * Math.sinh(4 * η);
double ξ3 = h3 * Math.sin(6 * ξ) * Math.cosh(6 * η);
double η3 = h3 * Math.cos(6 * ξ) * Math.sinh(6 * η);
double ξ4 = h4 * Math.sin(8 * ξ) * Math.cosh(8 * η);
double η4 = h4 * Math.cos(8 * ξ) * Math.sinh(8 * η);
double ξ0 = ξ - (ξ1 + ξ2 + ξ3 + ξ4);
double η0 = η - (η1 + η2 + η3 + η4);
// System.out.println("____________________________________________");
// System.out.println("ξ0="+ξ0+" η0="+η0);
double β = Math.asin(Math.sin(ξ0) / Math.cosh(η0));
// System.out.println("β="+β);
double Q1 = GeodeZ.asinh(Math.tan(β));
// System.out.println("Q1="+Q1);
double Q2 = Q1 + (e * GeodeZ.arcth(e * Math.tanh(Q1)));
for (int i = 0; i < 2; ++i) {
Q2 = Q1 + (e * GeodeZ.arcth(e * Math.tanh(Q2)));
}// выполняем 3 итерации
// System.out.println("Q2="+Q2);
ϕ = Math.atan(Math.sinh(Q2));
// System.out.println ("ϕ="+ ϕ );
λ = λo + Math.asin(Math.tanh(η0) / Math.cos(β));
XY[0] = Math.toDegrees(ϕ);
XY[1] = Math.toDegrees(λ);
}
 
return XY ;}
public static double arcth (double x) { return Math.log((1+x)/(1-x))/2.0;}
public static double asinh (double x) { return Math.log(x+Math.sqrt(x*x+1));}
public static double trGGMMSStoGG (double L) {
double G=(int)L/10000;
double M= (int)((L/10000-G)*100);
double S=(((L/10000)-G)*100-M)*100;
L=G+M/60+S/(60*60);
return L;}
//СК-42 преобразование в геодезических коордиеат в плоских прямоугольных 
public static double[] ck42uglToPlosc(double B, double L, double H) {
double[] Point = { 0, 0, 0 };
//вичисляем номер зоны пока L не перезаписалась в радианах
double n = (int) ((6 + L) / 6);
double l = (L - (3 + 6 * (n - 1))) / 57.29577951;
//
// System.out.println ("Проверка из геод. в плоские");
L = Math.toRadians(L);
//System.out.println("в радианах L="+L);
B = Math.toRadians(B);
double CosB = Math.cos(B);

double SinB = Math.sin(B);
// System.out.println("sinB=" +SinB);
double Sin2_B = Math.pow(SinB, 2);
double Sin4_B = Math.pow(SinB, 4);
double Sin6_B = Math.pow(SinB, 6);
double l2 = Math.pow(l, 2);
Point[0] = 6367558.4968
* B
- Math.sin(B * 2)
* (16002.8900 + 66.9607 * Sin2_B + 0.3515 * Sin4_B - l2
* (1594561.25 + 5336.535 * Sin2_B + 26.790 * Sin4_B
+ 0.149 * Sin6_B + l2
* (672483.4 - 811219.9 * Sin2_B + 5420.0
* Sin4_B - 10.6 * Sin6_B + l2
* (278194 - 830174 * Sin2_B + 572434
* Sin4_B - 16010 * Sin6_B + l2
* (109500 - 574700 * Sin2_B
+ 863700 * Sin4_B - 398600 * Sin6_B)))));
Point[1] = (5 + 10 * n)
* 100000.0
+ l
* CosB
* (6378245 + 21346.1415 * Sin2_B + 107.1590 * Sin4_B + 0.5977
* Sin6_B + l2
* (1070204.16 - 2136826.66 * Sin2_B + 17.98 * Sin4_B
- 11.99 * Sin6_B + l2
* (270806 - 1523417 * Sin2_B + 1327645 * Sin4_B
- 21701 * Sin4_B + l2
* (79690 - 866190 * Sin2_B + 1730360
* Sin4_B - 945460 * Sin6_B))));
return Point;
}
//СК-42 плоских прямоугольных преобразование в геодезические коордиеаты 
// ГОСТ 32453-2013 5.4.1
public static double[] PloscToCk42ug (double X, double Y, double H) {
double[] Point = { 0, 0, 0 };
double β = X / 6367558.4968;
double Sin2b = Math.sin(β * 2);
double Sin2_b = Math.pow(Math.sin(β), 2);
double Bo = β
+ Sin2b
* (0.00252588685 - 0.00001491860 * Sin2_b + 0.00000011904 // В ГОСТЕ ОШИБКА!!!!
* Sin2_b * Sin2_b);
// System.out.println("Bo="+Bo);
double n = (int) (Y * Math.pow(10, -6));
// System.out.println("n="+n);
double Sin_Bo = Math.sin(Bo);
double Sin2_Bo = Math.pow(Sin_Bo, 2);
double Sin4_Bo = Math.pow(Sin_Bo, 4);
double Sin6_Bo = Math.pow(Sin_Bo, 6);
double Cos_Bo = Math.cos(Bo);
double Sin2Bo = Math.sin(Bo * 2);
double Zo = (Y - (10 * n + 5) * 100000) / (6378245 * Cos_Bo);
//System.out.println("Zo="+Zo);
double Zo2 = Math.pow(Zo, 2);

double 
dB=-Zo2*Sin2Bo*(0.251684631-0.003369263*Sin2_Bo+0.000011276*Sin4_Bo-
Zo2*(0.10500614-0.04559916* Sin2_Bo+0.00228901* Sin4_Bo-0.00002987* Sin6_Bo-
Zo2*(0.042858-0.025318* Sin2_Bo+0.014346* Sin4_Bo-0.001264* Sin6_Bo-
Zo2*(0.01672-0.00630* Sin2_Bo+0.01186* Sin4_Bo-0.00328* Sin6_Bo))));
double l= Zo*(1-0.0033467108* Sin2_Bo-0.0000056002*Sin4_Bo-0.0000000187*Sin6_Bo-
Zo2*(0.16778975+0.16273586* Sin2_Bo-0.00052490* Sin4_Bo-0.00000846* Sin6_Bo-
Zo2*(0.0420025+0.1487407* Sin2_Bo+0.0059420* Sin4_Bo-0.0000150* Sin6_Bo-
Zo2*(0.01225+0.09477* Sin2_Bo+0.03282* Sin4_Bo-0.00034* Sin6_Bo-
Zo2*(0.0038+0.0524* Sin2_Bo+0.0482* Sin4_Bo+0.0032* Sin6_Bo)))));

//System.out.println("dB="+dB);
//System.out.println("l="+l);
// ГОСТ 32453-2013
Point[0]=Math.toDegrees(Bo+dB);
Point[1]=Math.toDegrees((6*(n-0.5)/57.29577951)+l);
//System.out.println("B="+Point[0]);
//System.out.println("L="+Point[1]);
return Point;
}
 
// возрашает строку вида GG°MM′SS″ 
public static String returnGMS (double Direct) {
double Gg = Direct;
int G = (int) Gg;
double Mg = (Gg -G) * 60;
int M = (int) Mg;
double S = (Mg - M) * 60;
String Out =((int) G+"°"+(int)M +"′"+S+"″" );
return Out;
}
 
 
}
 
 

 

Доброго времени суток.
И очередной вопрос.
Пытаюсь кратко сформулировать. Есть группа точек в пространственной системе координат A и системе координат Б.
Каким образом вычислить 7 параметров. Нужно грамотное описание этого алгоритма. "Из серии сделай сам"

Оффтоп

Т.к. я пошел по очень кривому пути.
1. Нашел центр масс для системы АиБ (определил dX, dY,dZ) 2. Из центра масс построил вектора на все точки. (Определил wZ (wX,wY пока не прописал), и на этом этапе понял, что я иду по очень долгому пути. Алгоритм (вернее его реализация) получился очень долгим сложным. А где то в просторах инета попадалась статься которую я читал, но не сохранил в закладках с красивым решением. А сейчас найти не могу.

 

никак, слишком много неизвестных: датум, проекция - каждый из которых описывается набором неизвестных переменных...
А вот если предпологается, что датум и проекция одинаковые, а отличаются только параметры проекции - то здесь уже можно, что то придумать, только для каждого типа проекции алгоритмы будут немного отличатся...