高斯坐标经纬度互相转换算法(Delphi)

减小字体 增大字体 作者:佚名  来源:转载  发布时间:2011-08-04 21:44:39

这个程序是根据网上找到的VC代码改写而成的Delphi库单元,经验算,比较准确,支持西安80及北京54。 
  本人原创,转载请保留本人信息。http://wallimn.javaeye.com 
  代码及使用方法如下(javaeye的博客代码高亮竟然不支持delphi,抗议): 

unit Translate; 
{ 
经纬度坐标与高斯-克吕格投影坐标的互算。 
时间:2009-05-11 
博客:http://wallimn.javaeye.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.javaeye.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.javaeye.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; 
{ 
引用此库单元,具体使用方法如下: 
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; 

} 
end. 

注:我的网络硬盘(http://wallimn.ys168.com)上放了这个转换工具,需要的请自行下载

下载了好几个源程序都不能用,在他们的基础上经过我与公式的核对.现在将代码改成了java

以下是完整代码:绝对能用,我现在的项目中就用的这个.

package tms.base.systemlib;
@SuppressWarnings("unchecked")
public class GaussXYDeal {
//  由高斯投影坐标反算成经纬度
public static double[] GaussToBL(double X, double Y)//, double *longitude, double *latitude)

{
  int ProjNo; int ZoneWide; ////带宽
  double[] output = new double[2];
  double longitude1,latitude1, longitude0, X0,Y0, xval,yval;//latitude0,
  double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;
  iPI = 0.0174532925199433; ////3.1415926535898/180.0;
  //a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数
  a=6378140.0; f=1/298.257; //80年西安坐标系参数
  ZoneWide = 6; ////6度带宽
  ProjNo = (int)(X/1000000L) ; //查找带号
  longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;
  longitude0 = longitude0 * iPI ; //中央经线
  
  
  X0 = ProjNo*1000000L+500000L;
  Y0 = 0;
  xval = X-X0; yval = Y-Y0; //带内大地坐标
  e2 = 2*f-f*f;
  e1 = (1.0-Math.sqrt(1-e2))/(1.0+Math.sqrt(1-e2));
  ee = e2/(1-e2);
  M = yval;
  u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));
  fai = u+(3*e1/2-27*e1*e1*e1/32)*Math.sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*Math.sin(
  4*u)
  +(151*e1*e1*e1/96)*Math.sin(6*u)+(1097*e1*e1*e1*e1/512)*Math.sin(8*u);
  C = ee*Math.cos(fai)*Math.cos(fai);
  T = Math.tan(fai)*Math.tan(fai);
  NN = a/Math.sqrt(1.0-e2*Math.sin(fai)*Math.sin(fai));
  R = a*(1-e2)/Math.sqrt((1-e2*Math.sin(fai)*Math.sin(fai))*(1-e2*Math.sin(fai)*Math.sin(fai))*(1-e2*Math.sin
  (fai)*Math.sin(fai)));
  D = xval/NN;
  //计算经度(Longitude) 纬度(Latitude)
  longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D
  *D*D*D*D/120)/Math.cos(fai);
  latitude1 = fai -(NN*Math.tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24
  +(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);
  //转换为度 DD
  output[0] = longitude1 / iPI;
  output[1] = latitude1 / iPI;
  return output;
  //*longitude = longitude1 / iPI;
  //*latitude = latitude1 / iPI;
}
////  由经纬度反算成高斯投影坐标
public void GaussToBLToGauss(double longitude, double latitude)
{
 int ProjNo=0; int ZoneWide; ////带宽
 double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;
 double a,f, e2,ee, NN, T,C,A, M, iPI;
 iPI = 0.0174532925199433; ////3.1415926535898/180.0;
 ZoneWide = 6; ////6度带宽
 a=6378245.0; f=1.0/298.3; //54年北京坐标系参数
 ////a=6378140.0; f=1/298.257; //80年西安坐标系参数
 ProjNo = (int)(longitude / ZoneWide) ;
 longitude0 = ProjNo * ZoneWide + ZoneWide / 2;
 longitude0 = longitude0 * iPI ;
 latitude0 = 0;
 System.out.println(latitude0);
 longitude1 = longitude * iPI ; //经度转换为弧度
 latitude1 = latitude * iPI ; //纬度转换为弧度
 e2=2*f-f*f;
 ee=e2*(1.0-e2);
 NN=a/Math.sqrt(1.0-e2*Math.sin(latitude1)*Math.sin(latitude1));
 T=Math.tan(latitude1)*Math.tan(latitude1);
 C=ee*Math.cos(latitude1)*Math.cos(latitude1);
 A=(longitude1-longitude0)*Math.cos(latitude1);
 M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2
 *e2/1024)*Math.sin(2*latitude1)
 +(15*e2*e2/256+45*e2*e2*e2/1024)*Math.sin(4*latitude1)-(35*e2*e2*e2/3072)*Math.sin(6*latitude1));
 xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
 yval = M+NN*Math.tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
 +(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
 X0 = 1000000L*(ProjNo+1)+500000L;
 Y0 = 0;
 xval = xval+X0; yval = yval+Y0;
 //*X = xval;
 //*Y = yval;
 System.out.println("x:"+xval);
 System.out.println("y:"+yval);

 }
}

Tags:

作者:佚名
  • 好的评价 如果您觉得此文章好,就请您
      0%(0)
  • 差的评价 如果您觉得此文章差,就请您
      0%(0)

文章评论评论内容只代表网友观点,与本站立场无关!

   评论摘要(共 0 条,得分 0 分,平均 0 分) 查看完整评论