高斯座標經緯度互相轉換演算法(Delphi)
阿新 • • 發佈:2019-01-04
這個程式是根據網上找到的VC程式碼改寫而成的Delphi庫單元,經驗算,比較準確,支援西安80及北京54。
程式碼及使用方法如下:
unit Translate; { 經緯度座標與高斯-克呂格投影座標的互算。 時間:2009-05-11 部落格:http://wallimn.iteye.com 轉載請保留此資訊 } interface uses Math; type TTranslate = class(TObject) protected a,f,e2,e12:double; A1,A2,A3,A4:double; private L0:double; // 中央子午線經度 public procedure BL2xy(B,L :double; var x,y :double); procedure xy2BL(x,y :double; var B,L :double); procedure SetL0(dL0:double); end; TTranslate_Krasovsky = class(TTranslate) public constructor Create; end; TTranslate_IUGG1975=class(TTranslate) public constructor Create; end; function Dms2Rad( Dms:double) : double ; function Rad2Dms( Rad:double) : double ; implementation { 將度、分、秒形式轉化成弧度 } function Dms2Rad( Dms:double) : double ; var Degree,Miniute,Second:Double; Rad:Double; Sign:Integer; begin if(Dms >= 0) then Sign := 1 else Sign := -1; Dms := abs(Dms); Degree := floor(Dms); Miniute := floor(Dms * 100) mod 100; Second := floor(Dms * 10000) mod 100; Rad := Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * PI / 180.0; result:= Rad; end; { 將弧度轉換成度、分、秒錶示形式; 轉換的結果是度, } function Rad2Dms( Rad:double) : double ; var Degree, Miniute, Second:double; Sign:integer; begin if(Rad >= 0) then Sign := 1 else Sign := -1; Rad := abs(Rad * 180.0 / PI); Degree := floor(Rad); Miniute := floor(Rad * 60) mod 60; Second := floor(Rad * 3600) mod 60; Result := Sign * (Degree + Miniute / 100.0 + Second / 10000.0); //Result:=Rad * 180/PI; end; { TTranslate } { B,L 為以度為單位的緯度及經度 x,y 為轉換結果,即投影座標,其中y不帶帶號 時間:2009-05-11 部落格:http://wallimn.iteye.com } procedure TTranslate.BL2xy(B, L: double; var x, y: double); var XX, N, t, t2, m, m2, ng2:double; sinB, cosB:double; begin B:= B*PI/180.0; L:= L*PI/180.0; XX := A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 * B); sinB := sin(B); cosB := cos(B); t := tan(B); t2 := t * t; N := a / sqrt(1 - e2 * sinB * sinB); m := cosB * (L - L0); m2 := m * m; ng2 := cosB * cosB * e2 / (1 - e2); //x,y的計算公式見孔祥元等主編武漢大學出版社2002年出版的《控制測量學》 x := XX + N * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2); y := N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0)); y := y + 500000; end; { 設定中央子午線的經度,以度為單位 } procedure TTranslate.SetL0(dL0: double); begin //L0:= Dms2Rad(dL0); L0:=dL0*PI/180.0; end; { x,y 投影座標,其中y不帶帶號 B,L 為轉換結果,以度為單位的緯度及經度 時間:2009-05-11 部落格:http://wallimn.iteye.com } procedure TTranslate.xy2BL(x, y: double; var B, L: double); var sinB, cosB, t, t2, N ,ng2, V, yN:double; preB0, B0:double; eta:double; begin y := y- 500000; B0 := x / A1; repeat begin preB0 := B0; B0 := B0 * PI / 180.0; B0 := (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0))) / A1; eta := abs(B0 - preB0); end until(eta <= 0.000000001); B0 := B0 * PI / 180.0; B := Rad2Dms(B0); sinB := sin(B0); cosB := cos(B0); t := tan(B0); t2 := t * t; N := a / sqrt(1 - e2 * sinB * sinB); ng2 := cosB * cosB * e2 / (1 - e2); V := sqrt(1 + ng2); yN := y / N; B := B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2; L := L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB; //B:=Rad2Dms(B); //L:=Rad2Dms(L); B:=B*180.0/PI; L:=L*180.0/PI; end; { TTranslate_Krasovsky } constructor TTranslate_Krasovsky.Create; begin a := 6378245; f := 298.3; e2 := 1 - ((f - 1) / f) * ((f - 1) / f); e12 := (f / (f - 1)) * (f / (f - 1)) - 1; A1 := 111134.8611; A2 := -16036.4803; A3 := 16.8281; A4 := -0.0220; end; { TTranslate_IUGG1975 } constructor TTranslate_IUGG1975.Create; begin a := 6378140; f := 298.257; e2 := 1 - ((f - 1) / f) * ((f - 1) / f); e12 := (f / (f - 1)) * (f / (f - 1)) - 1; A1 := 111133.0047; //這幾個A是什麼意思? A2 := -16038.5282; A3 := 16.8326; A4 := -0.0220; end; end.
引用此庫單元,具體使用方法如下:
procedure TFrmMain.Button1Click(Sender: TObject);
var
t:TTranslate;
L,B:Double;
begin
t :=TTranslate_IUGG1975.create;
t.SetL0(111);
t.xy2BL(strToFloat(edtX.text),strToFloat(edtY.text),L,B);
showmessage('L='+FloatToStr(L)+' B='+FloatToStr(B));
//執行結果:L=20,B=109.15
end;