Формулы пересчета данных GPS-измерений из WGS-84 в СК-42
-
- Интересующийся
- Сообщения: 21
- Зарегистрирован: 25 ноя 2011, 19:18
- Репутация: 3
Re: Формулы пересчета данных 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
Не совсем уверен в том, что ГОСТ прошитый в Мапинфо (с датумом 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
-
- Интересующийся
- Сообщения: 40
- Зарегистрирован: 19 янв 2012, 10:03
- Репутация: 3
Re: Формулы пересчета данных GPS-измерений из WGS-84 в СК-42
Не совсем понял, а что делает функция 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
-
- Новоприбывший
- Сообщения: 9
- Зарегистрирован: 04 июн 2009, 12:04
- Репутация: 1
Re: Формулы пересчета данных GPS-измерений из WGS-84 в СК-42
Попытался реализовать пересчет из WGS-84 в СК-95. Поправил параметры (там задействованы эллепсойд 42, повороты и масштабный коэффициент) получилась ошибка по широте 11 км. Обнулил масштабный коэффициент (примерно равный 0,5) - точность возросла до +-2см, но всё равно как-то боязно использовать в реальных расчетах.
С высоткой дела вообще плохо: погрешность 1.3 м. Может Cocodrilo или другие знатоки найдут ошибку. Буду благодарен хотя-бы за ссылки на первоисточники.
Мои мысли по теме: ввиду высокой (порядка 5-200м) деформации СК42 (некачественной закладки из-за отсутсвия спутникового оборудования в 50-70г прошлого века) точных параметров пересчета на территорию всей страны не существует. Можно рассчитать точные параметры для любого небольшого района (100-200 кв.км.), но расширяя этот район, мы начнем терять точность.Спасти положение могут матрицы типа Канадской NTv2, но насколько я знаю создавать их на всю страну никто не собирается.
С высоткой дела вообще плохо: погрешность 1.3 м. Может Cocodrilo или другие знатоки найдут ошибку. Буду благодарен хотя-бы за ссылки на первоисточники.
Мои мысли по теме: ввиду высокой (порядка 5-200м) деформации СК42 (некачественной закладки из-за отсутсвия спутникового оборудования в 50-70г прошлого века) точных параметров пересчета на территорию всей страны не существует. Можно рассчитать точные параметры для любого небольшого района (100-200 кв.км.), но расширяя этот район, мы начнем терять точность.Спасти положение могут матрицы типа Канадской NTv2, но насколько я знаю создавать их на всю страну никто не собирается.
-
- Гуру
- Сообщения: 3058
- Зарегистрирован: 19 май 2010, 19:44
- Репутация: 189
Re: Формулы пересчета данных GPS-измерений из WGS-84 в СК-42
Кто-то создает и торгует NTv2 (где-то как бесплатный образец встречалось на Новосибирскую область), но все равно точность лучше 0.4 - 0.5 м вряд ли будет...
-
- Новоприбывший
- Сообщения: 6
- Зарегистрирован: 14 апр 2016, 04:58
- Репутация: 0
Re: Формулы пересчета данных GPS-измерений из WGS-84 в СК-42
Добрый день!
Вроде сделал все правильно, но результат от WGS84_SK42_Lat(50;50;0) выходит совсем другой 49,9998280965696, а не 49.99980414. Формулы перепроверил 10 раз.
Вроде сделал все правильно, но результат от 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.
Кто сейчас на конференции
Сейчас этот форум просматривают: нет зарегистрированных пользователей и 10 гостей