Добро пожаловать!

Войдите или зарегистрируйтесь сейчас!

Войти

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

Тема в разделе "Исходные данные", создана пользователем RoBo, 22 апр 2015.

  1. RoBo

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

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

    Форумчанин

    Регистрация:
    5 янв 2008
    Сообщения:
    4.172
    Симпатии:
    11.938
    Адрес:
    Златоглавая и Белокаменная
    ???
    Нет.
    Да.
    То, что вы привели:
    это и есть параметры проекции Гаусса-Крюгера для
    Формулы можно и из ГОСТа, можно из Coordinate Conversions and Transformations including Formulas Revised - April 2015.
    Самые простые (компактные) и достаточно точные формулы даны в Transverse Mercator: Bowring series
    Посмотрите ссылки в сообщении
    http://geodesist.ru/forum/threads/Преобразование-координат-из-одной-плоской-системы-в-другую.42027/#post-479425
    А это-то вам зачем?
    Остаётся только «Главный вопрос жизни, вселенной и всего такого».
    Зачем вам это надо? Для какой такой задачи понадобился этот переход?
    А также, откуда у вас точные параметры связи из
    ? Как и кто это определял?
    Удивительное рядом.
     
    #23
    Stepan_S нравится это.
  4. Stepan_S

    Регистрация:
    25 авг 2010
    Сообщения:
    22
    Симпатии:
    5
    Адрес:
    СПб
    Это было банальное утверждение, что это все параметры перехода проекции Гаусса-Крюгера на геоиде с параметрами Красовского.


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

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

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

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

    Форумчанин

    Регистрация:
    5 янв 2008
    Сообщения:
    4.172
    Симпатии:
    11.938
    Адрес:
    Златоглавая и Белокаменная
    Оставим в стороне рениксы типа "геоид Красовского" и
    Ближе к теме.
    Во втором предложении непонятно зачем говорить о долготе, от которой длина дуги ну никак не зависит. А первое предложение справедливо лишь частично. Оно справедливо только для точек на осевом меридиане. Длину дуги меридиана S от экватора до параллели с геодезической широтой B можно записать как

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

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


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

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

    Вообще-то, умножение матриц обладает свойством ассоциативности, поэтому перемножив один раз матрицы (и запомнив/сохранив результат), можно напрямую переходить от WGS 84 к СК-42.
    Тогда посмотрите C++ and Java Code for Geodesic and Meridian Arc Computations -- by Klaus Hehl.
     
    #25
  6. Родичкин

    Форумчанин

    Регистрация:
    7 июл 2010
    Сообщения:
    2.074
    Симпатии:
    2.134
    Оффтоп

    Соседняя тема по непонятным причинам называцца "Помогите перевести координаты участка из МСК в 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 метров?
     
    #26
  7. trir

    Форумчанин

    Регистрация:
    25 ноя 2014
    Сообщения:
    3.253
    Симпатии:
    931
    Адрес:
    gnomtrir@mail.ru
     
    #27
  8. -=13=-

    Форумчанин

    Регистрация:
    26 июн 2013
    Сообщения:
    2.254
    Симпатии:
    3.320
    Адрес:
    Окраины Нерезиновска на немцеопасном направлении
    Оффтоп

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

    Форумчанин

    Регистрация:
    5 янв 2008
    Сообщения:
    4.172
    Симпатии:
    11.938
    Адрес:
    Златоглавая и Белокаменная
    Наличие развитой индустрии моделирования нисколько не мешает любителям клепать модельки чего-либо своими руками. Им это просто интересно.
    Для тех, у кого вместо карты пачка Беломора, этот вопрос не имеет смысла.::biggrin24.gif::
     
    #29
    Stepan_S нравится это.
  10. Родичкин

    Форумчанин

    Регистрация:
    7 июл 2010
    Сообщения:
    2.074
    Симпатии:
    2.134
    Оффтоп

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

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

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

    Форумчанин

    Регистрация:
    26 июн 2013
    Сообщения:
    2.254
    Симпатии:
    3.320
    Адрес:
    Окраины Нерезиновска на немцеопасном направлении
    Оффтоп

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

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

    Форумчанин

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

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

    X-Y-H Администратор
    Команда форума Форумчанин

    Регистрация:
    18 май 2007
    Сообщения:
    21.825
    Симпатии:
    7.104
    Адрес:
    Россия
    Оффтоп
    Родичкин, честно, создайте тему про сетки и СК и там все пожалуйста в одном месте. Человек пришел за конкретным, а вы демагогию опять разводите.
     
    #33
  14. Родичкин

    Форумчанин

    Регистрация:
    7 июл 2010
    Сообщения:
    2.074
    Симпатии:
    2.134
    Оффтоп

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

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

    Регистрация:
    25 авг 2010
    Сообщения:
    22
    Симпатии:
    5
    Адрес:
    СПб
    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 (раскрыть)


    // 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 ;}
     
    #35
    -=13=- нравится это.
  16. Stepan_S

    Регистрация:
    25 авг 2010
    Сообщения:
    22
    Симпатии:
    5
    Адрес:
    СПб
    Так и делаем.
     
    #36
  17. Yuri V.

    Yuri V. Только чтение
    Форумчанин

    Регистрация:
    31 мар 2009
    Сообщения:
    2.335
    Симпатии:
    2.085
    Адрес:
    Ивантеевка, РФ
    Родичкин, съездите в отпуск, например в Крым. Тут Роскартография длинные щупальца распустила. Полным ходом идёт обследование пунктов, спутниковые наблюдения на них. Программа предусм. перенаблюдать все пункты 1-2 классов. И это будет сделано. Какая для вас к чёрту разница сколько будет стоить?! Что вы заладили!
    -=13=-, всё будет хо-ро-шо! Работа идёт! ::smile24.gif::
     
    #37
    -=13=- нравится это.
  18. Stepan_S

    Регистрация:
    25 авг 2010
    Сообщения:
    22
    Симпатии:
    5
    Адрес:
    СПб
    Коли тема так называться. Привожу пример рабочих, проверенных методов преобразования координат на 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;
    }
     
     
    }
     
     
    
     
    #38
    Последнее редактирование: 11 июн 2015
    Фёдоров Роман и X-Y-H нравится это.
  19. Stepan_S

    Регистрация:
    25 авг 2010
    Сообщения:
    22
    Симпатии:
    5
    Адрес:
    СПб
    Доброго времени суток.
    И очередной вопрос.
    Пытаюсь кратко сформулировать. Есть группа точек в пространственной системе координат A и системе координат Б.
    Каким образом вычислить 7 параметров. Нужно грамотное описание этого алгоритма. "Из серии сделай сам"

    Оффтоп

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

    Форумчанин

    Регистрация:
    25 ноя 2014
    Сообщения:
    3.253
    Симпатии:
    931
    Адрес:
    gnomtrir@mail.ru
    никак, слишком много неизвестных: датум, проекция - каждый из которых описывается набором неизвестных переменных...
    А вот если предпологается, что датум и проекция одинаковые, а отличаются только параметры проекции - то здесь уже можно, что то придумать, только для каждого типа проекции алгоритмы будут немного отличатся...
     
    #40

Поделиться этой страницей

  1. Этот сайт использует файлы cookie. Продолжая пользоваться данным сайтом, Вы соглашаетесь на использование нами Ваших файлов cookie.
    Скрыть объявление