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

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

Войти

Как перевести широту и долготу от GPS в прямоугольные координаты?

Тема в разделе "Исходные данные", создана пользователем victor64, 8 ноя 2011.

  1. Ладиков Павел

    Форумчанин

    Регистрация:
    10 май 2011
    Сообщения:
    276
    Симпатии:
    80
    Адрес:
    Саранск
    Это смысл задачи? Тяжеловата...
    Я попробую переформулировать. Итак:
    Что Вы собираетесь сделать с полученными треками? До или после преобразования, в том числе? От ответа и будет зависеть метод получения координат.
    Если любоваться на красивый клубочек линий - это одно. Если использовать в дальнейшем как картооснову или ее элемент - это совсем другое. Там без понятия "система координат" обойтись не получится.

    Теперь о точке отсчета. Она может быть каждый раз новой для каждого мероприятия. В этом случае у Вас каждый раз будут создаваться новые системы координат и Вы не сможете так совместить разные треки - чтобы они между собой соответствовали правильным расстояниям.

    Если же Вы выбираете одну статическую точку отсчета, которая у Вас будет всегда... можно конечно ее придумать - тогда в приведенных формулах пересчета координат на основе госта - смещения в этой точке должны быть нулевыми, базовый меридиан равен нулю. Но если слегка помучиться и найти точку отсчета для принятой в вашей местности системы координат - тогда Вы мсожете использовать накапливаемые данные для совмещения с кадастровой информацией. Вдруг понадобится?
    Проще всего за такую точку принять реально существующий на местности геодезический или межевой опорный пункт или знак. Координаты межевых знаков и описания как их найти - можно получить, если заказать выписку кадастрового плата территории вашего кадастрового квартала. Они по идее должны их теперь указывать в конце выписки.
     
    #41
  2. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    да да да, именно так, больше ничего, просто я вижу в реальном времени на каком расстоянии я лечу от соседнего трека, потм пролетели и на следующий день забыли где летали ::biggrin24.gif::
     
    #42
  3. ЮС

    Форумчанин

    Регистрация:
    28 фев 2010
    Сообщения:
    4.567
    Симпатии:
    5.115
    Прямая равноугольная цилиндрическая проекция Меркатора и
    Равноугольная поперечно-цилиндрическая проекция Меркатора
    это не одно и то же.

    Почувствуйте разницу.

    На картинке два фрагмента из одного справочника,
    а по ссылке популярное изложение темы о различных проекциях:

    http://www.giscraft.ru/info/cs/cs.shtml
     

    Вложения:

    #43
  4. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    Теперь о точке отсчета. Она может быть каждый раз новой для каждого мероприятия. В этом случае у Вас каждый раз будут создаваться новые системы координат и Вы не сможете так совместить разные треки - чтобы они между собой соответствовали правильным расстояниям.
    Треки должны совмещаться только для одного поля
     
    #44
  5. Ладиков Павел

    Форумчанин

    Регистрация:
    10 май 2011
    Сообщения:
    276
    Симпатии:
    80
    Адрес:
    Саранск
    Один алгоритм? Нуда-нуда... вся музыка похожа - 12 одних и тех же нот.
    Все формулы преобразования координат похожи - все те же значи сложения... вычитания... деления и умножения...
    Вы не скажете, в таком случае - как это по одинаковому алгоритму получаются такие разные картинки? И как изменить "масштаб по осевому меридиану" чтобы картинки совпали? Растянуть по вертикали? Или сжать по вертикали?

    ЗЫ: Извините, это уже не смешно...

    [​IMG]

    [​IMG]
     
    #45
  6. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    карты полей есть на бумаге и др. носителях, например в обычном приемнике и существует куча разных программ для расчета площади, все это эсть, нужна решить только одну задачу - нарисовать "клубок" треков
     
    #46
  7. Ладиков Павел

    Форумчанин

    Регистрация:
    10 май 2011
    Сообщения:
    276
    Симпатии:
    80
    Адрес:
    Саранск
    Ага! У них даже название отличается! Кстати... я тоже умею пользоваться гуглом.
    (Добавление)
    Ну вот... уже теплее... через недельку можно будет переходить к делу 8)))
    Теперь вопрос:
    Нарисовать клубок треков... на чистом листе? Тогда зачем в метры, рисуйте в градусах!
    Или таки на карте? В этом случае все только начинается. И начинается не с хвоста и копыт, а носа: нужно абсолютно точно узнать - в какой системе координат и в какой проекции изготовлена Ваша карта.
    И работать либо в этой системе, лмбо преобразовать (трансформировать) карту в ту систему - в которой Вы планируете работать. Можно даже в меркатора - тогда не придется заниматься конвертацией координат треков. Заодно сможете накладывать космические снимки из гугла.
     
    #47
  8. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    вообще моё место работы в радиусе 30 км от точки N51гр 30,24мин E43гр 36,75мин. стоит ли вычислять dN и dE для данного радиуса действий? если эта разница до сотых долей, то меня это в полне устраивает. Иногда приходится работать на удалении 100-300км от данной точки.
    А можно алгоритм ::rolleyes24.gif:: ? или это интелектуальная собственность?
     
    #48
  9. ВасЯ

    Форумчанин

    Регистрация:
    27 фев 2008
    Сообщения:
    258
    Симпатии:
    48
    Адрес:
    Рязанская область
    Простите, пожалуйста, но я повторюсь. Ответ был дан выше:
    Если хотите просто видеть линии треков и при этом иметь возможность оценить в метрах (и даже долях метров) расстояние между треками (например, чтобы четко выдержать полосу полива или разброса удобрений), то для начала Вам будет достаточно этого.

    Например, берете первую точку трека и принимаете за начало координат. B - широта этой точки. dB и dL - приращения координат (широты и долготы) от начальной точки до текущей. dN и dE - то, что нужно в итоге нарисовать на экране ("север" и "восток").
    Какие-либо величины из формул надо прокомментировать (полуось эллипсоида, эксцентриситет, сжатие)?

    А все точки трека зачем рисовать? Для вывода на экран сформируйте упрощенную линию, а трек просто храните в памяти.
     
    #49
  10. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    Ды не нужны ни какие карты ::biggrin24.gif:: а в метрах, потому, что удобнее рисовать, например 1м это 4 пикселя на экране. все координаты пишутся в массив и прорисовываются 5 раз в секунду. если мне нужно отступить 30м от соседнего трека я рисую трек кистью в 120 пикселей (30м * 4пикселя) и все вижу на экране
    (Добавление)
    ДА спасибо, я этим уже занимаюсь - если точки лежат на одной прямой или на неком небольшом расстоянии от неё (1-2м), то они не беруться в расчет и не рисуются
    (Добавление)
    Если можно то прокомментируйте
    f(2-f) - что за функция?
    a - ?
    и это ведь часть программы ?

    dN = a*(1-e^2)/W^3*dB
    dE = a*cosB/W*dL
    W = sqrt(1-(e*sinB)^2); e^2 = f(2-f); f~1/300
    Еще бы остальное, или я что то не правильно понял?

     
    #50
  11. ВасЯ

    Форумчанин

    Регистрация:
    27 фев 2008
    Сообщения:
    258
    Симпатии:
    48
    Адрес:
    Рязанская область
    Не правильно понял. stout описал алгоритм действий, если нужно получить координаты точек в горизонтальной плоскости, касательной к данной точке эллипсоида.

    А вот формулки, процитированные мной - "длинна" 1 угловой секунды широты и долготы соответственно.

    В формулах f - сжатие эллипсоида, а - его большая полуось, e - его эксцентриситет, который можно вывести из сжатия по приведенной формуле.
     
    #51
  12. ZUCKtm

    Форумчанин

    Регистрация:
    21 май 2007
    Сообщения:
    1.832
    Симпатии:
    142
    Адрес:
    г. Мытищи Московской области
    2 Ладиков Павел
    Ну и что за картинки вы выложили? На первой - равноугольная цилиндрическая проекция Меркатора, на второй - коническая проекция. Мы о второй говорили? Она отношение к UTM имеет?
     
    #52
  13. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    т.е. достаточно подставить dB и dL из приёмника (широту и долготу) и вуаля, мы вычислим длину секунды широты и долготы? и всё ::blink.gif:: ?
    (Добавление)
    все, вроде разобрался, больше ни каких вычислений не надо, кроме вычисления dB и dL&?
    (Добавление)
    а эти параметры для какого элепсоида красовского или wgs84?
     
    #53
  14. ZUCKtm

    Форумчанин

    Регистрация:
    21 май 2007
    Сообщения:
    1.832
    Симпатии:
    142
    Адрес:
    г. Мытищи Московской области
    Навигатор выдает результат в WGS, его параметры и подставляйте.
     
    #54
  15. stout

    Форумчанин

    Регистрация:
    5 янв 2008
    Сообщения:
    4.172
    Симпатии:
    11.939
    Адрес:
    Златоглавая и Белокаменная
    Боротся с воинствующим невежеством не вижу смысла. На каждую твою поправку, уточнение сваливается пара новых ляпов.
    Это тот самый редкий случай, когда Павел прав. (Речь идёт только об этой теме, его мессаджи о GPS в подавляющем случае верны)
    Как справедливо отмечено в Wiki "The UTM system is not a single map projection. The system instead employs a series of sixty zones, each of which is based on a specifically defined secant transverse Mercator projection."

    Если вы смотрели DMA Technical Manual 8358.2 (полное название:"THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)"), то должны были заметить, что речь идет о UTM grid. И дано определение
    Т.е. речь здесь идет о системе (картографических) координат, но не только о проекции.
    Когда говорят только о проекции, то говорят просто о transverse Mercator projection
    Есть подозрение, что тот, кто первым употребил в русском языке термин "проекция UTM" имел ввиду родительный падеж, а не аббревиатуру UTM, как название проекции. Но здесь могу и ошибаться.
    (Приблизительно по аналогичной причине, я предпочитаю переводить IERS - International Earth Rotation and Reference Systems Service, не дословно, как Международная служба вращения Земли и референцных систем, а как "Международная служба референцных систем и вращения Земли". Иначе можно предположить, что референцные системы должны вращаться)

    А по поводу похожести там есть одно интересное место.
    В разделе 2-3 SPECIFICATIONS OF THE UTM. сказано:
    Прошу обратить внимание - по мнению авторов документа, поперечная прекция Меркатора принадлежит к проекциям типа Гаусса-Крюгера.

    Я уже неоднократно писал об этом. Но не грех и повториться.
    По моему, причиной ваших заблуждений является целый сонм неверных иллюстраций, коими пестрят не только популярные издания, но и некоторые монографии и более или менее солидные учебники. В интернете с этим просто беда.
    Достаточно посмотреть рисуночек из "ArcGIS® 9 Understanding Map Projections" с лампочкой в центре Земли и обёртывающим цилиндром.
    Очень редко можно найти правильные толкования, например такое как это:
    "The Mercator projection is often illustrated as a projection from the center of a globe onto a cylinder wrapped around the globe. However, that cylindrical projection is not the Mercator projection. If you see them side by side, you will see they dont look the same. The cylindrical projection does not have the property that Mercator was seeking; the rhumb lines are not straight lines. They do have the general form in common: meridians of longitude and parallels of latitude are straight, parallel lines. Also, since the meridians dont get closer together towards the poles (as they do on a globe), that means the scale must increase as the latitude increases." (http://babylon.acad.cai.cam.ac.uk/stude ... pirate.htm)
    Идея с цилиндром, на который происходит проектирование, это всего лишь иллюстрация, помогающая понять на интуитивном уровне, что происходит на самом деле. Вспомните о том огромном числе проекций, с помощью которых пытаются передать трехмерность на плоскости и сколько этих проекций впервые придумали художники.
    Лучше всего ао этому поводу сказал Морозов в своей книге "Курс сфероидической геодезии"
    Да без проблем.
    Текст прожки (раскрыть)
    Код:
    unit MainUnit;
    
    interface
    
    uses
      Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
      Dialogs, StdCtrls, JvExStdCtrls, JvEdit, JvValidateEdit
      , UnitTMercator, Buttons
      ;
    
    type
      TForm1 = class(TForm)
        memo: TMemo;
        editB: TJvValidateEdit;
        editDX: TJvValidateEdit;
        editL: TJvValidateEdit;
        editDY: TJvValidateEdit;
        txt1: TStaticText;
        txt2: TStaticText;
        txt3: TStaticText;
        txt4: TStaticText;
        txt5: TStaticText;
        editH: TJvValidateEdit;
        btnCalc: TBitBtn;
        procedure FormCreate(Sender: TObject);
        procedure btnCalcClick(Sender: TObject);
      private
        { Private declarations }
      public
        a, recF, f, e: Double;
        Bo, Lo, H, dX, dY: Double;
        B, L, X, Y, Z: Double;
        Xp, Yp, Zp, Xo, Yo, Zo: Double;
        M: array[1..3, 1..3] of Double;
        { Public declarations }
        TM: TransverseMercator;
      end;
    
    var
      Form1: TForm1;
    
    implementation
    
    {$R *.dfm}
    
    procedure BLH2XYZ(a, e, B, L, H: Double; var X, Y, Z: Double);
    var
      W, N: Double;
    begin
      W := Sqrt(1 - sqr(e * Sin(B)));
      N := a / W;
      X := (N + H) * Cos(B) * Cos(L);
      Y := (N + H) * Cos(B) * Sin(L);
      Z := (N * (1 - e * e) + H) * Sin(B);
    end;
    
    procedure TForm1.FormCreate(Sender: TObject);
    begin
      a := 6378245.0;
      recF := 298.3;
      f := 1.0 / recF;
      e := Sqrt(f * (2 - f));
      TM := TransverseMercator.Create(a, recF, 0, 0, 7, 1.0);
    
    end;
    
    procedure TForm1.btnCalcClick(Sender: TObject);
    var
      W : Double;
    begin
      Bo := editB.AsFloat / 180 * Pi;
      Lo := editL.AsFloat / 180 * Pi;
      H := editH.AsFloat;
      dX := editDX.AsFloat;
      dY := editDY.AsFloat;
      memo.Lines.Add('Bo = ' + FloatToStr(Bo * 180 / Pi));
      memo.Lines.Add('Lo = ' + FloatToStr(Lo * 180 / Pi));
    
      W := Sqrt(1 - sqr(e * Sin(Bo)));
      memo.Lines.Add('dN = ' + FloatToStr(a*(1-sqr(e))/(W*W*W)/206264.806247));
      memo.Lines.Add('dE = ' + FloatToStr(a*cos(Bo)/W/206264.806247));
    
      memo.Lines.Add('H  = ' + FloatToStr(H));
      memo.Lines.Add('dX = ' + FloatToStr(dX));
      memo.Lines.Add('dY = ' + FloatToStr(dY));
      TM.BL2XY(Bo, 0, X, Y);
      BLH2XYZ(a, e, Bo, Lo, H, Xo, Yo, Zo);
    
      M[1][1] := -Sin(Lo);
      M[1][2] := Cos(Lo);
      M[1][3] := 0;
      M[2][1] := -Sin(Bo) * Cos(Lo);
      M[2][2] := -Sin(Bo) * Sin(Lo);
      M[2][3] := Cos(Bo);
      M[3][1] := Cos(Bo) * Cos(Lo);
      M[3][2] := Cos(Bo) * Sin(Lo);
      M[3][3] := Sin(Bo);
    
      X := X + dX;
      Y := Y + dY;
      TM.XY2BL(X, Y, B, L);
      memo.Lines.Add('B = ' + FloatToStr(B * 180 / Pi));
      memo.Lines.Add('Lo = ' + FloatToStr(L * 180 / Pi));
      BLH2XYZ(a, e, B, L + Lo, H, Xp, Yp, Zp);
      x := M[1][1] * (Xp - Xo) + M[1][2] * (Yp - Yo) + M[1][3] * (Zp - Zo);
      y := M[2][1] * (Xp - Xo) + M[2][2] * (Yp - Yo) + M[2][3] * (Zp - Zo);
      z := M[3][1] * (Xp - Xo) + M[3][2] * (Yp - Yo) + M[3][3] * (Zp - Zo);
    
      memo.Lines.Add('x = ' + FloatToStr(x));
      memo.Lines.Add('y = ' + FloatToStr(y));
      memo.Lines.Add('z = ' + FloatToStr(z));
    
    end;
    
    end.
    
    
    
    unit UnitTMercator;
    
    interface
    uses
      Math
    , SysUtils
    ;
    
    const
      MinDegree = 4;
      MaxDegree = 7;
    type
      TransverseMercator = class
      private
        FNorthing: Double; // False Northing(Southing)
        FEasting: Double; // False Easting
        n: Double; // n - Third flattening
        Eccentricity: Double;
        Gm: Double;
        G: array[0..MaxDegree] of Double;
        F: array[0..MaxDegree] of Double;
        H: array[0..MaxDegree] of Double;
        CurrentDegree: Integer;
      protected
    
      public
        constructor Create(const SemiMajorAxis: Double; // SemiMajaorAxis
          const invFlattening: Double; // reciprocal of the flattening (1/Flattening ~300)
          const FalseNorthing: Double; // False Northing(Southing)
          const FalseEasting: Double;  // False Easting
          const degree: Integer;
          const scale: Double          // scale
          );
        procedure BL2XY(const Lat, dLon: Double; var X, Y: Double);
        procedure XY2BL(const X, Y: Double; var Lat, dLon: Double);
      end;
    
    //
    // =============================================================================
    //
    
    implementation
    
    const pseG: array[0..27] of Double = (
        +1.0/2.0, -2.0/3.0, +5.0/16.0, +41.0/180.0, -127.0/288.0, +7891.0/37800.0, +72161.0/387072.0,
        +13.0/48.0, -3.0/5.0, +557.0/1440.0, +281.0/630.0, -1983433.0/1935360.0, +13769.0/28800.0,
        +61/240.0, -103.0/140.0, +15061.0/26880.0, +167603.0/181440.0, -67102379.0/29030400.0,
        +49561.0/161280.0, -179.0/168.0, +6601661.0/7257600.0, +97445.0/49896.0,
        +34729.0/80640.0, -3418889.0/1995840.0, 14644087.0/9123840.0,
        +212378941.0/319334400.0, -30705481.0/10378368.0,
        +1522256789.0/1383782400.0);
    
    const pseH: array[0..27] of Double = (
        +2.0, -2.0/3.0, -2.0, +116.0/45.0, +26.0/45.0, -2854.0/675.0, +16822.0/4725.0,
        +7.0/3.0, -8.0/5.0, -227.0/45.0, +2704.0/315.0, +2323.0/945.0, -31256.0/1575.0,
        +56.0/15.0, -136.0/35.0, -1262.0/105.0, +73814.0/2835.0, +98738.0/14175.0,
        +4279.0/630.0, -332.0/35.0, -399572.0/14175.0, +11763988.0/155925.0,
        +4174.0/315.0, -144828.0/6237.0, -2046082.0/31185.0,
        +601676.0/22275.0, -115444544.0/2027025.0,
        +38341552.0/675675.0);
    
    const pseF: array[0..27] of Double = (
        -1.0/2.0, +2.0/3.0, -37.0/96.0, +1.0/360.0, +81.0/512.0, -96199.0/604800.0, +5406467.0/38707200.0,
        -1.0/48.0, -1.0/15.0, +437.0/1440.0, -46.0/105.0, +1118711.0/3870720.0, -51841.0/1209600.0,
        -17.0/480.0, +37.0/840.0, +209.0/4480.0, -5569.0/90720.0, -9261899.0/58060800.0,
        -4397.0/161280.0, +11.0/504.0, +830251.0/7257600.0, -466511.0/294800.0,
        -4583.0/161280.0, +108847.0/3991680.0, +8005831.0/63866880.0,
        -20648693.0/63866880.0, +16363163.0/518918400,
        -219941297.0/5535129600.0);
    
    //
    // =============================================================================
    //
    
    constructor TransverseMercator.Create(const SemiMajorAxis: Double; // SemiMajaorAxis
      const invFlattening: Double; // reciprocal of the flattening (1/Flattening ~300)
      const FalseNorthing: Double; // False Northing(Southing)
      const FalseEasting: Double;  // False Easting
      const degree: Integer;
      const scale: Double          // scale
      );
    var
      sF, sG, sH   : Double;
      np, n2       : Double;
      i, j, k, base: Integer;
    begin
      FNorthing := FalseNorthing;
      FEasting  := FalseEasting;
      CurrentDegree := min(MaxDegree, max(MinDegree, degree));
    
      n := 1/(2*invFlattening - 1); //  ThirdFlattening
      Eccentricity := 2*sqrt(n)/(1 + n);
      np := 1;
      base := 0;
    
      for k := 0 to MaxDegree - 1 do
      begin
        j := base + MaxDegree - 1 - k;
        sF := n*pseF[j];
        sG := n*pseG[j];
        sH := n*pseH[j];
        for i := MaxDegree - 2 - k downto 0 do
        begin
          j := base + i;
          sF := sF + pseF[j];
          sG := sG + pseG[j];
          sH := sH + pseH[j];
          sF := sF*n;
          sG := sG*n;
          sH := sH*n;
        end;
        F[k] := sF*np;
        G[k] := sG*np;
        H[k] := sH*np;
        np := n*np;
        base := base + (MaxDegree - k);
      end;
      n2 := n*n;
      Gm := scale*SemiMajorAxis/(n + 1) *
        (1. + n2*(1.0/4.0 + n2*(1.0/64.0 + n2*(1.0/256.0 + n2*25.0/16384.0))));
    end;
    
    //
    // =============================================================================
    //
    
    procedure TransverseMercator.BL2XY(const Lat, dLon: Double; var X, Y: Double);
    var
      sinB, esinB, zz: Double;
      tanBeta1, tanhBeta2, Beta1, Beta2: Double;
      tg2, th2: Double;
      sin2Beta1, cos2Beta1: Double;
      sinh2Beta2, cosh2Beta2: Double;
      ReSigma, ImSigma: Double;
      ReSn, ImSn, Resnp1, ImSnp1, Resnp2, ImSnp2: Double;
      k: Integer;
    begin
      sinB := sin(Lat);
      esinB := Eccentricity*sinB;
      zz := sqrt((1 + sinB)/(1 - sinB)*Power((1 - esinB)/(1 + esinB), Eccentricity));
    //   zz := tan(M_PI/4+Lat/2)*pow((1-esinB)/(1+esinB), Eccentricity/2.0);
      tanBeta1 := (zz - 1/zz)/(2*cos(dLon)); // tan(Beta1) := sh(H)/cos(l), H:=ln(z)
      tanhBeta2 := 2*sin(dLon)/(zz + 1/zz); // th(Beta2) := sin(l)/ch(H)
      Beta1 := ArcTan(tanBeta1);
      Beta2 := Ln((1 + tanhBeta2)/(1 - tanhBeta2))/2;
      tg2 := tanBeta1*tanBeta1;
      th2 := tanhBeta2*tanhBeta2;
      sin2Beta1 := 2*tanBeta1/(1 + tg2);
      cos2Beta1 := (1 - tg2)/(1 + tg2);
      sinh2Beta2 := 2*tanhBeta2/(1 - th2);
      cosh2Beta2 := (1 + th2)/(1 - th2);
      ReSigma := 2*cos2Beta1*cosh2Beta2;
      ImSigma := -2*sin2Beta1*sinh2Beta2;
    //  Clenshaw algorithm initialization
      ReSn := ReSigma*G[CurrentDegree - 1] + G[CurrentDegree - 2];
      ImSn := ImSigma*G[CurrentDegree - 1];
      ReSnp1 := G[CurrentDegree - 1];
      ImSnp1 := 0;
    //  loop of Clenshaw algorithm
      for k := CurrentDegree - 3 downto 0 do
      begin
        ReSnp2 := ReSnp1;
        ImSnp2 := ImSnp1;
        ReSnp1 := ReSn;
        ImSnp1 := ImSn;
        ReSn := -ReSnp2 + ReSigma*ReSnp1 - ImSigma*ImSnp1 + G[k];
        ImSn := -ImSnp2 + ReSigma*ImSnp1 + ImSigma*ReSnp1;
      end;
      ReSigma := sin2Beta1*cosh2Beta2;
      ImSigma := cos2Beta1*sinh2Beta2;
    //  end of Clenshaw algorithm
      X := ReSigma*ReSn - ImSigma*ImSn;
      Y := ReSigma*ImSn + ImSigma*ReSn;
      X := Gm*(Beta1 + X) + FNorthing;
      Y := Gm*(Beta2 + Y) + FEasting;
       //
    end;
    //
    // =============================================================================
    //
    
    procedure TransverseMercator.XY2BL(const X, Y: Double; var Lat, dLon: Double);
    var
      U, V           : Double;
      cos2U, cosh2V  : Double;
      sin2U, sinh2V  : Double;
      ReBeta, ImBeta : Double;
      ReSn, ImSn     : Double;
      ReSnp1, ImSnp1 : Double;
      ReSnp2, ImSnp2 : Double;
      Beta1, Beta2   : Double;
      eb, sinB, sin2B: Double;
      k: Integer;  
    begin
      U := (X - FNorthing)/Gm;
      V := (Y - FEasting)/Gm;
      cos2U := Cos(2*U); cosh2V := Cosh(2*V);
      sin2U := Sin(2*U); sinh2V := Sinh(2*V);
      ReBeta := 2*cos2U*cosh2V;
      ImBeta := -2*sin2U*sinh2V;
    
      ReSn := ReBeta*F[CurrentDegree - 1] + F[CurrentDegree - 2];
      ImSn := ImBeta*F[CurrentDegree - 1];
      ReSnp1 := F[CurrentDegree - 1];
      ImSnp1 := 0;
      for k := CurrentDegree - 3 downto 0 do
      begin
        ReSnp2 := ReSnp1;
        ImSnp2 := ImSnp1;
        ReSnp1 := ReSn;
        ImSnp1 := ImSn;
        ReSn := -ReSnp2 + ReBeta*ReSnp1 - ImBeta*ImSnp1 + F[k];
        ImSn := -ImSnp2 + ReBeta*ImSnp1 + ImBeta*ReSnp1;
      end;
      ReBeta := sin2U*cosh2V;
      ImBeta := cos2U*sinh2V;
      Beta1 := ReBeta*ReSn - ImBeta*ImSn + U;
      Beta2 := ReBeta*ImSn + ImBeta*ReSn + V;
    
      eb := Exp(Beta2);
      sinB := 2*Sin(Beta1)/(eb + 1/eb);
      dLon := Arctan(0.5*(eb - 1/eb)/Cos(Beta1));
    
      ReBeta := 2 - 4*sqr(sinB);
      ReSn := ReBeta*H[MaxDegree - 1] + H[MaxDegree - 2];
      ReSnp1 := H[MaxDegree - 1];
      for k := MaxDegree - 3 downto 0 do
      begin
        ReSnp2 := ReSnp1;
        ReSnp1 := ReSn;
        ReSn := ReBeta*ReSnp1 - ReSnp2 + H[k];
      end;
      sin2B := 2*sinB*sqrt(1 - sqr(sinB));
      Lat := Arcsin(sinB) + sin2B*ReSn;
    end;
    //
    // =============================================================================
    //
    end.
    

    Чёрт! Тег code так и не работает!
     
    #55
    Geoshaman нравится это.
  16. Романыч

    Форумчанин

    Регистрация:
    5 ноя 2009
    Сообщения:
    1.892
    Симпатии:
    849
    Адрес:
    Россия, г. Иваново
    Попробуйте программа "Панорама". Я вообще не парюсь-создаю карту в СК-42, в Гармине получаю координаты в WGS-84. В программе переставляю галочку в создание объекта в координатах WGS, воожу данные с Гармина, переставляю галочку отображение координат в прямоугольных СК-42 и всё.
     
    #56
  17. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    stout попробовал Ваш алгоритм
    два вопроса:1. Правильно ли, что нетребуется ввод L0?
    2. Программа выдаёт 7значные числа до запятой?

    в excele перевожу координаты (гр,мин,сек) в градусы с десятичными долями по формуле:
    гр+мин/60+сек/3600
    затем в радианы
    db = PI / 180 * db
    dl = PI / 180 * dl
    b0 = PI / 180 * b0
    у меня правильный ход мыслей?
    ввожу данные: b0=52, db=52,00027778(смещение по b соответствует угловой секунде, по моим расчетам), dl=43


    Imports System.Math
    Public Class Form1
    Const a As Int32 = 6378245
    Const b As Int32 = 6356863
    Const f As Double = 0.00335233
    Const et As Double = 0.08181337

    Dim dn, de, db, dl, b0, w As Double
    Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
    db = TextBox2.Text
    dl = TextBox3.Text
    b0 = TextBox1.Text
    db = PI / 180 * db
    dl = PI / 180 * dl
    b0 = PI / 180 * b0

    w = Sqrt(1 - (et * Sin(b0)) ^ 2)
    dn = a * (1 - et ^ 2) / w ^ 3 * db
    de = a * Cos(b0) / w * dl
    Label4.Text = dn
    Label5.Text = de

    End Sub
    End Class
     
    #57
  18. chnav

    Форумчанин

    Регистрация:
    5 янв 2011
    Сообщения:
    1.003
    Симпатии:
    944
    Адрес:
    Москва
    victor64
    У вас программа что - на Excel+VBA работает ?

    Если нет, то пользуйтесь нормальными формулами и не надо ничего упрощать. Эта математика настолько проста для современного процессора в плане потраченных тактов, что выгода от оптимизации будет максимум 1 мкс.

    И еще. Обычно программисты при работе с проекциями предварительно расчитывают параметры, общие на весь проект (сжатие, полуоси, эксцентриситет - в C++ это делается на стадии инициализации объекта класса). В функциях перевода подставляются готовые значения. Например величины b0, w в вашем примере относится к "постоянным" величинам.

    PS: кстати ответ на ваш вопрос содержится в вашем коде
    dn = a * (1 - et ^ 2) / w ^ 3 * db
    de = a * Cos(b0) / w * dl

    т.е.
    - в одном [радиане] по широте содержится a*(1-et^2)/w^3 [метров]
    - в одном [радиане] по долготе содержится a*Cos(b0)/w [метров]
    w рассчитывается один раз на центр площадки (или произвольный угол, разницы вы не заметите).
     
    #58
  19. victor64

    Регистрация:
    6 май 2012
    Сообщения:
    0
    Симпатии:
    1
    нет конечно, и навигатор секунды не выдаёт, просто опрабирование данной формулы. Все оптимизации будут потом, это не долго зделать. Процедура не работает - вот в чем проблема и мне нужен совет непосредственно по расчетам, а не по программированию.
    (Добавление)
    поэтому и получается 7 значное число?
     
    #59
  20. chnav

    Форумчанин

    Регистрация:
    5 янв 2011
    Сообщения:
    1.003
    Симпатии:
    944
    Адрес:
    Москва
    db должно быть 0.0027778, а не 52.00027778
     
    #60

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

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