下载文档原格式
![[转载]高斯正反算](https://uimg.taocdn.com/d13bd7fbe109581b6bd97f19227916888486b9df.webp)
[转载]⾼斯正反算⼤地坐标向笛卡尔坐标转换⾼斯正反算采⽤不同椭球实现⾼斯克⾥格投影，将经纬度坐标转换为⾼斯平⾯坐标：正算⾼斯平⾯坐标转换为不同椭球下的经纬度坐标：反算1void GaussProjectDirect(double a,double efang,double B,double L,double L0,double& x,double &y,double& R)//⾼斯投影正算克⽒2 {34double b=aefangtob(a,efang);5double e2=seconde(a,b);6double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f",W);7double N=a/W;printf("N=%f",N);8double M=a*(1-efang)/pow(W,3);printf("M=%f",M);9double t=tee(B);10double eitef=eitefang(a,b,B);11double l=L-L0;12//主曲率半径计算13double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;14 m0=a*(1-efang); n0=a;15 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;16 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;17 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;18 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;19//⼦午线曲率半径20double a0,a2,a4,a6,a8;21 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;22 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;23 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;24 a6=m6/32+m8/16;25 a8=m8/128;2627double X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8*sin(8*B);28 x=X+N/2*t*cos(B)*cos(B)*l*l+N/24*t*(5-t*t+9*eitef+4*pow(eitef,2))*pow(cos(B),4)*pow(l,4)+N/720*t*(61-58*t*t+pow(t,4))*pow(cos(B),6)*pow(l,6);29 y=N*cos(B)*l+N/6*(1-t*t+eitef)*pow(cos(B),3)*pow(l,3)+N/120*(5-18*t*t+pow(t,4)+14*eitef-58*eitef*t*t)*pow(cos(B),5)*pow(l,5);30 R=sqrt(M*N);31 }323334//⾼斯投影反算353637void GaussProjectInvert(double a,double efang,double x,double y,double L0,double &B,double& L,double& R)38 {39double b=aefangtob(a,efang);404142double m0,m2,m4,m6,m8,n0,n2,n4,n6,n8;43 m0=a*(1-efang); n0=a;44 m2=3.0/2.0*efang*m0; n2=1.0/2.0*efang*n0;45 m4=5.0/4.0*efang*m2; n4=3.0/4.0*efang*n2;46 m6=7.0/6.0*efang*m4; n6=5.0/6.0*efang*n4;47 m8=9.0/8.0*efang*m6; n8=7.0/8.0*efang*n6;484950//⼦午线曲率半径51double a0,a2,a4,a6,a8;52 a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;53 a2=m2/2+m4/2+15*m6/32+7.0/16.0*m8;54 a4=m4/8.0+3.0/16.0*m6+7.0/32.0*m8;55 a6=m6/32+m8/16;56 a8=m8/128;575859double X=x;60double FBf=0;61double Bf0=X/a0,Bf1=0;62while((Bf0-Bf1)>=0.0001)63 { Bf1=Bf0;64 FBf=a0*Bf0-a2/2*sin(2*Bf0)+a4/4*sin(4*Bf0)-a6/6*sin(6*Bf0)+a8/8*sin(8*Bf0);65 Bf0=(X-FBf)/a0;66 }67double Bf=Bf0;68double Vf=bigv(a,b,Bf);69double tf=tee(Bf);70double Nf=bign(a,b,Bf);71double eiteffang=eitefang(a,b,Bf);72double Bdu=rad_deg(Bf)-1/2.0*Vf*Vf*tf*(pow((y/Nf),2)-1.0/12*(5+3*tf*tf+eiteffang-9*eiteffang*tf*tf)*pow((y/Nf),4)+1.0/360.0*(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6))*180/PI; 73double ldu=1.0/cos(Bf)*(y/Nf+1.0/6.0*(1+2*tf*tf+eiteffang)*pow((y/Nf),3)+1.0/120.0*(5+28*tf*tf+24*tf*tf+6*eiteffang+8*eiteffang*tf*tf)*pow((y/Nf),5))*180.0/PI;747576 B=deg_int(Bdu);77 L=L0+deg_int(ldu);78double W=sqrt(1-efang*sin(B)*sin(B));printf("W=%f\n",W); 79double N=a/W;printf("N=%f\n",N);80double M=a*(1-efang)/pow(W,3);printf("M=%f\n",M);81 R=sqrt(M*N);828384 }。
GAUSSLE坐标正反算fx－5800程序源程序1.正算主程序Lbi 0：“G”？K：“Z=-1，Q=0，Y=1”？AIf A=-1:Then Prog “ZX”:Prog “GSZS”:IfEnd If A=1:Then Prog “YX”:Prog “GSZS”:IfEnd If A=0:Then Prog “QX”:Prog “GSZS”:IfEnd “X=”:X◢”Y=”:Y◢Goto 0说明：K 正算时所求点的里程A 选择线路，左幅=-1，右幅=1，整体式=0 正算子程序GSZS((P-R)÷(2(H-O)PR))→D：“JIAODU”?M:”JULI(-Z +Y)” ?L(Abs(K-O)) →J：Prog"SUB1"：(F-M) →FReturn2. 反算主程序 GSFSLbi 0：?X：?Y：X→Z[2]：Y→Z[3]：“QDXO”？I："QDY0"?S:"QDLC"?O:"QDFWJ "?G："ZDLC"?H："QDR"?P："ZDR"? R：”Q(Z=-1 ZX=0 Y=1)” ?Q:( (P-R)÷(2(H-O)PR)) →D:（Abs((Y-S)cos(G-90)－(X-I)sin(G-90)) ) →J：0→L：90→M:Lbl 1：Prog "SUB1"：((Z[3]-Y)cos(G-90+QJ(1÷P+JD)×180÷π)-(Z[2]-X)sin(G-90+QJ(1÷P +JD) ×180÷π)) →L：If:AbsL<10-6 Then Goto2：Else J+L→J：Go to 1:←┘Lbl 2：0→L：Prog "SUB1"：((Z[3]-Y)÷sinF)→L:”K=”: O+J→k◢”L=”:L→L◢Goto 03. 反算，正算子程序(SUB1)0.1184634425→A:0.2393143352→B: 0.2844444444→Z[4]:0.0469100770→C: 0.2307 653449→E: 0.5→Z[1]:(I+J(Acos(G+QCJ(1÷P+CJD)×180÷π)+Bcos(G+QEJ(1÷P+EJD)×180÷π)+Z[4]cos(G+Q Z[1]J(1÷P+Z[1]JD)×180÷π)+Bcos(G+Q(1-E)J(1÷P+(1-E)JD)×180÷π)+Acos(G+Q (1-C)J(1÷P+(1-C)JD) ×180÷π))) →X:(S+J(Asin(G+QCJ(1÷P+CJD)×180÷π)+Bsin(G+QEJ(1÷P+EJD)×180÷π)+Z[4]sin(G+QZ [1]J(1÷P+Z[1]JD)×180÷π)+Bsin(G+Q(1-E)J(1÷P+(1-E)JD)×180÷π)+Asin(G+Q (1-C)J (1÷P+(1-C)JD) ×180÷π))) →Y：(G+QJ(1÷P+JD) ×180÷π+M) →F：(X+LcosF)→X：(Y+LsinF) →Y4. 曲线元要素数据库：ZX/YX/QXIf K<(起点里程):Then Goto 2：IfEndIf K<（ZDLC）: Then QDXO →I：QDY0→S:QDLC→O:QDFWJ →G：ZDLC→H：QDR→P：ZDR→R：Q(Z=-1 ZX=0 Y=1)→Q: Goto 3：IfEnd……….(注：如有多个曲线元要素继续添加入数据库ZX/YX/QX中)Lbl 2:”NO”Lbl 3:Return说明：一、程序功能及原理1.功能说明：本程序由两个主程序——正算主程序(GSZS)、反算主程序(GSFS)和两个子程——正算子程序(SUB1)、线元数据库(DAT-M)构成，可以根据曲线段——直线、圆曲线、缓和曲线（完整或非完整型）的线元要素（起点坐标、起点里程、起点切线方位角、终点里程、起点曲率半径、止点曲率半径）及里程边距或坐标，对该曲线段范围内任意里程中边桩坐标进行正反算。
C#⾼斯正算⾼斯反算⾼斯换带等⾸先，你要确定椭球参数：C#代码1. a = 6378140; //西安80椭球 IGA752. e2 = 0.006694384999588;3. m0 = a * (1 - e2);4. m2 = 3.0 / 2 * e2 * m0;5. m4 = 5.0 / 4 * e2 * m2;6. m6 =7.0 / 6 * e2 * m4;7. m8 = 9.0 / 8 * e2 * m6;8. a0 = m0 + m2 / 2 + (3.0 / 8.0) * m4 + (5.0 / 16.0) * m6 + (35.0 / 128.0) * m8;9. a2 = m2 / 2 + m4 / 2 + 15.0 / 32 * m6 + 7.0 / 16 * m8;10. a4 = m4 / 8 + 3.0 / 16 * m6 + 7.0 / 32 * m8;11. a6 = m6 / 32 + m8 / 16;12. a8 = m8 / 128;13. xx = 0;14. yy = 0;15. _x = 0;16. _y = 0;17. BB = 0;18. LL = 0;下⾯才开始正题：⾼斯正算：把经纬度坐标转换为平⾯坐标C#代码1. void GaussPositive(double B, double L, double L0)2. {3. double X, t, N, h2, l, m, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;5. Bdu = (int)B;6. Bfen = (int)(B * 100) % 100;7. Bmiao = (B - Bdu - Bfen * 0.01) * 10000.0;8. B = Bdu * PI / 180.0 + (Bfen / 60.0) * PI / 180.0 + Bmiao / 3600.0 * PI / 180.0;9. Ldu = (int)L;10. Lfen = (int)(L * 100) % 100;11. Lmiao = (L - Ldu - Lfen * 0.01) * 10000.0;12. L = Ldu * PI / 180.0 + (Lfen / 60.0) * PI / 180 + Lmiao / 3600.0 * PI / 180.0;13. l = L - L0 * PI / 180;14. X = a0 * B - Math.Sin(B) * Math.Cos(B) * ((a2 - a4 + a6) + (2 * a4 - 16.0 / 3.0 * a6) * Math.Sin(B) * Math.Sin(B) + 16.0 / 3.0 * a6 * Math.Pow(Math.Sin(B), 4)) + a8 / 8.0 * Math.Sin(8 * B);15. t = Math.Tan(B);16. h2 = e2 / (1 - e2) * Math.Cos(B) * Math.Cos(B);17. N = a / Math.Sqrt(1 - e2 * Math.Sin(B) * Math.Sin(B));18. m = Math.Cos(B) * l;19. xx = X + N * t * ((0.5 + (1.0 / 24.0 * (5 - t * t + 9 * h2 + 4 * h2 * h2) + 1.0 / 720.0 * (61 - 58 * t * t + Math.Pow(t, 4)) * m * m) * m * m) * m * m);20. yy = N * ((1 + (1.0 / 6.0 * (1 - t * t + h2) + 1.0 / 120.0 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * h2 - 58 * h2 * t * t) * m * m) * m * m) * m);21. yy = yy + 500000;22. }⾼斯反算：把平⾯坐标转换为经纬度坐标：C#代码1. void GaussNegative(double x, double y, double L0)2.3. double Bf, Vf, l, tf, hf2, Nf, Bmiao, Lmiao;4. int Bdu, Bfen, Ldu, Lfen;4. int Bdu, Bfen, Ldu, Lfen;5. y = y - 500000;6. Bf = hcfansuan(x);7. Vf = Math.Sqrt(1 + e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf));8. tf = Math.Tan(Bf);9. hf2 = e2 / (1 - e2) * Math.Cos(Bf) * Math.Cos(Bf);10. Nf = a / Math.Sqrt(1 - e2 * Math.Sin(Bf) * Math.Sin(Bf));11. BB = (Bf - 0.5 * Vf * Vf * tf * (Math.Pow(y / Nf, 2) - 1.0 / 12 * (5 + 3 * tf * tf + hf2 - 9 * hf2 * tf * tf) * Math.Pow(y / Nf, 4) + 1.0 / 360 * (61 + 90 * tf * tf + 45 * tf * tf) * Math.Pow(y / Nf, 6))) * 180 / PI;12. Bdu = (int)BB;13. Bfen = (int)((BB - Bdu) * 60);14. Bmiao = ((BB - Bdu) * 60 - Bfen) * 60;15. BB = Bdu + 0.01 * Bfen + 0.0001 * Bmiao;16. l = 1.0 / Math.Cos(Bf) * (y / Nf - 1.0 / 6.0 * (1 + 2 * tf * tf + hf2) * Math.Pow(y / Nf, 3) + 1.0 / 120.0 * (5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * hf2 + 8 * hf2 * tf * tf) * Math.Pow(y / Nf, 5)) * 180.0 / PI;17. LL = L0 + l;18. Ldu = (int)LL;19. Lfen = (int)((LL - Ldu) * 60);20. Lmiao = ((LL - Ldu) * 60 - Lfen) * 60;21. LL = Ldu + 0.01 * Lfen + 0.0001 * Lmiao;⾥⾯涉及到的弧长反算：C#代码1. double hcfansuan(double pX)2. {3. double Bf0 = pX / a0;4. double Bf1, Bf2;5. Bf1 = Bf0;6. Bf2 = (pX - F(Bf1)) / a0;7. while ((Bf2 - Bf1) > 1.0E-11)8. {9. Bf1 = Bf2;10. Bf2 = (pX - F(Bf1)) / a0;11. }12. return Bf1;13. }⾼斯换带就⽐较简单了：C#代码1. void GaussZone(double x, double y, double L0, double L0new)2. {3. GaussNegative(x, y, L0);4. GaussPositive(BB, LL, L0new);5. }。
高斯坐标正反matlab算程序：function [x,y]=zgauss(B,L,B0)%B是弧度为单位%L B0都是一度为单位e2=0.006693421622966; %第一偏心率ee2=0.006738525414683; %第二偏心率a=6378245; %长半轴X=111134.8611*B0-16036.4803*sin(2*B)+16.8281*sin(4*B)-0.022*sin(6*B); %子午弧长p=3600*180/pi; %ρ"W=sqrt(1-e2*(sin(B)^2)); %WN=a/W; %卯酉圈半径t=tan(B);yita=ee2*cos(B)*cos(B); %η的平方l=L-117; %Lll=l*3600;x=X+(N*sin(B)*cos(B)*ll^2)/(2*p^2)+N*sin(B)*(cos(B)^3)*(5-t^2-9*yita) *(ll^4)/(24*p^4);%x坐标y=N*cos(B)*ll/p+N*(cos(B)^3)*(1-t^2+yita)*ll^3/(6*p^3)+N*(cos(B)^5)*( 5-18*t^2+t^4)*ll^5/(120*p^5);%坐标function [BB,l]=fgauss(x,y)%在克拉索夫椭球上e2=0.006693421622966; %第一偏心率ee2=0.006738525414683; %第二偏心率a=6378245; %长半轴a0=111134.8611;Bf(1)=x/a0;Bf(1)=Bf(1)*pi/180; %把度换算成弧度%以下用迭代法解出 Bffor i=1:5Bf(i+1)=(x-(-16036.4803*sin(2*Bf(i))+16.8281*sin(4*Bf(i))-0.022*sin(6 *Bf(i))))/a0*pi/180;end%代入公式计算BBf=Bf(5)；tf=tan(BBf);yitaf2=ee2*(cos(BBf)^2);Wf=sqrt(1-e2*(sin(BBf)^2));Mf=a*(1-e2)/Wf;Nf=a/Wf;BB=BBf-tf*y^2/(2*Mf*Nf)+tf^3*(5+3*tf^2+yitaf2-9*yitaf2*tf^2)*y^4/(24* Mf*Nf^3);BB=BB*180/pi;l=y/(Nf*cos(BBf))-(1+2*tf^2+yitaf2)*y^3/(6*Nf^3*cos(BBf))+(5+28*tf^2+ 24*tf^4)*y^5/(120*Nf^5*cos(BBf));我是一名在校大学生，大家在实际生产生活中若有什么关于matlab，和测量有什么问题欢迎来邮件与我联系。
#include<math.h>#include<stdio.h>#define PI 3.1415926535897932#define P 206264.806247096355void main(){long double AngleToRadian(long double alfa);long double RadianToAngle(long double alfa);long double Bf(long double a,long double b,long double x);void GSZS(long double a,long double b,long double l,long double B,long double *x,long double *y);void GSFS(long double a,long double b,long double Bf,long double y,long double *B,long double *l);char XZ,XZZF,XZSL; int N,num;long double a,b,B,L,x,y;char YN='O';long double L1,B1,L0,l;long double *pointer_x,*pointer_y;long double FSB,FSL;long double *pointer_B,*pointer_L;long double DH;pointer_B=&FSB;pointer_L=&FSL;pointer_x=&x;pointer_y=&y;printf("\n============================欢迎使用xx投影计算程序============================\n");printf(" 说明\n");printf("1.程序可实现在任意椭球上进行高斯三度带和六度带正算和反算;\n");printf("2.经纬度的输入输出格式为：度分秒,例如：113°05\'13.68466\"输入时为:113.051368466;\n");printf("3.y坐标的输入输出格式为：带号*10E6+500公里+自然坐标,例如：y=18637682.388,18为带号,自然坐标为137682.388;\n");printf("4.其余可按照提示完成,如有疑问可发送Email至gys_long@,我们将尽快为您解答。
高斯投影正反算与换带计算True BASIC程序欧龙;陈性义;欧阳平【摘要】实际测量工作中,经常需要进行高斯投影正、反算与换带计算,如果用手工完成,则计算量庞大,且容易出错.为此,利用True BASIC编写了计算程序,只需输入高斯平面坐标,就能在克拉索夫斯基和IUGG-1975参考椭球面上进行高斯投影正、反算和换带计算,还可将计算成果导入CASS中展点,或上传到全站仪内存中使用.【期刊名称】《铁道勘察》【年(卷),期】2006(032)005【总页数】4页(P12-15)【关键词】高斯投影;正反算;换带计算;True BASIC【作者】欧龙;陈性义;欧阳平【作者单位】中国地质大学,湖北武汉,430074;中国地质大学,湖北武汉,430074;中国地质大学,湖北武汉,430074【正文语种】中文【中图分类】U21 数学模型1.1 参考椭球面的高斯投影参考椭球面的高斯投影是指将地表的观测元素先投影到参考椭球面上(称为高斯归化[1])，再投影到高斯平面上(称为高斯投影改化[1])，这样就可以在高斯平面直角坐标系中进行测量平差计算。
在控制测量学中，由控制点的大地经纬度(L，B)计算其高斯平面坐标(x，y)，称为高斯投影正算；由高斯平面坐标(x，y)计算其大地经纬度(L，B)，称为高斯投影反算；由一个投影带的高斯平面坐标(x1，y1)计算其在另一个投影带的高斯平面坐标(x2，y2)，称为高斯投影换带计算。
1.2 参考椭球面的高斯投影正算公式设投影带的主子午线经度为L0，地表P点的经纬度为(L，B)，其高斯平面坐标为(x，y)，子午线收敛角为γ，经度差为l=L-L0，则有高斯投影正算公式[2](1)γ0为以度为单位的子午线收敛角。
式中它们均与高程修正值ΔH无关。
X为子午圈弧长，计算公式为(2)式中，B0是以度为单位的纬度值，5个系数的计算公式为1.3 参考椭球面的高斯投影反算公式高斯投影反算就是已知地面点P的高斯坐标x，y，求其大地坐标L，B。
摘要本设计主要阐述了高斯投影分带以及高斯投影坐标正、反算的推导公式，从而根据公式来编写基于VB语言基础上的换带及坐标转换程序。
作者系统介绍了测量中经常使用的坐标系以及地图投影的概念和高斯投影的具体含义，叙述了换带和临带计算的原因以及它们在运算时的原理、过程，详细叙述了在VB语言中实现的原理基础以及代码的编写设计。
在设计中根据高斯的正反算公式写出了基于VB语言的程序设计，其程序设计任务完成了由地理坐标向54平面坐标系和80平面坐标系转换的功能，以及由54坐标系和80坐标系向地理坐标系转换的功能，同时也有同一平面坐标系不同投影带之间的换带计算和同一平面坐标系相同投影带临带计算等相互转换的功能。
关键词：高斯投影、坐标正反算、换带计算、临带换算、程序设计5程序设计5.1界面设计本程序要实现的功能是根据所选择的椭球参数和指定的分带情况，将已知地理坐标或高斯投影坐标经正算和反算求得相应的高斯坐标和地理坐标，以及相应的换带计算和临带计算。
因此需要用一个框架控件来组织椭球参数、两个框架分别组织分带选择和换算方式选择，两个框架组织地理坐标和高斯坐标，三个命令按钮分别执行投影计算、换带和临带计算。
程序设计界面如图5-1[9]图5-1 高斯投影计算程序设计界面命令按钮属性设置表如表5-1表5-1 命令按钮属性设置表对象属性值Command1 Caption BL->xy Command1 Name cmdCalc Command2 Caption 6->3 Command2 Name cmdChange Command3 Caption 临带计算Command3 Name cmdNear选择椭球框架内控件的属性值表5-2表5-2 择椭球框架内控件的属性值单选按钮控件属性设置表5-35-3 单选按钮控件属性设置表5.2程序代码设计在这里主要介绍高斯投影坐标转换的正反算代码设计，完整的代码见附录1所示。
5.2.1投影计算过程的正算子过程代码设计①54系高斯投影正算子过程Public Sub Pro54()Dim ll#, N#, a0#, a4#, a6#, a3#, a5#, cosB#cosB = Cos(B)ll = L - DoToHu(L0)N = 6399698.902 - (21562.267 - (108.973 - 0.612 * cosB * cosB) * cosB * cosB) * cosB * cosBa0 = 32140.404 - (135.3302 - (0.7092 - 0.004 * cosB * cosB) * cosB * cosB) * cosB * cosBa4 = (0.25 + 0.00252 * cosB * cosB) * cosB * cosB - 0.04166a6 = (0.166 * cosB * cosB - 0.084) * cosB * cosBa3 = (0.3333333 + 0.001123 * cosB * cosB) * cosB * cosB - 0.1666667a5 = 0.0083 - (0.1667 - (0.1968 + 0.004 * cosB * cosB) * cosB * cosB) * cosB * cosBX = 6367558.4969 * B - (a0 - (0.5 + (a4 + a6 * ll * ll) * ll * ll) * ll * ll * N) * Sin(B) * cosBY = (1 + (a3 + a5 * ll * ll) * ll * ll) * ll * N * cosBEnd Sub②80系高斯投影正算子过程Public Sub Pro80()Dim ll#, N#, a0#, a4#, a6#, a3#, a5#, cosB#cosB = Cos(B)ll = L - DoToHu(L0)N = 6399596.652 - (21565.045 - (108.996 - 0.603 * cosB * cosB) * cosB * cosB) * cosB * cosBa0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * cosB * cosB) * cosB *cosB) * cosB * cosBa4 = (0.25 + 0.00253 * cosB * cosB) * cosB * cosB - 0.04167a6 = (0.167 * cosB * cosB - 0.083) * cosB * cosBa3 = (0.3333333 + 0.001123 * cosB * cosB) * cosB * cosB - 0.1666667a5 = 0.00878 - (0.1702 - 0.20382 * cosB * cosB) * cosB * cosBX = 6367452.1328 * B - (a0 - (0.5 + (a4 + a6 * ll * ll) * ll * ll) * ll * ll * N) * Sin(B) * cosBY = (1 + (a3 + a5 * ll * ll) * ll * ll) * ll * N * cosBEnd Sub5.2.2投影计算过程的反算子过程代码设计①54系高斯投影反算子过程[12]Public Sub ConPro54()Dim Bf#, bet#, Z#, Nf#, b2#, b3#, b4#, b5#, cos2B#, cos2Bf#bet = X / 6367558.4969cos2B = Cos(bet) * Cos(bet)Bf = bet + (50221746 + (293622 + (2350 + 22 * cos2B) * cos2B) * cos2B) * 0.0000000001 * Sin(bet) * Cos(bet)cos2Bf = Cos(Bf) * Cos(Bf)Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * cos2Bf) * cos2Bf) * cos2BfZ = Y / (Nf * Cos(Bf))b2 = (0.5 + 0.003369 * cos2Bf) * Sin(Bf) * Cos(Bf)b3 = 0.333333 - (0.166667 - 0.001123 * cos2Bf) * cos2Bfb4 = 0.25 + (0.16161 + 0.00562 * cos2Bf) * cos2Bfb5 = 0.2 - (0.1667 - 0.0088 * cos2Bf) * cos2BfB = Bf - (1 - (b4 - 0.12 * Z * Z) * Z * Z) * Z * Z * b2L = DoToHu(L0) + (1 - (b3 - b5 * Z * Z) * Z * Z) * ZEnd Sub②80系高斯投影反算子过程Public Sub ConPro80()Dim Bf#, bet#, Z#, Nf#, b2#, b3#, b4#, b5#, cos2B#, cos2Bf#bet = X / 6367558.4969cos2B = Cos(B) * Cos(B)Bf = bet + (50221746 + (293622 + (2350 + 22 * cos2B) * cos2B) * cos2B) * 0.0000000001 * Sin(bet) * Cos(bet)cos2Bf = Cos(Bf) * Cos(Bf)Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * cos2Bf) * cos2Bf) * cos2BfZ = Y / (Nf * Cos(Bf))b2 = (0.5 + 0.00336975 * cos2Bf) * Sin(Bf) * Cos(Bf)b3 = 0.333333 - (0.166667 - 0.001123 * cos2Bf) * cos2Bfb4 = 0.25 + (0.161612 + 0.005617 * cos2Bf) * cos2Bfb5 = 0.2 - (0.16667 - 0.00878 * cos2Bf) * cos2BfB = Bf - (1 - (b4 - 0.147 * Z * Z) * Z * Z) * Z * Z * b2L = DoToHu(L0) + (1 - (b3 - b5 * Z * Z) * Z * Z) * ZEnd Sub5.3程序的操作介绍下面以实例来介绍程序的操作步骤。
#include <iostream>#include <iomanip>#include <manth.h>#define PI 3.141592653#define a 6378137 //CGCS2000椭球长半轴using namespace std;double abs(double x){return x<0? -1*x;x;}double max(double x,double y){x=abs(x);y=abs(y);return x>y?x;y;}void main(){double e=sqrt(1-pow(1-1/298.257222101),2.0);double X,Y,Z；cout<<"请输入该点在直角坐标系中的X坐标值";cin>>X;cout<<"请输入该点在直角坐标系中的Y坐标值";cin>>Y;cout<<"请输入该点在直角坐标系中的Z坐标值";cin>>Z;double L;//定义大地经度LL=atan(X/Y)*180/PI;int i=0; //用于记录迭代次数double s=sqrt(pow(X,2.0)+pow(Y,2.0));double B,H;//定义大地维度B，大地高Hdouble B0,H0;//定义大地维度和大地高的初值double N0,M0,Z0,S0;//定义公式中的一些间接量double ds,dz,db,dh;//定义修正量ds=1.0,dz=1.0;//设定初始值while(max(ds,dz)>10e-5)//使精度大于10^（-5）时停止迭代{i++;B0=atan(Z/(sqrt(pow(X,2.0)+pow(Y,2.0))*(1-pow(e,2.0)));H0=sqrt(pow(s,2.0)+pow(Z+N0*pow(e,2.0)*sin(B0*PI/180),2.0))-N0;N0=a/sqrt(1-pow(e,2.0)*pow(sin(B0*PI/180),2.0));M0=a*(1-pow(e,2.0))/pow(sqrt(1-pow(e,2.0)*pow(sin(B0*PI/180),2.0)),3.0);S0=(N0+H0)*cos(B0*PI/180);Z0=(N0*(1-pow(e,2.0))+H0)*sin(B0*PI/180);ds=s-S0;dz=Z-Z0;db=-sin(B0*PI/180)/(M0+H0)*ds+cos(B0*PI/180)/(M0+H0)*dz;dh=cos(B0*PI/180)*ds+sin(B0*PI/180)*dz;B=B0=B0+db;H=H0=H0+dh;};cout<<"大地坐标（B,H,L）为：("<<setprecision(10)<<B<<","<<setprecision(10)<<H<","<<setprecision(10)<<H<<")"<<"共迭代"<<i<<"次"<<endl;}#include <iostream>#include <iomanip>#include <math.h>#define PI 3.141592653#define a 6378137#define b 6356752using namespace std;void main(){double e=sqrt(1-pow(1-1/298.257222101,2.0));double X,Y,Z;cout<<"请输入该点在直角坐标系中的X坐标值";cin>>X;cout<<"请输入该点在直角坐标系中的Y坐标值";cin>>Y;cout<<"请输入该点在直角坐标系中的Z坐标值";cin>>Z;double L;//定义大地经度LL=atan(Y/X)*180/PI;double N0,H0,B0;//迭代开始的初值double N,H,B;N=N0=a;H=H0=sqrt(pow(X,2.0)+pow(Y,2.0)+pow(Z,2.0))-sqrt(a)*sqrt(b);B=B0=atan(Z/(1-pow(e,2.0)*N0/(N0+H0))/sqrt(pow(X,2.0)+pow(Y,2.0)))*180/PI;int i=0;//记录迭代次数while(i<10000){N0=N;N=a/sqrt(1-pow(e,2.0)*pow(sin(B0*180/PI),2.0));H=sqrt(pow(X,2.0)+pow(Y,2.0))/cos(B0)-N0;B=atan(Z/(1-pow(e,2.0)*N0/(N0+H0))/sqrt(pow(X,2.0)+pow(Y,2.0)))*180/PI;i++;if((H-H0<=0.001||H-H0>=-0.001)&&(B-B0<=0.00001||B-B0>=-0.00001))break;}if(i==10000)cout<<"迭代失败"<<endl;elsecout<<"解为B="<<setprecision(20)<<B<<",H="<<setprecision(20)<<H<<",L="<<setprecision(20)<<L<<endl ;}。
高斯正反算源代码using System;using System.Collections.Generic;using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Windows.Forms;namespace 高斯正反算{public partial class Form1 : Form{public Form1(){InitializeComponent();}void a0a2a4a6a8(double a, double e, double[] Coeficient_a0){double m0, m2, m4, m6, m8;m0 = a * (1 - e * e);m2 = 3 * e * e * m0 / 2;m4 = 5 * e * e * m2 / 4;m6 = 7 * e * e * m4 / 6;m8 = 9 * e * e * m6 / 8;/*计算a0 a2 a4 a6 a8*/Coeficient_a0[0] = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128;Coeficient_a0[1] = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16;Coeficient_a0[2] = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32;Coeficient_a0[3] = m6 / 32 + m8 / 16;Coeficient_a0[4] = m8 / 128;}private void button1_Click(object sender, EventArgs e){int h=0,k=0;double a=0,Alfa=0;double dmslat,dmslon,dmsl0;double a0 ,a2 ,a4, a6, a8;double radlat,radlon,radl0,l;double b,t,sb,cb,ita,e1,ee;double X,l0;double N,c,v;double coor_x,coor_y;ouble Bf0,Bf1,dB,FBf,bf;double itaf, tf;double Nf,Mf;double B,L,dietaB,dietal;double sinBf,cosBf;double[] Coeficient_a0 = new double[5];try{if (radioButton4.Checked) { k = 1; }if (radioButton5.Checked) { k = 2; }if (k == 1) //正算{richTextBox1.Text = "";if (radioButton1.Checked == true) { h = 1; }if (radioButton2.Checked == true) { h = 2; }if (radioButton3.Checked == true) { h = 3; }// if (radioButton4.Checked == true) { h = 4; }if (h == 1) { a = 6378137; Alfa = 0.00335281066474; }//WGS-84if (h == 2) { a = 6378245; Alfa = 1.0 / 298.3; }//BJ54if (h == 3) { a = 6378140; Alfa = 1.0 / 298.257; }//gdz80/*输入已知数据:经度\纬度\ 中央子午线*/dmslat = double.Parse(textBox1.Text);dmslon = double.Parse(textBox2.Text);dmsl0 = double.Parse(textBox3.Text);/*将角度转化为弧度*/radlat = ANGLERAD.DMSTORAD(dmslat);radlon = ANGLERAD.DMSTORAD(dmslon);radl0 =ANGLERAD. DMSTORAD(dmsl0);l = radlon - radl0;/*计算椭球的基本参数和中间变量*/b = a * (1 - Alfa);sb = Math.Sin(radlat);cb = Math.Cos(radlat);t = sb / cb;e1 = Math.Sqrt((a / b) * (a / b) - 1);ee = Math.Sqrt(1 - (b / a) * (b / a));ita = e1 * cb;/*计算a0 a2 a4 a6 a8*/a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a4 = Coeficient_a0[2];a6 = Coeficient_a0[3];a8 = Coeficient_a0[4];/*计算X*/X = a0 * radlat - sb * cb * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3) * sb * sb + 16 * a6 * sb * sb * sb * sb / 3.0);/*计算卯酉圈半径N*/c = a * a / b;v = Math.Sqrt(1 + e1 * e1 * cb * cb);N = c / v;/*计算未知点的坐标*/coor_x = X + N * sb * cb * l * l / 2 + N * sb * cb * cb * cb * (5 - t * t + 9 * ita * ita + 4 * ita * ita * ita * ita) * l * l * l * l / 24 + N * sb * cb * cb * cb * cb * cb * (61 - 58 * t * t + t * t * t * t) * l * l * l * l * l * l / 720;coor_y = N * cb * l + N * cb * cb * cb * (1 - t * t + ita * ita) * l * l * l / 6 + N * cb * cb * cb * cb * cb * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * l * l * l * l * l / 120;/*输出未知点坐标*/richTextBox1.Text = "x:" + coor_x.ToString() + "\r\n" + "y:" + coor_y.ToString();}if (k == 2)//反算{richTextBox1.Text = "";if (radioButton1.Checked == true) { h = 1; }if (radioButton2.Checked == true) { h = 2; }if (radioButton3.Checked == true) { h = 3; }if (radioButton4.Checked == true) { h = 4; }if (h == 1) { a = 6378137; Alfa = 1.0 / 298.257223563; }if (h == 2) { a = 6378245; Alfa = 1.0 / 298.3; }if (h == 3) { a = 6378140; Alfa = 1.0 / 298.257; }/*输入l0，已知点坐标*/l0 = double.Parse(textBox6.Text); l0 = l0 *Math. PI / 180.0;coor_x = double.Parse(textBox4.Text);coor_y = double.Parse(textBox5.Text);/*计算b,e1,e*/b = a * (1 - Alfa);e1 = Math.Sqrt((a / b) * (a / b) - 1);ee = Math.Sqrt(1 - (b / a) * (b / a));/*计算a0 a2 a4 a6 a8*/a0a2a4a6a8(a, ee, Coeficient_a0);a0 = Coeficient_a0[0];a2 = Coeficient_a0[1];a4 = Coeficient_a0[2];a6 = Coeficient_a0[3];a8 = Coeficient_a0[4];X = coor_x;Bf0 = X / a0;do{sinBf = Math.Sin(Bf0); cosBf = Math.Cos(Bf0);FBf = -sinBf * cosBf * ((a2 - a4 + a6) + (2.0 * a4 - 16.0 * a6 / 3.0) * sinBf * sinBf +(16.0 / 3.0) * a6 * (sinBf * sinBf * sinBf * sinBf));Bf1 = (X - FBf) / a0;dB = Bf1 - Bf0;Bf0 = Bf1;}while (Math.Abs(dB * 180.0 / Math.PI * 3600.0) > 0.00001);bf = Bf1;/*计算c,v,N,M*/c = a * a / b;v = Math.Sqrt(1 + e1 * e1 * Math.Cos(bf) * Math.Cos(bf));Nf = c / v;Mf = c / (v * v * v);tf = Math.Sin(bf) / Math.Cos(bf);itaf = e1 * Math.Cos(bf);/*计算dietaB,dietal*/dietaB = tf * coor_y * coor_y / (2 * Mf * Nf) - tf * (5 + 3 * tf * tf + itaf * itaf - 9 * tf * tf * itaf * itaf) * coor_y * coor_y * coor_y * coor_y / (24 * Mf * Nf * Nf * Nf) + (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * coor_y * coor_y * coor_y * coor_y * coor_y * coor_y / (720 * Mf * Nf * Nf * Nf * Nf * Nf);dietal = coor_y / (Nf * Math.Cos(bf) + (1 + 2 * tf * tf + itaf * itaf) * Math.Cos(bf) * coor_y * coor_y / (6 * Nf)) + (5 + 44 * tf * tf + 32 * tf * tf * tf * tf - 2 * itaf * itaf - 16 * itaf * itaf * tf * tf) / (360 * Nf * Nf * Nf * Mf * Mf * Math.Cos(bf));B = bf - dietaB;L = l0 + dietal;dmslat = ANGLERAD.RADTODMS(B);dmslon = ANGLERAD.RADTODMS(L);richTextBox1.Text = "已知点的纬度:" + dmslat.ToString() + "\r\n" + "已知点的经度:" + dmslon.ToString();}}catch{MessageBox.Show("请检查！", "提示框");return;}}private void textBox5_TextChanged(object sender, EventArgs e){}}}。
高斯投影正、反算代码下载作者：boy转贴文章来源：/ShowPost.asp?id=25701 点击数：3295 更新时间：2006-7-10//高斯投影正、反算//////6度带宽54年北京坐标系//高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y){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;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*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)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*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*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;}//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude){int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;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-sqrt(1-e2))/(1.0+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)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*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)/cos(fai);latitude1 = fai -(NN*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*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}豆豆提供的程序里几个重要参数的意义NN卯酉圈曲率半径，测量学里面用N表示M为子午线弧长，测量学里用大X表示fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示R为底点所对的曲率半径，测量学里用Nf表示本篇文章来源于GIS空间站转载请以链接形式注明出处网址：/Article/168.htm。
