Формулы пересчета данных GPS-измерений из WGS-84 в СК-42

[Родичкин]

grig!

Не совсем уверен в том, что ГОСТ прошитый в Мапинфо (с датумом 9999, 3, 23.9, -141.3, -80.9, 0, -0.35, -0.86, -0.12, 0) дает точность 3-4 сантиметра. Более того, совсем не уверен.

Во первых, само Мапинфо без режима районирования в строчке файла проекций, дает ошибки координат 7.5 см.
Во вторых, сам ГОСТ дает по территории страны ошибки до 30 метров и более.
Проверено неоднократно.
Рекомендую посмотреть http://rutracker.org/forum/viewtopic.ph ... pmode=full

[Cocodrilo]

Не совсем понял, а что делает функция WGS84Alt ?

Ясно - недоделанная функция по определению высоты
Вместо нее надо дописать что-то типа:

Код: Выделить всё

Function WGS84_SK42_Alt(Bd, Ld, H) As Double
    WGS84_SK42_Alt = H - dH(Bd, Ld, H) 
End Function

Function SK42_WGS84_Alt(Bd, Ld, H) As Double
    SK42_WGS84_Alt = H + dH(Bd, Ld, H)
End Function

Function dH(Bd, Ld, H) As Double
    Dim B, L, N As Double
    B = Bd * Pi / 180
    L = Ld * Pi / 180
    N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5
    dH = -a / N * da + N * Sin(B) ^ 2 * de2 / 2 + (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B) - N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L)) + (a ^ 2 / N + H) * ms
End Function

Код не проверял, писал наугад

[Олег_7]

Попытался реализовать пересчет из WGS-84 в СК-95. Поправил параметры (там задействованы эллепсойд 42, повороты и масштабный коэффициент) получилась ошибка по широте 11 км. Обнулил масштабный коэффициент (примерно равный 0,5) - точность возросла до +-2см, но всё равно как-то боязно использовать в реальных расчетах.
С высоткой дела вообще плохо: погрешность 1.3 м. Может Cocodrilo или другие знатоки найдут ошибку. Буду благодарен хотя-бы за ссылки на первоисточники.

Мои мысли по теме: ввиду высокой (порядка 5-200м) деформации СК42 (некачественной закладки из-за отсутсвия спутникового оборудования в 50-70г прошлого века) точных параметров пересчета на территорию всей страны не существует. Можно рассчитать точные параметры для любого небольшого района (100-200 кв.км.), но расширяя этот район, мы начнем терять точность.Спасти положение могут матрицы типа Канадской NTv2, но насколько я знаю создавать их на всю страну никто не собирается.

[Донецков]

Кто-то создает и торгует NTv2 (где-то как бесплатный образец встречалось на Новосибирскую область), но все равно точность лучше 0.4 - 0.5 м вряд ли будет...

[tarabukinivan]

Добрый день!
Вроде сделал все правильно, но результат от WGS84_SK42_Lat(50;50;0) выходит совсем другой 49,9998280965696, а не 49.99980414. Формулы перепроверил 10 раз.

Код: Выделить всё

unit files;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Menus;
const
pi=3.14159265358979;
ro=206264.8062;
//эллипсоид красовского
ap=6378245;
alp=1/298.3; //0.00335232986925914
e2p=2*alp-alp*alp;//0.00669342162296594
//эллипсод WGS84
aw=6378137;
alw=1/298.257223563; //0.00335281066474748
e2w=2*alw-alw*alw;//0.0066943799901413

a=(ap+aw)/2;   //6378191
e2=(e2p+e2W)/2;     //0.0066939
da=aw-ap; //-108
de2=e2w-e2p;//9.58367175373768e-7
//линейные элементы трансформирования
dx=23.92;
dy=-141.27;
dz=-80.9;
//угловые элементы трансформирования
wx=0;
wy=0;
wz=0;
ms=0;//масштабный элемент трансформирования
type
  TForm1 = class(TForm)
    Button1: TButton;
    Edit1: TEdit;
    Button2: TButton;
    Memo1: TMemo;
    OpenDialog1: TOpenDialog;
    MainMenu1: TMainMenu;
    N1: TMenuItem;
    N2: TMenuItem;
    N4: TMenuItem;
    Label1: TLabel;
    N3: TMenuItem;
    Memo2: TMemo;
    Button3: TButton;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
    procedure N2Click(Sender: TObject);
    procedure N4Click(Sender: TObject);
    procedure N3Click(Sender: TObject);
    procedure Button3Click(Sender: TObject);


  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;filename:string;
implementation
{$R *.dfm}

Function dB(Bd,Ld,H:double):double;
var   B, L, M, N,skb1:double;
begin
    B:= Bd * Pi/180;
    L:= Ld * Pi/180;
    skb1:=n/a*e2*sin(b)*cos(b)*da+(sqr(n)/sqr(a)+1)*n*sin(b)*cos(b)*de2/2-(dx*cos(l)+dy*sin(l))*sin(b)+dz*cos(b);
    M:= a * (1 - e2) / sqrt((1 - e2 * Sin(b)*sin(b))*(1 - e2 * Sin(b)*sin(b))*(1 - e2 * Sin(b)*sin(b))) ;
    N:= a /sqrt(1 - e2 * sqr(Sin(B)));
    dB:= ro / (M + H) * skb1 - wx * Sin(L) * (1 + e2 * Cos(2 * B))+ wy * Cos(L) * (1 + e2 * Cos(2 * B))- ro * ms * e2 * Sin(B) * Cos(B);
End;
Function dL(Bd, Ld, H:double):double;
    Var B, L, N:double;
begin
    B:=Bd * Pi / 180;
    L:= Ld * Pi / 180;
    N:= a / sqrt(1 - e2 * Sin(B)*sin(b));
    dL:= ro / ((N + H) * Cos(B)) * (-dx * Sin(L) + dy * Cos(L))+ sin(B)/cos(b) * (1 - e2) * (wx * Cos(L) + wy * Sin(L)) - wz;
End;
Function WGS84Alt(Bd, Ld, H:double):double;
var B, L, N, dH:double;
begin
    B:=Bd * Pi / 180;
    L:= Ld * Pi / 180;
    N:= a / sqrt(1 - e2 * Sin(B)*sin(b));
    dH:= -a / N * da + N * Sin(B) *Sin(B) * de2 / 2 + (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B)- N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L))+ (a *a/ N + H) * ms;
    WGS84Alt:= H + dH;
End;
Function WGS84_SK42_Lat(Bd, Ld, H:double):double;  //Из wgs84 в Краковскую (широта)
begin    WGS84_SK42_Lat:=Bd - dB(Bd, Ld, H) / 3600;
End;
Function SK42_WGS84_Lat(Bd, Ld, H:double):double;
begin SK42_WGS84_Lat:= Bd + dB(Bd, Ld, H) / 3600;
End;
Function WGS84_SK42_Long(Bd, Ld, H:double):double;  //Из wgs84 в Краковскую (долгота)
begin    WGS84_SK42_Long:=Ld - dL(Bd, Ld, H) / 3600;
End;
Function SK42_WGS84_Long(Bd, Ld, H:double):double;
begin    SK42_WGS84_Long:= Ld + dL(Bd, Ld, H) / 3600;
End;
procedure TForm1.Button1Click(Sender: TObject);
var
f1: TextFile; // файл fName: String[80]; // имя файла
buf: string;
fname: String[80]; // буфер для чтения из файла
begin
fName := Edit1.Text;
Memo1.Lines.Clear; 
AssignFile(f1, fName);

{$I-}
Reset(f1); // открыть для чтения
{$I+}

if (IOResult <> 0)or(fname='') then
begin
  MessageDlg('Ошибка доступа к файлу ' + fName,
  mtError,[mbOk],0); exit;
end;

  while not EOF(f1) do begin
  readln(f1, buf); // прочитать строку из файла
//writeln(f2,buf);
  Memo1.Lines.Add( Utf8ToAnsi(buf));
end;
CloseFile(f1); // закрыть файл
end;  
procedure TForm1.Button2Click(Sender: TObject);
var  f2: TextFile; fname:string; i:integer;
begin
if filename='' then
begin
  ShowMessage('файл пуст');
end
else begin
assignfile(f2,filename+'.txt');
Rewrite(f2); // открыть для перезаписи
// запись в файл
for i:=0 to Memo2.Lines.Count do // строки нумеруются с нуля
writeln(f2, Memo2.Lines[i]);
CloseFile(f2); // закрыть файл
MessageDlg('Данные записаны в файл:'+filename+'.txt',mtInformation,[mbOk],0); end;
end;

procedure TForm1.N2Click(Sender: TObject);
begin
if opendialog1.Execute then
  begin
  filename:=opendialog1.filename;
  Edit1.Text:=filename;
  end;
end;

procedure TForm1.N4Click(Sender: TObject);
begin
close;
end;

procedure TForm1.N3Click(Sender: TObject);
begin
ShowMessage('Программа ');
end;
procedure TForm1.Button3Click(Sender: TObject);
var i,j,n:integer; s,s1,s2:string;s1f,s2f,r1,r2,r3,x,y,sk,lb:real;
b,l:real;
begin
i:=0; Memo2.Lines.Clear;
for i:=0 to Memo1.Lines.Count do // строки нумеруются с нуля
if (pos('lat="',Memo1.Lines[i])<>0)and(pos('lon="',Memo1.Lines[i])<>0)and(pos('maxlat',Memo1.Lines[i])=0)
then begin
//Memo2.Lines.Add(Memo1.Lines[i]); 
      s:= Memo1.Lines[i]; s1:='';s2:='';
      for j:= Pos('lat="',s)-1 to Pos('" ',s)-1 do s1:=s1+s[j];
      Delete(S1,1,6);
s1 := StringReplace(s1, '.', DecimalSeparator, [rfReplaceAll]);
s1 := StringReplace(s1, ',', DecimalSeparator, [rfReplaceAll]);
      s1f:=strtofloat(s1);
      for j:= Pos('lon="',s)-1 to Pos('">',s)-1 do s2:=s2+s[j];
      Delete(S2,1,6);
s2 := StringReplace(s2, '.', DecimalSeparator, [rfReplaceAll]);
s2 := StringReplace(s2, ',', DecimalSeparator, [rfReplaceAll]);
      s2f:=strtofloat(s2);
      b:=s1f;lb:= s2f;
sk:=(6+lb)/7;
n:=trunc(sk);

l:=(lb-(3+6*(n-1)))/57.29577951;
x:=6367558.4968*b-2*sin(b)*cos(b)*(16002.8900+66.9607*sin(b)*sin(b)+0.3515*sin(b)*sin(b)*sin(b)*sin(b)-l*l*(1594561.25+5336.535*sin(b)*sin(b)+26.790*sin(b)*sin(b)*sin(b)*sin(b)+0.149*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(672483.4-811219.9*sin(b)*sin(b)+5420*sin(b)*sin(b)*sin(b)*sin(b)-10.6*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(278194-830174*sin(b)*sin(b)+572434*sin(b)*sin(b)*sin(b)*sin(b)-16010*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(109500-574700*sin(b)*sin(b)+863700*sin(b)*sin(b)*sin(b)*sin(b)-398600*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b))))));
y:=(5+10*n)*100000+l*cos(b)*(6378245+21346.1415*sin(b)*sin(b)+107.1590*sin(b)*sin(b)*sin(b)*sin(b)+0.5977*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(1070204.16-2136826.66*sin(b)*sin(b)+17.98*sin(b)*sin(b)*sin(b)*sin(b)-11.99*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(270806-1523417*sin(b)*sin(b)+1327645*sin(b)*sin(b)*sin(b)*sin(b)-21701*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)+l*l*(79690-866190*sin(b)*sin(b)+1730360*sin(b)*sin(b)*sin(b)*sin(b)-945460*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)*sin(b)))));

      Memo2.Lines.Add(s1+' '+s2+' '+ 's1f='+floattostr(x)+' s2f='+floattostr(y));
      end;


s1:=floattostr(WGS84_SK42_Lat(50, 50, 0));
s2:=floattostr(y);

showmessage('s1:');
end;
end.