高斯投影正、反算 代码

 

[cpp] view plain copy

print?

1. <span style="font-size:13px;">   
2. //高斯投影正、反算   
3. //////6度带宽 54年北京坐标系   
4. //高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)   
5.   
6. void GaussProjCal(double longitude, double latitude, double *X, double *Y)   
7. {  
8.  int ProjNo=0; int ZoneWide; ////带宽   
9.  double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;  
10.  double a,f, e2,ee, NN, T,C,A, M, iPI;  
11.  iPI = 0.0174532925199433; ////3.1415926535898/180.0;   
12.  ZoneWide = 6; ////6度带宽   
13.  a=6378245.0; f=1.0/298.3; //54年北京坐标系参数   
14.  ////a=6378140.0; f=1/298.257; //80年西安坐标系参数   
15.  ProjNo = (int)(longitude / ZoneWide) ;   
16.  longitude0 = ProjNo * ZoneWide + ZoneWide / 2;   
17.  longitude0 = longitude0 * iPI ;  
18.  latitude0=0;   
19.  longitude1 = longitude * iPI ; //经度转换为弧度  
20.  latitude1 = latitude * iPI ; //纬度转换为弧度  
21.  e2=2*f-f*f;  
22.  ee=e2*(1.0-e2);  
23.  NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));  
24.  T=tan(latitude1)*tan(latitude1);  
25.  C=ee*cos(latitude1)*cos(latitude1);  
26.  A=(longitude1-longitude0)*cos(latitude1);   
27.   
28.   
29.  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  
30.   *e2/1024)*sin(2*latitude1)  
31.   +(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*sin(6*l  
32.   atitude1));  
33.  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);  
34.  yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24  
35.   +(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);  
36.  X0 = 1000000L*(ProjNo+1)+500000L;   
37.  Y0 = 0;   
38.  xval = xval+X0; yval = yval+Y0;   
39.  *X = xval;  
40.  *Y = yval;  
41. }   
42.   
43.    
44.   
45.   
46. //高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)  
47. void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) 字串9   
48.   
49. {  
50.  int ProjNo; int ZoneWide; ////带宽   
51.  double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;  
52.  double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;  
53.  iPI = 0.0174532925199433; ////3.1415926535898/180.0;   
54.  a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数   
55.  ////a=6378140.0; f=1/298.257; //80年西安坐标系参数   
56.  ZoneWide = 6; ////6度带宽   
57.  ProjNo = (int)(X/1000000L) ; //查找带号  
58.  longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;   
59.  longitude0 = longitude0 * iPI ; //中央经线  
60.  X0 = ProjNo*1000000L+500000L;   
61.  Y0 = 0;   
62.  xval = X-X0; yval = Y-Y0; //带内大地坐标  
63.  e2 = 2*f-f*f;  
64.  e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));  
65.  ee = e2/(1-e2);  
66.  M = yval;  
67.  u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));  
68.  fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(  
69.   4*u)  
70.   +(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);  
71.  C = ee*cos(fai)*cos(fai);  
72.  T = tan(fai)*tan(fai);  
73.  NN = a/sqrt(1.0-e2*sin(fai)*sin(fai)); 字串1   
74.   R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin  
75.   (fai)*sin(fai)));  
76.  D = xval/NN;  
77.  //计算经度(Longitude) 纬度(Latitude)  
78.  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  
79.   *D*D*D*D/120)/cos(fai);  
80.  latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24  
81.   +(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);   
82.  //转换为度 DD  
83.  *longitude = longitude1 / iPI;   
84.  *latitude = latitude1 / iPI;  
85. }   
86.   
87.    
88.   
89.   
90. //NN卯酉圈曲率半径，测量学里面用N表示   
91.   
92. //M为子午线弧长，测量学里用大X表示 字串2   
93.   
94. //fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示 字串4   
95.   
96. //R为底点所对的曲率半径，测量学里用Nf表示  
97. </span>