高斯投影坐标正反算VB程序
坐标反算(VB编程代码)
计算△X，△Y 的值
若△X，△Y 有为零 的，则直接判断方位 角。
若△X，△Y 不为零， 则计算象限角 R。
由象限角判断方 位角
输出方位角α
坐标反Sub Command3_Click() Dim Xa, Ya, Xb, Yb, M, N, Rab, R, F, R1, R2, R3 As Single pi = 3.1415926 Xa = Text1.Text Ya = Text4.Text Xb = Text2.Text Yb = Text5.Text M = Xb - Xa '求纵坐标增量 N = Yb - Ya '求横坐标增量 Rab = Math.Atn(Abs(N / M)) * 180 / pi '计算象限角 If M > 0 And N > 0 Then F = Rab '由象限角判断方位角 If M < 0 And N > 0 Then F = 180 - Rab If M < 0 And N < 0 Then F = 180 + Rab If M > 0 And N < 0 Then F = 360 - Rab If M = 0 And N > 0 Then F = 90 If M = 0 And N < 0 Then F = 270 If N = 0 And M > 0 Then F = 0 If N = 0 And M < 0 Then F = 180 R1 = Fix(F) '把弧度化为度 R2 = Fix((F - R1) * 60) R3 = Fix((((F - R1) * 60) - R2) * 60) Text3.Text = R1 & "°" & R2 & "′" & R3 & "″" Text6 = Sqr(M ^ 2 + N ^ 2) End Sub Private Sub Command2_Click() Text1 = "" Text2 = "" Text3 = "" Text4 = "" Text5 = "" Text6 = "" End Sub
VB实现高斯正反算
••Βιβλιοθήκη •• •• •
•
高斯－克吕格投影分带规定:该投影是国家基本比例尺地形图的数学基础,为控 制变形,采用分带投影的方法,在比例尺 1:2.5万-1:50万图上采用6°分带,对比例 尺为 1:1万及大于1:1万的图采用3°分带。 6°分带法:从格林威治零度经线起,每6°分为一个投影带,全球共分为60个 投影带,东半球从东经0°-6°为第一带,中央经线为3°,依此类推,投影带号为130。其投影代号n和中央经线经度L0的计算公式为:L0=(6n-3)°;西半球投影带从 180°回算到0°,编号为31-60,投影代号n和中央经线经度L0的计算公式为 L0=360-(6n-3)°。 3°分带法:从东经1°30′起,每3°为一带,将全球划分为120个投影带,东经 1°30′-4°30′,...178°30′-西经178°30′,...1°30′-东经1°30′。 东半球有60个投影带,编号1-60,各带中央经线计算公式:L0=3°n ,中央经线 为3°、6°...180°。西半球有60个投影带,编号1-60,各带中央经线计算公 式:L0=360°-3°n ,中央经线为西经177°、...3°、0°。 我国规定将各带纵坐标轴西移500公里,即将所有y值加上500公里,坐标值前再加 各带带号以18带为例,原坐标值为y=243353.5m,西移后为y=743353.5,加带号通 用坐标为y=18743353.5 。
昆明冶金高等专科学校2010届 毕业设计
VB实现高斯正反算
系 部：测计学院 专业班级：测绘0729班 学生姓名：张智宁 学生学号：0700001955 指导教师：欧阳老师 完成时间：2010年6月5日
前言
• 在日益发展的社会中，能很好很快的完成各种任 务才能取得更多的发展机会。用VB实现高斯正反 算，就能很快实现高斯坐标系和大地坐标系的转 换，为我们的计算和测量工作提供高效的服务。 • 此次设计要求我们利用大学所学的VB编程和控制 测量中的高斯平面坐标系为基础，用VB实现高斯 正反算，做到能够完整写出VB实现高斯正反算的 代码，而且还能生成运行文件进行高斯正反算。 同时，也让我们体会到编程在测量中应用之广， 在以后的工作中应该博学多思 。
vb课程课件测绘程序设计8(八)解析教学提纲
W 1e2si2nB
Ma(1W3e2) Vc3
M─子午圈曲率半径
t tan B 2 e 2 cos 2 B
e2
a2 b2 b2
N a 2 (1 e 2 sin 2 B )
四、工程测量投影面与投影带选择
3、地图投影的分类
等角投影——投影后角度不变，保持小范围内图形相似。 等面积投影——用于某些专题地图，投影后面积不变。 平面投影——投影平面与椭球面在某一点相切，按数学投影建立函数关系。 圆锥面投影——圆锥面与椭球体在某一纬圈相切或某两纬圈相割，按数学 投影。 圆柱面投影——圆柱面或椭圆柱面与椭球面在赤道或某一子午面上相切， 按数字投影。 正轴投影——圆柱面中心轴与椭球短轴重合，圆柱面与赤道相切。 横轴投影——圆柱面中心轴与椭球长轴重合，圆柱面与某一子午圈相切。 斜轴投影——圆柱面中心轴与椭球长、短轴都不重合，位于两者之间。
式中:
A 1 3 e 2 45 e 4 175 e 6 11025 e 8 43659 e10 693693 e12
4 64 256 16384
65536
1048576
B 3 e 2 15 e 4 525 e 6 2205 e 8 72765 e10 297297 e12
6°带带号N和中央子午线经度 LN的关系式：LN=6N-3 3°带带号n和中央子午线经度 Ln的关系式：Ln=3n
6°带与3°带带号之间的关系为 n2N1
三、高斯投影正算和反算
高斯投影坐标正算——由（B，L）求（x，y） 高斯投影坐标反算——由（x，y）求（B，L）
1、由（B，L）计算（x，y）－－正算
1048576
E
315 e 8 3465 e10 99099 e12
高斯UTM投影坐标正、反算
高斯UTM投影坐标正、反算测量每天不厌其烦的发招聘信息，图文教程给你一、高斯/UTM投影坐标正反算L、B——大地经、纬度坐标X、Y——高斯/UTM投影平面直角坐标高斯/UTM投影坐标正算：LB→XY高斯/UTM投影坐标反算：XY→LB序号点名大地坐标（度分秒）高斯平面直角坐标（m）L B X Y1 二塘大圆113 05 13.68466 22 38 49.52623 2507005.7794 714567.78812 渠黎医院113 04 55.00186 22 31 54.96727 2494239.9972 714212.25703 硝铵113 03 20.21577 22 38 39.70442 2506658.3680711330.6495厂4 扶绥中学113 02 42.90165 22 38 03.52353 2505530.1986 710280.04165 氮肥厂113 02 29.79765 22 36 54.35915 2503396.5128 709934.87696 渠黎中学113 04 02.83240 22 36 52.28043 2503369.2088 712593.96247 龙须塘113 03 43.97260 22 38 09.06116 2505724.6996 712022.36108 滑石公司113 05 42.00255 22 38 00.70728 2505514.7026 715397.8721 【思考和计算】工程中直角坐标和高程值通常取位到1mm或0.1mm，按照同样的精度要求，大地坐标经纬度应如何取位？取地球平均半径6371km，容易得出以下角度对应的弧长：1°→111.2km1′→1.85km0.01°→1.11km1″→31m0.00001°→1m0.1″→3m0.0000001°→1cm0.001″→3cm0.00000001°→1mm0.0001″→3mm0.000000001°→0.1mm0.00001″→0.3mm故：（1）以度分秒为单位，秒数值保留4~5位小数（2）以度为单位，度数值保留8~9位小数一、高斯投影坐标正反算程序高斯投影坐标正反算计算公式非常复杂，手工计算非常繁琐，一般使用程序来计算。
高斯正反算代码 vb
a4 = m4 / 8 + 3 / 16 * m6 + 7 / 32 * m8
a6 = m6 / 32 + m8 / 16
r = x / a0
π = 4 * Atn(1)
If Option1.Value = True And Option3.Value = True Then '国家2000坐标正算
p2 = 0 - a2 / (2 * a0)
p4 = a4 / (a0 * 4)
p6 = 0 - a6 / (6 * a0)
q2 = 0.5 * p2 ^ 3 - p2 - p2 * p4
x3 = 1 / 24 * x1 + 1 / 720 * x2 * m ^ 2
y1 = 1 - t ^ 2 + G ^ 2
y2 = 5 - 18 * t ^ 2 + t ^ 4 + 14 * G ^ 2 - 58 * G ^ 2 * t ^ 2
m = Cos(B) * L
m0 = a * (1 - e ^ 2)
m2 = 1.5 * e ^ 2 * m0
m4 = 5 / 4 * e ^ 2 * m2
m6 = 7 / 6 * e ^ 2 * m4
m2 = 1.5 * e ^ 2 * m0
m4 = 5 / 4 * e ^ 2 * m2
m6 = 7 / 6 * e ^ 2 * m4
m8 = 9 / 8 * e ^ 2 * m6
a0 = m0 + m2 / 2 + 3 / 8 * m4 + 5 / 16 * m6 + 35 / 128 * m8
Gauss投影正反算程序
#include<stdio.h>#include<math.h>#include<stdlib.h>double trans1() //由度分秒到弧度{double B1,B11,B12,B13,B111;scanf("%lf°%lf′%lf″",&B11,&B12,&B13);B111=fabs(B11); //B11可能为负值B1=B111+B12/60.0+B13/3600.0;B1=B1*atan(1)/45.0;if(B11<0)B1=-B1;return B1;}void trans2(double B) //由弧度到度分秒并输出角度值{int a,b;double B0;B0=fabs(B); //B可能为负值double c;B0=B0*45.0/atan(1);a=int(B0);b=int((B0-a)*60);c=(B0-a)*3600-b*60;if((int)(c)==60) //为了避免出现59′60″这种情况，不过好像不起作用，不知道为什么,原来是int没有加括号{b=b+1;c=0.0;}if(b==60){b=0;a=a+1;}if(B<0)a=-a;printf("%d°%d′%.4f″\n",a,b,c);}void Gauss1() //高斯投影坐标正算，克拉索夫斯基椭球{double B,L; //输入经纬度并转换成弧度printf("请输入大地纬度B=\n");B=trans1();printf("请输入大地经度L=\n");L=trans1();double L0,l;printf("请以度分秒的形式（如10°10′10″）输入中央子午线L0=\n");L0=trans1();l=L-L0; //经度与中央子午线之差……坑爹的书上居然又错了!!114°20′0″-117°0′0″=3°20′0″FUCK! 原来是我误会了……人家的L0=111°double R; //平面子午线收敛角γ（仅与B和l有关）R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l*l)*l*l*cos(B)*cos(B))* l*sin(B);double N,a0,a4,a6,a3,a5; //求解一些常数以及系数N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);double x,y; //求解Gauss投影平面坐标想，x,yx=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*cos(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);printf("平面坐标系中的坐标\nx=%.4fM\n",x);printf("平面坐标系中的坐标\ny=%.4fM\n",y);printf("平面子午线收敛角R=\n");trans2(R);}void Gauss2() //Gauss投影坐标反算，克拉索夫斯基椭球{double x,y,L0; //输入平面坐标x,y及中央子午线L0=printf("请输入平面坐标系中的坐标x=\n");scanf("%lf",&x);printf("请输入平面坐标系中的坐标y=\n");scanf("%lf",&y);printf("请以度分秒的形式输入（如10°10′10″）输入中央子午线L0=\n");L0=trans1();double beta,Bf,Nf,Z;beta=x/6367558.4969;Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(beta)*co s(beta))*1e-10*cos(beta)*sin(beta);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*cos(B f);Z=y/(Nf*cos(Bf));double b2,b3,b4,b5;b2=(0.5+0.003369*cos(Bf)*cos(Bf))*cos(Bf)*sin(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);double B,L,l,R;B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;l=(1-(b3-b5*Z*Z)*Z*Z)*Z;L=L0+l;R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l*l)*l*l*cos(B)*cos(B))* l*sin(B);printf("大地纬度B=\n");trans2(B);printf("大地经度L=\n");trans2(L);printf("平面子午收敛角R=\n");trans2(R);}void Gauss3() //换带计算，克拉索夫斯基椭球{double x1,y1,L10;printf("请输入平面坐标系中的原坐标x1=\n");scanf("%lf",&x1);printf("请输入平面坐标系中的原坐标y1=\n");scanf("%lf",&y1);printf("请以度分秒的形式输入（如10°10′10″）原中央子午线L10=\n");L10=trans1();double beta,Bf,Nf,Z;beta=x1/6367558.4969;Bf=beta+(50221746+(293622+(2350+22*cos(beta)*cos(beta))*cos(beta)*cos(beta))*cos(beta)*co s(beta))*1e-10*cos(beta)*sin(beta);Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf))*cos(Bf)*cos(B f);Z=y1/(Nf*cos(Bf));double b2,b3,b4,b5;b2=(0.5+0.003369*cos(Bf)*cos(Bf))*cos(Bf)*sin(Bf);b3=0.333333-(0.166667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);b5=0.2-(0.1667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);double B,L,l1;B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2;l1=(1-(b3-b5*Z*Z)*Z*Z)*Z;L=L10+l1;double L20,l2,R;printf("请以度分秒的形式（如10°10′10″）输入新中央子午线L20=\n");L20=trans1();l2=L-L20; //R=(1+((0.33333+0.00674*cos(B)*cos(B))+(0.2*cos(B)*cos(B)-0.0067)*l2*l2)*l2*l2*cos(B)*cos (B))*l2*sin(B);//子午线收敛角γ（仅与B和l有关）double N,a0,a4,a6,a3,a5; //求解一些常数以及系数N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.1666667;a5=0.0083-(0.1667-(0.1968+0.004*cos(B)*cos(B))*cos(B)*cos(B))*cos(B)*cos(B);double x2,y2;x2=6367558.4969*B-(a0-(0.5+(a4+a6*l2*l2)*l2*l2)*l2*l2*N)*sin(B)*cos(B);y2=(1+(a3+a5*l2*l2)*l2*l2)*l2*N*cos(B);printf("平面坐标系中的坐标x2=%.4fM\n",x2);printf("平面坐标系中的坐标y2=%.4fM\n",y2);printf("平面子午线收敛角R=\n");trans2(R);printf("大地纬度B=\n");trans2(B);printf("大地经度L=\n");trans2(L);}void main(){int k;printf("请输入k=\n");printf("注：k=1表示正算（BL→xy）k=2表示反算（xy→BL）k=3表示换带（x1y1→x2y2）\n");scanf("%d",&k);if(k==1)Gauss1();if(k==2)Gauss2();if(k==3)Gauss3();}。
基于VB6.0的高斯投影坐标转换的实现
收 稿 日期 ：０ ０— ２ ２ ２１Leabharlann １ — ２作 者简 介 ： 安
卫 （９４一 ）男 ， １８ ， 彝族 ， 州大 方人 ， 贵 助理工 程 师 ， 士 ， 学 主要从 事城镇 地籍 测量 及 Ｃ Ｄ与 ＧＳ 面 的二次 开发 工作 。 Ａ Ｉ方
我国规定按经差 ６ 和 ３进行分 带 ， 。 。 一般 情况下 ， 在进行
ｎ—ｂ
（） ９
采用式 （ ２）～式 （ ９）即 － 计 算 出 高 斯 投 影 的 正 算 ． ｊ
大 比例尺测 图和工程测量时 采用 ３带投影 。在特殊 的情况 。 下， 工程测量控制 网 也可 以采 用 １５ 带或 任 意带 。但 是 为 ．。
Ｏ 引 言
在 日常 的 测 量 工 作 中 ， 业 测 量 的 坐 标 一 般 为 外 ＷＧ Ｓ一８ ４坐标 ， 也称 经纬 度坐 标 ， 而最终 的测 图 成果 均为 平面 直 角坐标 系 ， 因此 这 就 涉 及 经 纬 度 坐标 与平 面 直 角
㈩
＝ Ｙ （ 曰） ， Ｊ
Ｐ ｏｅｔ ｎ Ｕｓ ｇＶＢ ６ ０ ｒｊｃｉ ｉ ． ｏ ｎ
ＡＮ Ｗ ｅ 。 ＧＥ Ｙａ ｇ ， ＡＯ Ｗ ｅ ， ＯＮＧ ｉ ｎ Ｃ ｉ Ｓ Ｂｏ
，
（ ． ｉｎｉ ｎｔ ｕｅｏ ｕ ｖｙｎ ｎ ａ ｐｎ ， ｉｎｉ ０ ３ １ Ｃ ｉａ １ Ｔａｊ Ｉｓｉ ｔ ｆ ｒｅｉｇａ ｄＭ ｐ ｉｇ Ｔａ ｊ ３０ ８ ， ｈｎ ； ｎ ｔ Ｓ ｎ ２ Ｔ ｅＸ ｑｎ ｒｎｈ ｏ ｉｎｉ ｎｃ ａ Ｂ ｒａ ｆ ａ ｄ Ｒ ｓｕ ｃｓａ ｄＨｏ ｎ ａ ａ ｅ ｎ ， ｉｎｉ ０ ３ １ Ｃ ｎ ） ． ｈ ｉｉｇＢ ａ ｃ ｆ ａ ｊ Ｍｕ ｉｉ ｌ ｕｅ ｕｏ ｎ ｅｏ ｒｅ ｍｉｇＭ ｎ ｇ ｍｅ ｔＴａ ｊ ３０ ８ ， ｈ ａ Ｔ ｎ ｐ Ｌ ｎ ｎ ｉ
测绘程序设计(VS2008)实验报告--高斯投影正反算
《测绘程序设计（）》上机实验报告（Visual C++.Net）班级：学号：姓名：序号：二零一一年五月实验7 常用测量程序设计1.实验目的：1.1 巩固类的创建与使用;1.2 掌握数组参数的传递；1.3 掌握常用测绘程序设计的技巧。
1.42.实验内容：编写高斯投影正反算程序。
3.设计思路：这次的实验目的是实现高斯正反算。
需要考虑投影方式即分带的方式，又要考虑椭球参数的类型，所以我添加了两个函数来完成此功能。
分别是int SetProjectType(int m)和void SetParameter(int m,double &a,double &b)。
4.界面设计：界面设计很简单，具体见运行结果。
5.主要代码：文件名：GaussProjectDlg.cpp代码：const double PI=4*atan(1.0);//获得分带方式返回中央子午线经度int CGaussProjectDlg::SetProjectType(int m){UpdateData(TRUE);int n; //记录分带带号double L; //经度L=iDegreeL+iMinL/60+dSecondL/3600;if(m==1) //6度带{n=int(L/6)+1;L0=6*n-3;}else if(m==2) //3度带{n=int((L+1.5)/3);L0=3*n;}else if(m==3) //自主分带L0=L0;return L0;}//获取椭球参数void CGaussProjectDlg::SetParameter(int m,double &a,double &b) {if(m==1) //克拉索夫斯基椭球{a=6378245.0;b=6356863.0187730473;//e=sqrt(0.006693421622966);}else if(m==2) //1975国际协议椭球{a=6378140.0;b=6356755.2881575287;//e=sqrt(0.006694384999588);}else if(m==3) //WGS-84椭球{a=6378137.0;b=6356752.3142;//e=sqrt(0.0066943799013);}}void CGaussProjectDlg::OnBnClickedButtonpositivecal(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double N;double t;double Eta;double X;double A0,A2,A4,A6,A8;double RadB;double Rou;Rou=180*3600/PI;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double l;L0=SetProjectType(iProjectType);double L;L=iDegreeL+double(iMinL)/60+dSecondL/3600;l=(L-L0)*3600;RadB=(iDegreeB+double(iMinB)/60+dSecondB/3600)*PI/180;N=a/sqrt(1-e1*e1*sin(RadB)*sin(RadB));t=tan(RadB);Eta=e2*cos(RadB);A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);A2=-1.0/2*(3.0/4*e1*e1+60.0/64*pow(e1,4)+525.0/512*pow(e1,6)+17640.0/16384*pow(e1,8));A4=1.0/4*(15.0/64*pow(e1,4)+210.0/512*pow(e1,6)+8820.0/16384*pow(e1,8));A6=-1.0/6*(35.0/512*pow(e1,6)+2520.0/16384*pow(e1,8));A8=1.0/8*(315.0/16384*pow(e1,8));X=a*(1-e1*e1)*(A0*RadB+A2*sin(2*RadB)+A4*sin(4*RadB)+A6*sin(6*RadB)+A8*sin(8*RadB));x=X+N/(2*Rou*Rou)*sin(RadB)*cos(RadB)*l*l+N/(24*pow(Rou,4))*sin(RadB)*pow(cos(RadB),3)*(5-t*t+9*Eta*Eta+4*pow(Eta,4))*pow(l,4)+ N/(720*pow(Rou,6))*sin(RadB)*pow(cos(RadB),5)*(61-58*t*t+pow(t,4))*pow(l,6);y=N/Rou*cos(RadB)*l+N/(6*pow(Rou,3))*pow(cos(RadB),3)*(1-t*t+Eta*Eta)*pow(l,3)+N/(120*pow(Rou,5))*pow(cos(RadB),5)*(5-18*t*t+pow(t,4)+14*Eta*Eta-58*Eta*Eta*t*t)*pow( l,5);UpdateData(FALSE);}void CGaussProjectDlg::OnBnClickedButtonantical(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double t_f;double Eta_f;double B_f;double N_f;double M_f;double X=x;double B0;double K0,K2,K4,K6;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double A0;A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);B0=X/(a*(1-e1*e1)*A0);K0=1.0/2*(3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8));K2=-1.0/3*(63.0/64*pow(e1,4)+1108.0/512*pow(e1,6)+58239.0/16384*pow(e1,8));K4=1.0/3*(604.0/512*pow(e1,6)+68484.0/16384*pow(e1,8));K6=-1.0/3*(26328.0/16384*pow(e1,8));B_f=B0+sin(2*B0)*(K0+sin(B0)*sin(B0)*(K4+K6*sin(B0)*sin(B0)));t_f=tan(B_f);Eta_f=e2*cos(B_f);N_f=a/sqrt(1-e1*e1*sin(B_f)*sin(B_f));M_f=N_f/(1+e2*e2*cos(B_f)*cos(B_f));double B;B=B_f-t_f/(2*M_f*N_f)*y*y+t_f/(24*M_f*pow(N_f,3))*(5+3*t_f*t_f+Eta_f*Eta_f-9*Eta_f*Eta_f*t_f*t_f)*pow(y,4)- t_f/(720*M_f*pow(N_f,5))*(61+90*t_f*t_f+45*pow(t_f,4))*pow(y,6);double l;l=1.0/(N_f*cos(B_f))*y-1.0/(6*pow(N_f,3)*cos(B_f))*(1+2*t_f*t_f+Eta_f*Eta_f)*pow(y,3)+1.0/(120*pow(N_f,5)*cos(B_f))*(5+28*t_f*t_f+24*pow(t_f,4)+6*Eta_f*Eta_f+8*Eta_f*Eta_f* t_f*t_f)*pow(y,5);//将B转化为度分秒的形式double dDegB;dDegB=B*180/PI;iDegreeB=int(dDegB);iMinB=int((dDegB-iDegreeB)*60);dSecondB=((dDegB-iDegreeB)*60-iMinB)*60;double dDegL;dDegL=l*180/PI+L0;iDegreeL=int(dDegL);iMinL=int((dDegL-iDegreeL)*60);dSecondL=((dDegL-iDegreeL)*60-iMinL)*60;UpdateData(FALSE);}6.运行结果：实验的运行结果如下图所示：正算：反算：7.实验总结这次实验是实现高斯投影的正反算，方法很多，实现并不复杂，但是计算公式复杂，变量繁多，稍有不慎，就会造成计算错误。
高斯投影正反算程序
printf("高斯投影反算\nunder CGCS2000,六度带\n测试数据：x=2562083.2708,y=19512837.2851\n");
double a=6378137;
double f=1/298.257222101;
double e2=-f*f+2*f;
double x=2562083.2708;
double y=19512837.2851;
int Y=int(y);
int count=0;
while (Y>=100)
{
Y=Y/10;
count++;
}
y=y-Y*pow(10.0,count)-500000;
{
int model;
printf("高斯投影正算请输入1，反算请输入2\n");
scanf("%d",&model);
if(model==1)
{
printf("高斯投影正算\nunder WGS-84,六度带\n测试数据：B=22.155898294，L=111.285215387\n");
double W=sqrt(1-e2*sin(B)*sin(B));
double N=a/W;
double ep2=e2/(1-e2);
double ita2=ep2*cos(B)*cos(B);
double t=tan(B);
double temp1=5-t*t+9*ita2+4*ita2*ita2;
fBf=-a2*sin(2*Bf)/2+a4*sin(4*Bf)/4-a6*sin(6*Bf)/6+a8*sin(8*Bf)/8;
高斯投影VB源码
高斯投影VB源码坐标转换代码double a, f, e2, e12; // 基本椭球参数double A1, A2, A3, A4; // 用于计算X的椭球参数Const PI = 3.14159265358979Const Y0 = 500000#Dim A1#, A2#, A3#, A4#, a#, f#'高斯正算求XPublic Function X(B#, L#, L0#, zbx%) As Double // B纬度L经度L0中央子午线zbcs zbxe2# = 1 - ((f - 1) / f) * ((f - 1) / f)e12# = (f / (f - 1)) * (f / (f - 1)) - 1B = D_R(B): L = D_R(L): L0 = D_R(L0)Dim n, t, t2, m, m2, ng2Dim sinB, cosBX = A1 * B * 180# / PI + A2 * Sin(2 * B) + A3 * Sin(4 * B) + A4 * Sin(6 * B)sinB = Sin(B)cosB = Cos(B)t = Tan(B)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)m = cosB * (L - L0)m2 = m * mng2 = cosB * cosB * e2 / (1 - e2)X = X + n * t * ((0.5 + ((5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24# + (61 - 58 * t2 + t2 * t2) * m2 / 720#) * m2) * m2)End Function'高斯正算求YPublic Function Y(B#, L#, L0#, zbx%) As Doublezbcs zbxe2# = 1 - ((f - 1) / f) * ((f - 1) / f)e12# = (f / (f - 1)) * (f / (f - 1)) - 1B = D_R(B): L = D_R(L): L0 = D_R(L0)Dim n, t, t2, m, m2, ng2Dim sinB, cosBsinB = Sin(B)cosB = Cos(B)t = Tan(B)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)m = cosB * (L - L0)m2 = m * mng2 = cosB * cosB * e2 / (1 - e2)Y = n * m * (1 + m2 * ((1 - t2 + ng2) / 6# + m2 * (5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2) / 120#))Y = Y + Y0End Function'高斯反算求B（纬度）Public Function B(X#, Y#, L0#, zbx%) As Doublezbcs zbxe2# = 1 - ((f - 1) / f) * ((f - 1) / f)e12# = (f / (f - 1)) * (f / (f - 1)) - 1L0 = D_R(L0)Dim sinB, cosB, t, t2, n, ng2, V, yNDim preB0, B0Dim etaY = Y - Y0B0 = X / A1Do While TruepreB0 = B0B0 = B0 * PI / 180#B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0))) / A1eta = Abs(B0 - preB0)If eta < 0.0000001 Then Exit DoLoopB0 = B0 * PI / 180#sinB = Sin(B0)cosB = Cos(B0)t = Tan(B0)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)ng2 = cosB * cosB * e2 / (1 - e2)V = Sqr(1 + ng2)yN = Y / nB = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12# + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360#) * V * V * t / 2B = R_D(B)End Function'高斯反算求L（经度）Public Function L(X#, Y#, L0#, zbx%) As Doublezbcs zbxe2# = 1 - ((f - 1) / f) * ((f - 1) / f)e12# = (f / (f - 1)) * (f / (f - 1)) - 1L0 = D_R(L0)Dim sinB, cosB, t, t2, n, ng2, V, yNDim preB0, B0Dim etaY = Y - Y0B0 = X / A1Do While TruepreB0 = B0B0 = B0 * PI / 180#B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0))) / A1eta = Abs(B0 - preB0)If eta < 0.0000001 Then Exit DoLoopB0 = B0 * PI / 180#sinB = Sin(B0)cosB = Cos(B0)t = Tan(B0)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)ng2 = cosB * cosB * e2 / (1 - e2)V = Sqr(1 + ng2)yN = Y / nL = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6# + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120#) / cosBL = R_D(L)End Function'把弧度转化为DD.FFMM格式Public Function R_D(R) As DoubleR = R * 180 / PIDD = Int(R)FF = Int((R - DD) * 60#)MM = ((R - DD) * 60# - FF) * 60#If Abs(MM - 60) < 0.0001 ThenMM = 0: FF = FF + 1End IfIf FF = 60 ThenFF = 0: DD = DD + 1End IfR_D = DD + FF / 100 + MM / 10000End Function'把DD.FFMM转化为弧度Public Function D_R(DFM) As Double'DFM = Round(DFM, 4)Dim DD%, FF%, MM#DD = Int(DFM)FF = Int((DFM - Int(DFM)) * 100 + 0.001)MM = ((DFM * 100) - Int(Round(DFM * 100, 3))) * 100D_R = (DD + FF / 60 + MM / 3600) * PI / 180End FunctionSub zbcs(DaiHao As Integer)Select Case DaiHaoCase 1A1 = 111134.8611A2 = -16036.4803A3 = 16.8281A4 = -0.022a = 6378245#f = 298.3Case 2A1 = 111133.0047A2 = -16038.5282A3 = 16.8326A4 = -0.022a = 6378140 1975 year，IUUG，长轴6378140f = 298.257(扁率)Case ElseA1 = 111134.8611A2 = -16036.4803A3 = 16.8281A4 = -0.022a = 6378245#f = 298.3End SelectEnd Sub。
高斯投影正反算及换带计算VB程序设计
摘要本设计主要阐述了高斯投影分带以及高斯投影坐标正、反算的推导公式，从而根据公式来编写基于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程序的操作介绍下面以实例来介绍程序的操作步骤。
高斯正反算VB实验代码
高斯正反算实验报告一、实验目的：通过对高斯正反算的软件编程，增强学生对《大地测量学基础》课程的理解，使学生牢固掌握高斯正反算的基本原理和公式，熟悉高斯正反算的基本技能和计算方法。
培养学生利用计算机编程的能力、综合运用知识的能力及理论联系实际的能力二、实验要求：要求学生综合运用测绘知识、大地测量知识、数学知识和计算机知识，设计数学模型和程序算法，编制程序实现数据的自动化处理。
三、实验环境Visual studio 2008界面如下：输入数据计算:代码如下：Public Class高斯投影正反算Private Sub高斯投影正反算_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load'自动化控件显示初始值RadioButton2.Checked = TrueRadioButton3.Checked = TrueTextBox1.Focus()TextBox1.Text = "0"TextBox2.Text = "0"TextBox3.Text = "0"TextBox4.Text = "0"TextBox5.Text = "0"TextBox6.Text = "0"RichTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标"End SubPrivate Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click'获取输入数据Dim bb1 As Double, bb2 As Double, bb3 As Double, ll1 As Double, ll2 As Double, ll3 As Doublebb1 = TextBox1.Textbb2 = TextBox2.Textbb3 = TextBox3.Textll1 = TextBox4.Textll2 = TextBox5.Textll3 = TextBox6.TextDim B As Double, L As Double, x As Double, y As Double, X1 As Double, Y1 As Double'X1 Y1Dim N As IntegerDim L0 As Double, l1 As DoubleDim p As Doublep = 206264.80625'检查输入格式的正确性If (bb1 >= 0 And bb1 < 90 And bb2 >= 0 And bb2 < 60 And bb3 >= 0 And bb3 < 60) Then B = bb1 * 3600 + bb2 * 60 + bb3ElseRichTextBox1.Text = "纬度输入格式不正确！"ReturnEnd IfIf (ll1 >= 0 And ll1 < 360 And ll2 >= 0 And ll2 < 60 And ll3 >= 0 And ll3 < 60) Then L = ll1 * 3600 + ll2 * 60 + ll3ElseRichTextBox1.Text = "经度输入格式不正确！"ReturnEnd IfDim b1 As Double'b1b1 = B / p'获取用户选项值Dim tyd As IntegerIf (RadioButton2.Checked = True) Thentyd = 6Elsetyd = 3End IfDim ty As IntegerIf (RadioButton3.Checked = True) Thenty = 1Elsety = 2End If'按带求带号以及中央子午线经度()If (tyd = 6) ThenN = Int(L / (6 * 3600.0)) + 1L0 = (6 * N - 3) * 3600.0l1 = (L - L0) / pElseIf (tyd = 3) ThenN = Int(L / (3 * 3600.0) + 0.5)L0 = 3 * N * 3600.0l1 = (L - L0) / pEnd If'采用克氏椭球参数的计算公式If (ty = 1) ThenDim c As DoubleDim c1 As DoubleDim c2 As DoubleDim l2 As DoubleDim n1 As DoubleDim a0 As DoubleDim a3 As DoubleDim a4 As DoubleDim a5 As DoubleDim a6 As DoubleDim i As Integeri = 1c = Math.Pow(Math.Cos(b1), 2)c1 = Math.Sin(b1) * Math.Cos(b1)c2 = Math.Cos(b1)l2 = Math.Pow(l1, 2)n1 = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c 'n1a0 = 32140.404 - (135.3302 - (0.7092 - 0.004 * c) * c) * ca4 = (0.25 + 0.00252 * c) * c - 0.04166a6 = (0.166 * c - 0.084) * ca3 = (0.3333333 + 0.001123 * c) * c - 0.1666667a5 = 0.0083 - (0.1667 - (0.1968 + 0.004 * c) * c) * cx = 6367558.4969 * b1 - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n1) * c1y = (1 + (a3 + a5 * l2) * l2) * l1 * n1 * c2X1 = xDim y11 As Doubley11 = y + 500000.0If y11 / i > 1 ThenY1 = N * i * 10 + y11 'Y坐标加代号i = i * 10End IfDim tuoqiu As Stringtuoqiu = "采用克氏椭球参数,"RichTextBox1.Text = "按经差" + tyd.ToString("0") + "°进行投影分带，" + tuoqiu + "其计算结果为：" + vbCrLf + "x=" + x.ToString("0.00000000") + vbCrLf + "y=" +y.ToString("0.00000000") + vbCrLf + "国家统一坐标为："+ vbCrLf + "X="+ X1.ToString("0.00000000") + vbCrLf + "Y=" + N.ToString + Y1.ToString("0.00000000")End If'采用1975国际椭球参数的计算公式()If (ty = 2) ThenDim c As DoubleDim c1 As DoubleDim c2 As DoubleDim l2 As DoubleDim n1 As DoubleDim a0 As DoubleDim a3 As DoubleDim a4 As DoubleDim a5 As DoubleDim a6 As DoubleDim i As Integeri = 1c = Math.Cos(b1) * Math.Cos(b1)c1 = Math.Sin(b1) * Math.Cos(b1)c2 = Math.Cos(b1)l2 = Math.Pow(l1, 2)n1 = 6399596.652 - (21565.045 - (108.996 - 0.603 * c) * c) * ca0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * c) * c) * ca4 = (0.25 + 0.00253 * c) * c - 0.04167a6 = (0.167 * c - 0.083) * ca3 = (0.3333333 + 0.001123 * c) * c - 0.1666667a5 = 0.00878 - (0.1702 - 0.20382 * c) * cx = 6367452.1328 * b1 - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n1) * c1y = (1 + (a3 + a5 * l2) * l2) * l1 * n1 * c2X1 = xDim y11 As Doubley11 = y + 500000.0If y11 / i > 1 ThenY1 = N * i * 10 + y11 'Y坐标加代号i = i * 10End IfDim tuoqiu As Stringtuoqiu = "采用1975国际椭球参数,"RichTextBox1.Text = "按经差" + tyd.ToString("0") + "°进行投影分带，" + tuoqiu + "其计算结果为：" + vbCrLf + "x=" + x.ToString("0.00000000") + vbCrLf + "y=" +y.ToString("0.00000000") + vbCrLf + "国家统一坐标为："+ vbCrLf + "X="+ X1.ToString("0.00000000") + vbCrLf + "Y=" + N.ToString + Y1.ToString("0.00000000")End IfEnd SubPrivate Sub Button2_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button2.ClickDim y As Double, X1 As Double, Y1 As Double, B As Double, L As DoubleX1 = TextBox7.Text '获取输入数据Y1 = TextBox8.TextDim L0 As DoubleDim i As Integeri = 1Do While Y1 / i >= 10y = Y1 - Int(Y1 / i) * i - 500000L0 = 6 * Int(Y1 / i) - 3i = i * 10Loop'对Y坐标处理并求出中央子午线经度'按6°带克氏椭球反算Dim bt As Double, BT1 As Double, c3 As Double, c4 As Double, Bf As Double, c5 As Double, c6 As Double, Nf As Double, Z As DoubleDim b2 As Double, b3 As Double, b4 As Double, b5 As Double, z2 As DoubleDim p As Doublep = 206264.80625bt = X1 / 6367558.4969 * pBT1 = X1 / 6367558.4969c3 = Math.Cos(BT1) * Math.Cos(BT1)c4 = Math.Sin(BT1) * Math.Cos(BT1)Bf = (bt + (50221746 + (293622 + (2350 + 22 * c3) * c3) * c3) * c4 * Math.Pow(10, -10) * p) / pc5 = Math.Pow(Math.Cos(Bf), 2)c6 = Math.Sin(Bf) * Math.Cos(Bf)Nf = 6399698.902 - (21562.267 - (108.973 - 0.612 * c5) * c5) * c5Z = y / (Nf * Math.Cos(Bf))b2 = (0.5 + 0.003369 * c5) * c6b3 = 0.333333 - (0.166667 - 0.001123 * c5) * c5b4 = 0.25 + (0.16161 + 0.00562 * c5) * c5b5 = 0.2 - (0.1667 - 0.0088 * c5) * c5z2 = Math.Pow(Z, 2)B = (Bf * p - (1 - (b4 - 0.12 * z2) * z2) * z2 * b2 * p) / 3600.0L = L0 + ((1 - (b3 - b5 * z2) * z2) * Z * p) / 3600.0RichTextBox2.Text = "按6°带克氏椭球反算后，结果为：" + vbCrLf + "B=" +B.ToString("0.0000000000000") + vbCrLf + "L=" + L.ToString("0.0000000000000")End SubPrivate Sub Button3_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button3.ClickTextBox1.Focus()TextBox1.Text = "0"TextBox2.Text = "0"TextBox3.Text = "0"TextBox4.Text = "0"TextBox5.Text = "0"TextBox6.Text = "0"RichTextBox1.Text = " "End SubPrivate Sub Button4_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button4.ClickTextBox7.Text = " "TextBox8.Text = " "RichTextBox2.Text = "说明：输入的坐标需为按6°带投影且采用克氏椭球参数所得的国家统一坐标"End SubPrivate Sub Button5_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button5.ClickEndEnd SubEnd Class。
基于Visual C++ 的高斯投影坐标反算
术 Applied Technology 应 用 技
基于 Visual C++的高斯投影坐标反算
许大亮
渊安徽理工大学测绘学院袁 安徽 淮南 232001冤
摘 要院文中利用 Visual C+ + 语言编写高斯投影坐标反算的程序,采用的数据是西安 80 椭球参数和高斯坐标曰 首先 根据子午线弧长公式迭代计算求出垂足纬度袁 其次算出大地纬度和经差袁 最后可以得到大地坐标遥 程序计算出的坐 标精度为 0.000 1义遥 结果表明院 用转换程序解算出的坐标存在一定的误差袁 但是可以满足低精度的工程建设遥 关键词院Visual C++语言曰坐标反算曰高斯坐标曰垂足纬度 中图分类号院TP311 文献标志码院A DOI院10.3969/j.issn.1674-9146.2015.08.059
根据式 渊3冤 可写出迭代式
Bi+1 = X/c 茁0 -
蓸 蔀 3
5
7
茁2cosBi + 茁4cos Bi + 茁6cos Bi + 茁8cos Bi
sinBi 茁0
.
渊4冤
1.3 Mf袁 Nf袁 tf袁 浊f等关于Bf 的公式
地球椭球的长半径为 a袁 短半径为 b遥
第一偏心率
e
=
姨a2
-
2
b
a.
换算后坐标的精度遥 在垂足纬度的计算中袁 常采用
迭代计算的方法遥 笔者参考了一些文献尧 著作袁 利
用 Visual C++语言设计了一个程序袁 实现了高斯投
影坐标转换成西安 80 大地坐标袁 计算结果精确可
VB报告(坐标正反算使用说明)
一、输入界面1.单击VB中的“运行”快捷键，弹出图1所示的运行界面图12.在图1中选择按钮，进入坐标正算模式，如图2图2在图2中输入已知点坐标、已知点至未知点的边长和坐标方位角、保留的小数位数，点名可不输入，输入完后，选择“计算”按钮。
如要计算多个点坐标，则计算完一个点后用鼠标单击“刷新”按钮，重复以上操作即可。
若想退出坐标正算功能，用鼠标左键单击“退出”按钮，在弹出的提示框中选择“是”，如图3图33. 在图1中选择按钮，进入坐标反算模式，如图4图4在图4中输入两个已知点坐标、保留的小数位数，点名可不输入，输入完后，选择“计算”按钮。
如要计算多个点坐标，则计算完一个点后用鼠标单击“刷新”按钮，重复以上操作即可。
若想退出坐标反算功能，用鼠标左键单击“退出”按钮，在弹出的提示框中选择“是”，如图3二、输入，输出数据表格中的方位角以小数的形式输入，例如“321°18′56″“的输入格式为321.1856，若度数为321°18′56.5″的输入格式为321.18565，坐标的保留位数保留三位，别的以此类推方位角以小数的形式输入，例如“164°02′02″”的输入格式为164.0202，若度数为164°02′02.5″的输入格式为164.02025，表格中的方位角的小数保留位数为0位三、源代码Private Sub Command1_Click()If Option1.Value = True And Option2.Value = False Thenzhengsuan '当选择坐标正算按钮时调用坐标正算程序End IfIf Option1.Value = False And Option2.Value = True Thenfansuan '当选择坐标反算按钮时调用坐标反算程序End IfEnd Sub'坐标正算程序Private Sub zhengsuan()Dim s As DoubleDim a As DoubleDim x As DoubleDim y As DoubleDim m As DoubleDim n As DoubleDim degree As DoubleDim minute As DoubleDim second As DoubleDim rad As DoubleDim bbt As IntegerDim result As DoubleDim p As Doubles = Val(Text5.Text)a = Val(Text6.Text)Rad_do a, degree, minute, second, rad '调用将度分秒转化为弧度的程序x = Val(Text2.Text) + s * Cos(rad)y = Val(Text3.Text) + s * Sin(rad)p = Fix(x)x = x - pbbt = Val(Text4.Text)Sheru m, n, bbt, x, result '调用奇进偶舍程序Label12.Caption = result + pp = Fix(y)y = y - pSheru m, n, bbt, y, result '调用奇进偶舍程序Label13.Caption = result + pEnd Sub'坐标反算程序Private Sub fansuan()Dim s As DoubleDim x As DoubleDim y As DoubleDim sing As DoubleDim m As DoubleDim n As DoubleDim f As DoubleDim result As DoubleDim bbt As IntegerDim degree As DoubleDim minute As DoubleDim second As DoubleDim rad As DoubleDim p As Doublex = Val(Text5.Text)y = Val(Text6.Text)m = y - Val(Text3.Text)n = x - Val(Text2.Text)If n <> 0 Thenm = m / nrad = Atn(m)sing = Sgn(rad)End Ifdegree_m_s degree, minute, second, rad '调用将弧度转化为度分秒的程序bbt = Val(Text8.Text)m = second * 10 ^ bbtp = Fix(second)second = second - pSheru m, n, bbt, second, result '调用奇进偶舍程序second = result + pfangweijiao degree, minute, second, (x - Val(Text2.Text)), (y - Val(Text3.Text)), sing '调用计算坐标方位角的程序zhuanhua minute, second, degreesecond = second * 10 ^ bbtLabel13.Caption = degree & "." & minute & secondIf minute < 10 ThenLabel13.Caption = degree & "." & "0" & minute & secondEnd IfIf second < 10 ^ (bbt + 1) ThenLabel13.Caption = degree & "." & minute & "0" & secondEnd IfIf minute < 10 And second < 10 ^ (bbt + 1) ThenLabel13.Caption = degree & "." & "0" & minute & "0" & secondEnd If'计算边长Ss = (y - Val(Text3.Text)) ^ 2 + (x - Val(Text2.Text)) ^ 2s = Sqr(s)bbt = Val(Text4.Text)p = Fix(s)s = s - pSheru m, n, bbt, s, result '调用奇进偶舍程序Label12.Caption = result + pEnd Sub'将度分秒化为弧度Private Sub Rad_do(ByVal a As Double, ByVal degree As Double, ByVal minute As Double, ByVal second As Double, ByRef rad As Double)degree = a \ 1a = a - Fix(a)minute = Fix(a * 100)second = (a * 100 - minute) * 100rad = (3600 * degree + 60 * minute + second) / 206264.8063End Sub'将弧度转化为度分秒的程序Private Sub degree_m_s(ByRef degree As Double, ByRef minute As Double, ByRef second As Double, ByRef rad As Double)rad = rad * 180 / 3.1415926535degree = rad \ 1rad = (rad - degree) * 60minute = rad \ 1second = (rad - minute) * 60End Sub'奇进偶舍程序Private Sub Sheru(ByVal m As Double, ByVal n As Double, ByVal bbt As Double, ByVal x As Double, ByRef result As Double)m = x * 10 ^ bbtn = m - m \ 1If n < 0.5 Thenresult = (m \ 1) / 10 ^ bbtElseIf n > 0.5 Thenresult = (m \ 1 + 1) / 10 ^ bbtElseIf (m \ 1) Mod 2 Thenresult = (m \ 1 + 1) / 10 ^ bbtElseresult = (m \ 1) / 10 ^ bbtEnd IfEnd IfEnd Sub'计算坐标方位角的程序Private Sub fangweijiao(ByRef degree As Double, ByRef minute As Double, ByRef second As Double, ByVal m, ByVal n As Double, ByVal sing As Integer)Dim i As IntegerIf sing = 1 ThenIf n < 0 And m < 0 Thendegree = degree + 180End IfElseIf sing = -1 ThenIf n > 0 And m < 0 Thendegree = degree + 179minute = minute + 59second = second + 60ElseIf n < 0 And m > 0 Thendegree = degree + 359minute = minute + 59second = second + 60End IfEnd Ifi = Sgn(n)If m <> 0 ThenIf m < 0 And n = 0 Thendegree = "180"ElseIf m > 0 And n = 0 Thendegree = "0"End IfElseIf i = 1 Thendegree = "90"ElseIf i = -1 Thendegree = "270"ElseMsgBox "您输入了两个相同的点，请重新输入！"End IfEnd IfEnd Sub'当分秒超过60时须向上一级进位及方位角度数超过360°须减360°的程序Private Sub zhuanhua(ByRef minute As Double, ByRef second As Double, ByRef degree As Double)If second >= 60 Thenminute = minute + 1second = second - 60ElseIf second < 0 Thenminute = minute - 1second = second + 60End IfIf minute >= 60 Thendegree = degree + 1minute = minute - 60ElseIf minute < 0 Thendegree = degree - 1minute = minute + 60End IfIf degree >= 360 Thendegree = degree - 360End IfEnd Sub'退出程序Private Sub Command3_Click()If MsgBox("是否退出？", vbYesNo, "提示") = vbYes ThenUnload MeEnd IfEnd Sub'刷新程序Private Sub Command4_Click()Text2.Text = ""Text3.Text = ""Text4.Text = 3Text5.Text = ""Text6.Text = ""Text8.Text = 1Label12.Caption = ""Label13.Caption = ""End Sub'设置窗体的大小，使窗体充满整个屏幕并对label6和label7赋值Private Sub Form_Load()Me.Height = Screen.HeightMe.Width = Screen.WidthMe.Left = 0Me.Top = 0Label6.Caption = "边长（S）"Label7.Caption = "方位角（a）"Label17.Caption = "陈亮编程"Label18.Caption = " 2011.09.05"End Sub'设置坐标正算时Label6.Caption和Label7.Caption的值及设计小数点位数Private Sub Option1_Click()If Option1.Value = True And Option2.Value = False ThenLabel14.Caption = Text1.Text & "—>未知点" & Text7.TextLabel6.Caption = "边长（S）"Label7.Caption = "方位角（a）"Label10.Caption = "X"Label11.Caption = "Y"Text4.Text = 3Text8.Text = 1End IfEnd Sub''设置坐标反算时Label0.Caption和Label1.Caption的值及设计小数点位数Private Sub Option2_Click()If Option1.Value = False And Option2.Value = True ThenLabel6.Caption = "X"Label7.Caption = "Y"Label14.Caption = ""Label10.Caption = "边长（S）"Label11.Caption = "方位角（a）"Text4.Text = 3Text8.Text = 1End IfEnd Sub。
高斯消去法 矩阵运算 坐标正反算等vb程序设计
高斯消去法Private Sub Command1_Click()Dim a(1 To 10, 1 To 11) As DoubleDim x(1 To 10) As DoubleDim Sum As DoubleDim n As IntegerDim i As IntegerDim j As IntegerDim k As Integern = Val(InputBox("输入未知量个数(最多10个)")) For i = 1 To nFor j = 1 To n + 1a(i, j) = Val(InputBox("输入增广矩阵")) Next jNext iFor i = 1 To nFor j = 1 To n + 1Print a(i, j);Next jPrintNext iRem 消元过程For k = 1 To n - 1For i = k + 1 To nFor j = k + 1 To n + 1a(i, j) = a(i, j) - a(i, k) * a(k, j) / a(k, k) Next jNext iNext kRem 回代过程x(n) = a(n, n + 1) / a(n, n)Print "x"; n; "="; x(k);PrintFor k = n - 1 To 1 Step -1Sum = 0For j = k + 1 To nSum = a(k, j) * X(j) + SumNext jx(k) = (a(k, n + 1) - Sum) / a(k, k)Print "x"; k; "="; x(k);PrintNext kEnd Sub'将度.分秒形式化为弧度：输入为度.分秒形式，输出为弧度Public Function DoToHu(ByVal DoFenMiao As Double) As Single Dim du%, fen%, miao%, angle#du = Fix(DoFenMiao)DoFenMiao = (DoFenMiao - du) * 100fen = Fix(DoFenMiao)miao = (DoFenMiao - fen) * 100angle = du + fen / 60 + miao / 3600DoToHu = angle * PI / 180End Function'弧度化为度.分秒的形式：输入弧度值，输出度.分秒（各占两位）Public Function HuToDo(ByVal Hu As Double) As SingleDim du%, fen%, miao%Hu = Hu * 180 / PIdu = Fix(Hu)Hu = (Hu - du) * 60fen = Fix(Hu)Hu = (Hu - fen) * 60miao = Fix(Hu + 0.5)If miao = 60 Thenfen = fen + 1miao = 0End IfIf fen = 60 Thendu = du + 1fen = 0End IfHuToDo = du + fen / 100 + miao / 10000End Function矩阵运算Dim a() As LongDim b() As LongDim c() As LongDim m As LongDim n As LongDim p As LongDim q As LongPrivate Sub Command1_Click()On Error GoTo lab:m = InputBox("请输入A矩阵行数", "提示")n = InputBox("请输入A矩阵列数", "提示")Picture1.ClsPicture3.ClsReDim a(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To na(i, j) = InputBox("请输入矩阵a(" & i & "," & j & ")数值", "提示") Picture1.Print a(i, j);Next jPicture1.PrintNext ilab:End SubPrivate Sub Command10_Click()EndEnd SubPrivate Sub Command2_Click()On Error GoTo lab:p = InputBox("请输入B矩阵行数", "提示")q = InputBox("请输入B矩阵列数", "提示")Picture2.ClsPicture3.ClsReDim b(1 To p, 1 To q) As LongFor i = 1 To pFor j = 1 To qb(i, j) = InputBox("请输入矩阵b(" & i & "," & j & ")数值", "提示") Picture2.Print b(i, j);Next jPicture2.PrintNext ilab:End SubPrivate Sub Command3_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf m <> p Or n <> q ThenMsgBox "请输入行数和列数相同的矩阵才可相加", vbOKOnly, "提示" End IfIf m = p And n = q ThenLabel1.Caption = "+"ReDim c(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To nc(i, j) = a(i, j) + b(i, j)Picture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End SubPrivate Sub Command4_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf m <> p Or n <> q ThenMsgBox "请输入行数和列数相同的矩阵才可相减", vbOKOnly, "提示" End IfIf m = p And n = q ThenLabel1.Caption = "-"ReDim c(1 To m, 1 To n) As LongFor i = 1 To mFor j = 1 To nc(i, j) = a(i, j) - b(i, j)Picture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End SubPrivate Sub Command5_Click()On Error GoTo lab:Picture3.ClsIf m = 0 Or n = 0 Or p = 0 Or q = 0 ThenMsgBox "请先输入矩阵", vbOKOnly, "提示"End IfIf n <> p ThenMsgBox "请输入A矩阵列数和B矩阵行数相等的矩阵再做乘积", vbOKOnly, "提示" End IfIf n = p ThenLabel1.Caption = "x"ReDim c(1 To m, 1 To q) As LongFor i = 1 To mFor j = 1 To qFor k = 1 To nc(i, j) = c(i, j) + a(i, k) * b(k, j)Next kPicture3.Print c(i, j);Next jPicture3.PrintNext iEnd Iflab:End Sub逆矩阵ReDim a(1 To n, 1 To n * 2)ReDim x(1 To n)For i = 1 To nFor j = 1 To n * 2a(i, j) = Val(InputBox("输入增广矩阵第" & i & "行第" & j & "列元素"))Picture1.Print a(i, j);Next jPicture1.PrintNext iFor k = 1 To n - 1 '消元For i = k + 1 To nFor j = k + 1 To n * 2a(i, j) = a(i, j) - a(i, k) * a(k, j) / a(k, k) Next jNext iNext kFor j = n + 1 To 2 * na(n, j) = a(n, j) / a(n, n)Next jFor k = n To 1 Step -1 '反消元For i = k - 1 To 1 Step -1For j = k + 1 To n * 2a(i, j) = (a(i, j) - a(i, k) * a(k, j)) / a(i, i) Next jNext iNext kFor i = 1 To nFor j = n + 1 To 2 * nPicture2.Print a(i, j); Space(3);Next jPicture2.PrintNext ilab:End Sub坐标正反算Private Sub Command1_Click()XX = Val(Text3.Text) - Val(Text1.Text)YY = Val(Text4.Text) - Val(Text2.Text)If XX = 0 And YY = 0 ThenMsgBox ("A与B点不能为同一点，请重新输入") Text1.Text = ""Text2.Text = ""Text3.Text = ""Text4.Text = ""ElseIf YY = 0 And XX > 0 ThenText5.Text = "00"Text6.Text = "00"Text7.Text = "00"ElseIf XX = 0 And YY > 0 ThenText5.Text = 90Text6.Text = "00"Text7.Text = "00"ElseIf XX < 0 And YY = 0 ThenText5.Text = 180Text6.Text = "00"Text7.Text = "00"ElseIf XX = 0 And YY < 0 ThenText5.Text = 270Text6.Text = "00"Text7.Text = "00"Elsek = Abs(YY / XX)R = Atn(k) * 180 / 3.1415926If XX > 0 And YY > 0 Thena = RElseIf XX < 0 And YY > 0 Thena = 180 - RElseIf XX < 0 And YY < 0 Thena = 180 + RElseIf XX > 0 And YY < 0 Then ‘AB在第四象限a = 360 – REnd Ifb = Format(a, "0.00000") ‘取小数后5位Text5.Text = Fix(b) ‘求度数c = Fix((b - Text5.Text) * 60)Text6.Text = c ‘求分d = Round(((b - Text5.Text) * 60 - c) * 60)Text7.Text = d ‘求秒End IfEnd SubMatlabX=[107563.8100,107620.9521,109989.2770,111411.7664,109584.1244,109506.4700];Y=[571684.5200,568999.8520,568164.3993,569885.3537,572397.4802,570099.6000];a=[0,1,0,0,1,1;1,0,1,0,0,1;0,1,0,1,0,1;0,0,1,0,1,1;1,0,0,1,0,1;1,1,1,1,1,0][m,n]=size(a)plot(Y(1),X(1),'r^',Y(2),X(2),'bo',Y(3),X(3),'bo',Y(4),X(4),'bo',Y(5),X(6),'bo',Y(6),X(6),'r^') hold onfor i=1:mfor j =1:nif a(i,j)==1line([Y(i),Y(j)],[X(i),X(j)])hold onendendendtext(Y(1),X(1),'AA')text(Y(2),X(2),'11')text(Y(3),X(3),'22')text(Y(4),X(4),'33')text(Y(5),X(6),'44')text(Y(6),X(5),'BB')grid。
高斯投影正反算编程
高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式（2）高斯反算基本公式以上主要通过大地测量学基础课程得到.这不进行详细的推导.只是列出基本公式指导编程的进行。
二．编程的基本方法和流程图（1）编程的基本方法高斯投影正反算基本上运用了所有的编程基本语句.本文中是利用C++语言进行基本的设计。
高斯正算中对椭球参数和带宽的选择主要运用了选择语句。
而高斯反算中除了选择语句的应用.在利用迭代算法求底点纬度还应用了循环语句。
编程中还应特别注意相关的度分秒和弧度之间的相互转换.这是极其重要的。
（2）相关流程图1）正算2）反算三．编程的相关代码（1）正算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))int i;struct jin{double B;double L;double L0;};struct jin g[100];main(int argc, double *argv[]) {FILE *r=fopen("a.txt","r"); assert(r!=NULL);FILE *w=fopen("b.txt","w"); assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].B,&g[i].L,&g[i].L 0)!=EOF){double a,b;int zuobiao;printf("\n请输入坐标系：北京54=1.西安80=2.WGS84=3：");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=6378245;b=6356863.0187730473;}if(zuobiao==2){a=6378140;b=6356755.2881575287;}if(zuobiao==3){a=6378137;b=6356752.3142;} //选择坐标系//double f=(a-b)/a;double e,e2;e=sqrt(2*f-f*f);e2=sqrt((a/b)*(a/b)-1);//求椭球的第一.第二曲率//double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;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=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;double Bmiao,Lmiao, L0miao;Bmiao=(int)(g[i].B)*3600.0+(int)((g[i].B-(int)(g[i].B) )*100.0)*60.0+(g[i].B*100-(int)(g[i].B*100))*100.0;Lmiao=(int)(g[i].L)*3600.0+(int)((g[i].L-(int)(g[i].L) )*100.0)*60.0+(g[i].L*100-(int)(g[i].L*100))*100.0;L0miao=(int)(g[i].L0)*3600.0+(int)((g[i].L0-(int)(g[i] .L0))*100.0)*60.0+(g[i].L0*100-(int)(g[i].L0*100))*100 .0;double db;db=pi/180.0/3600.0;double B1,L1,l;B1=Bmiao*db;L1= Lmiao*db;l=L1-L0miao*db;//角度转化为弧度//double T=tan(B1)*tan(B1);double n=e2*e2*cos(B1)*cos(B1);double A=l*cos(B1);double X,x,y;X=a0*(B1)-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6 +a8*sin(8*B1)/8;//求弧长//double N=a/sqrt(1-e*e*sin(B1)*sin(B1));int Zonewide;int Zonenumber;printf("\n请输入带宽：3度带或6度带Zonewide=");scanf("%d",&Zonewide);getchar();if(Zonewide==3){Zonenumber=(int)((g[i].L-Zonewide/2)/Zonewide+1);}else if(Zonewide==6){Zonenumber=(int)g[i].L/Zonewide+1;}else{printf("错误");exit(0);}//选择带宽//doubleFE=Zonenumber*1000000+500000;//改写为国家通用坐标//y=FE+N*A+A*A*A*N*(1-T*T+n*n)/6+A*A*A*A*A*N*(5-18*T*T+T *T*T*T+14*n*n-58*n*n*T*T)/120;x=X+tan(B1)*N*A*A/2+tan(B1)*N*A*A*A*A*(5-T*T+9*n*n+4*n *n*n*n)/24+tan(B1)*N*A*A*A*A*A*A*(61-58*T*T+T*T*T*T)/7 20;printf("\n所选坐标系的转换结果：x=%lf y=%lf\n",x,y);fprintf(w,"%lf %lf\n",x,y);//输出结果到文本文件//}fclose(r);fclose(w);system("pause");return 0;}（2）反算# include "stdio.h"# include "stdlib.h"# include "math.h"# include "assert.h"#define pi (4*atan(1.0))double X,Y,B1,B2,B3,F,t;double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8,a1,b1;double BB,LL,Bf;double e,e1;int d,m,s,i,zuobiao;double sort(double,double);struct jin{double x;double y;double L0;};struct jin g[100];//x,y,L0为输入量：x,y坐标和中央子午线经度//main(int argc, double *argv[]){FILE *r=fopen("c.txt","r");assert(r!=NULL);FILE *w=fopen("d.txt","w");assert(r!=NULL);int i=0;while(fscanf(r,"%lf %lf %lf",&g[i].x,&g[i].y,&g[i].L 0)!=EOF)//文件为空.无法打开//{double a1=6378245.0000000000;//克拉索夫斯基椭球参数//double b1=6356863.0187730473;double a75=6378140.0000000000;//1975国际椭球参数//double b75=6356755.2881575287;double a84=6378137.0000000000;//WGS-84系椭球参数//double b84=6356752.3142000000;double M,N;//mouyou圈曲率半径.子午圈曲率半径// double t,n;double A,B,C;double BB,LL,Bf,LL0,BB0;double a,b;printf("\n选择参考椭球：1=克拉索夫斯基椭球.2=1975国际椭球.3=WGS-84系椭球:");scanf("%d",&zuobiao);getchar();if(zuobiao==1){a=a1;b=b1;}if(zuobiao==2){a=a75;b=b75;}if(zuobiao==3){a=a84;b=b84;}//选择参考椭球.求解第一偏心率e,第二偏心率e1// Bf=sort(a,b);//调用求解底点纬度的函数//double q=sqrt(1-e*e*sin(Bf)*sin(Bf));double G=cos(Bf);M=a*(1-e*e)/(q*q*q);N=a/q;double H,I;A=g[i].y/N;H=A*A*A;I=A*A*A*A*A;t=tan(Bf);n=e1*cos(Bf);B=t*t;C=n*n;BB0=Bf-g[i].y*t*A/(2*M)+g[i].y*t*H/(24*M)*(5+3*B+C-9*B *C)-g[i].y*t*I/(720*M)*(61+90*B+45*B*B);LL0=g[i].L0*pi/180.0+A/G-H/(6*G)*(1.0+2*B+C)+I/(120*G)*(5.0+28*B+24*B*B+6*C+8*B*C);//利用公式求解经纬度// int Bdu,Bfen,Ldu,Lfen;double Bmiao,Lmiao;Ldu=int(LL0/pi*180);Lfen=int((LL0/pi*180)*60-Ldu*60);Lmiao=LL0/pi*180*3600-Ldu*3600-Lfen*60;Bdu=int(BB0/pi*180);Bfen=int((BB0/pi*180)*60-Bdu*60);Bmiao=BB0/pi*180*3600-Bdu*3600-Bfen*60;//将弧度转化为角度//printf("\n所选坐标系的转换结果：%d度%d分%lf秒 %d 度%d分%lf秒 \n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);fprintf(w,"%d°%d’%lf”%d°%d’%lf”\n",Bdu,Bfen,Bmiao,Ldu,Lfen,Lmiao);//将结果输出到文本文件//}fclose(r);fclose(w);system("pause");return 0;}double sort(double a,double b){double e,e1;e=sqrt(1-(b/a)*(b/a));e1=sqrt((a/b)*(a/b)-1);double m0,m2,m4,m6,m8;double a0,a2,a4,a6,a8;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=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7*m8/16;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;B1=g[i].x/a0;do{F=-a2*sin(2*B1)/2+a4*sin(4*B1)/4-a6*sin(6*B1)/6+a8*s in(8*B1)/8;B2=(g[i].x-F)/a0;B3=B1;B1=B2;} while(fabs(B3-B2)>10e-10);//利用迭代算法求解底点纬度//return B2;}。
vb课程课件测绘程序设计8(八)解析
直到
Bi1小 B于i 某一个指定数值，即可停止迭代。
F(B) a(1 e2 )[AarcB Bsin 2B Csin 4B Dsin 6B Esin 8B F sin10B Gsin12B]
F(B) a(1 e2 )[A 2Bcos2B 4Ccos4B 6Dcos6B 8Ecos8B 10Fcos10B 12Gcos12B]
2 投影变形的处理方法
➢通过改变 Hm从而选择合适的高程参考面，将抵偿分 带投影变形，这种方法通常称为抵偿投影面的高斯正 形投影；
➢通过改变ym，从而对中央子午线作适当移动，来抵
偿由高程面的边长归算到参考椭球面上的投影变形， 这就是通常所说的任意带高斯正形投影；
➢通过既改变 Hm（选择高程参考面），又改变 ym（移
测绘程序设计(八)
第八讲 高斯投影变换
主要内容 一、地图投影和正形投影 二、高斯投影与分带 三、高斯投影正算和反算 四、工程测量投影面与投影带选择 五、高斯投影换带计算
一、地图投影和正形投影
1、地图投影及其变形
几何投影——又叫透视投影，有中心投影、平行投影、垂直投影等。 数学投影——是数学的投影，建立椭球面大地坐标（B、L）与投影平 面上对应的坐标（x、y）之间的函数关系。
式中:
A 1 3 e2 45 e4 175 e6 11025 e8 43659 e10 693693 e12 4 64 256 16384 65536 1048576
B 3 e2 15 e4 525 e6 2205 e8 72765 e10 297297 e12 8 32 1024 4096 131072 524288
高斯投影除了在中央子午线上没有长度变形外，不在中央子午线
上的各点，其长度比都大于1，且离开中央子午线愈远，长度变形
高斯投影计算Vlisp程序开发
北京测绘 ＢｅｉｊｉｎｇＳｕｒｖｅｙｉｎｇａｎｄＭａｐｐｉｎｇ
Ｖｏｌ．３２ Ｎｏ．９ Ｓｅｐｔｅｍｂｅｒ２０１８
引文格式：岳江潮，欧阳永龙，靳宝．高斯投影计算 Ｖｌｉｓｐ程序开发［Ｊ］．北京测绘，２０１８，３２（９）：１０５９１０６３． 犇犗犐：１０．１９５８０／ｊ．ｃｎｋｉ．１００７３０００．２０１８．０９．０１３
四种坐标系椭球参数如表１所示。
［收稿日期］ ２０１８ ０５ １８ ［作者简介］ 岳江潮（１９７６－），男，陕西咸阳人，大学本科，高级工程师，从事工程测量工作。 犈犿犪犻犾：ｘｐ＿ｙｊｃ＠１６３．ｃｏｍ
１０６０
北京测绘
第３２卷 第９期
坐标系名称 １９５４北京坐标系 １９８０西安坐标系 ２０００国家大地坐标系 ＷＧＳ８４坐标系
［关键词］ Ｖｌｉｓｐ；高斯投影正反算；坐标转换 ［中图分类号］ Ｐ２０８ ［文献标识码］ Ａ ［文章编号］ １００７－３０００（２０１８）０９－１０５９－６
０ 引言
当前，测绘 工 作 常 用 的 坐 标 系 统 有 北 京 ５４、 西安８０、ＷＧＳ８４、国 家 ２０００ 坐 标 系，以 及 城 市、 矿区根据自身需要建立的地方独立坐标系、矿区 坐标系。随着测绘事业的迅速发展、计算机及信 息技术的 进 步，测 绘 ３Ｓ技 术、４Ｄ 产 品 的 应 用， 测绘 成 果 、地 理 信 息 数 据 的 重 要 性 也 日 益 凸 显 ， 在城市规 划、矿 业 权 核 查、城 镇 土 地 调 查、地 理 国情普查、压 覆 矿 产 报 告 计 算、矿 区 兼 并 重 组、 编制矿井 关 闭 报 告 的 工 作 中，经 常 需 要 进 行 坐 标转 换 ，形 成 统 一 的 地 理 信 息 数 据 ，满 足 不 同 任 务的需求。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
高斯投影坐标正反算V B程序Jenny was compiled in January 2021高斯投影坐标正反算学院：班级：学号：姓名：课程名称:指导老师：实验目的：1.了解高斯投影坐标正反算的基本思想；2.学会编写高斯正反算程序，加深了解。
实验原理：高斯投影正算公式中应满足的三个条件：1. 中央子午线投影后为直线；2. 中央子午线投影后长度不变；3. 投影具有正形性质，即正形投影条件。
高斯投影反算公式中应满足的三个条件：1. x坐标轴投影成中央子午线，是投影的对称轴；2. x轴上的长度投影保持不变；3. 正形投影条件，即高斯面上的角度投影到椭球面上后角度没有变形，仍然相等。
操作工具：计算机中的代码：Dim a As Double, b As Double, x As Double, y As Double, y_#Dim l_ As Double, b_ As Double, a0#, a2#, a4#, a6#, a8#, m2#, m4#, m6#, m8#, m0#, l0#, e#, e1#Dim deg1 As Double, min1 As Double, sec1 As Double, deg2 As Double, min2 As Double, sec2 As DoublePrivate Sub Command1_Click()Dim x_ As Double, t#, eta#, N#, W#, k1#, k2#, ik1%, ik2%, dh% deg1 = Valmin1 = Valsec1 = Valdeg2 = Valmin2 = Valsec2 = Vall_ = (deg1 * 3600 + min1 * 60 + sec1) / 206265b_ = (deg2 * 3600 + min2 * 60 + sec2) / 206265dh = Valk1 = ((l_ * 180 / + 3) / 6)k2 = (l_ * 180 / / 3)ik1 = Round(k1, 0)ik2 = Round(k2, 0)If dh = 6 Thenl0 = 6 * ik1 - 3ElseIf dh = 3 Thenl0 = 3 * ik2ElseMsgBox "error", 48, "error": Exit SubEnd IfEnd Ifl = l_ - l0 * / 180e = Sqr(a * a - b * b) / am0 = a * (1 - e * e)m2 = e * e * m0 * 3 / 2m4 = e * e * m2 * 5 / 4m6 = m4 * e * e * 7 / 6m8 = e * e * m6 * 9 / 8a0 = m0 + m2 / 2 + m4 * 3 / 8 + m6 * 5 / 16 + m8 * 35 / 128a2 = m2 / 2 + m4 / 2 + m6 * 15 / 32 + m8 * 7 / 16a4 = m4 / 8 + m6 * 3 / 16 + m8 * 7 / 32a6 = m6 / 32 + m8 / 16a8 = m8 / 128x_ = a0 * b_ - a2 * Sin(2 * b_) / 2 + a4 * Sin(4 * b_) / 4 - Sin(6 * b_) * a6 / 6 + Sin(8 * b_) * a8 / 8t = Tan(b_)e1 = Sqr((a * a - b * b) / (b * b))eta = Sqr(e1 * e1 * Cos(b) * Cos(b))W = Sqr(1 - e * e * Sin(b_) * Sin(b_))N = a / Wx = x_ + N * Sin(b_) * Cos(b_) * l * l / 2 + N * Sin(b_) *Cos(b_) ^ 3 * (5 - t * t + 9 * eta * eta + 4 * eta ^ 4) * l ^ 4 / 24 + N * Sin(b_) * Cos(b_) ^ 5 * (61 - 58 * t * t + t ^ 4) * l ^ 6 / 720y = N * Cos(b_) * l + N * Cos(b_) ^ 3 * (1 - t * t + eta * eta) * l * l * l / 6 + N * Cos(b_) ^ 5 * (5 - 18 * t * t + t ^ 4 + 14 * eta * eta - 58 * eta * eta * t * t) * l ^ 5 / 120Text7 = xIf dh = 6 Theny_ = y + 500000 + 1000000 * ik1ElseIf dh = 3 Theny_ = y + 500000 + 1000000 * ik2ElseMsgBox "error", 48, "error": Exit SubEnd IfEnd IfText8 = y_End SubPrivate Sub Command2_Click()Dim bf#, j%, Wf#, Vf#, Nf#, Mf#, c#, tf#, etaf#, dh%, ik%x = Valy_ = Vale = Sqr((a * a - b * b) / (a * a))m0 = a * (1 - e * e)m2 = e * e * m0 * 3 / 2m4 = e * e * m2 * 5 / 4m6 = m4 * e * e * 7 / 6m8 = e * e * m6 * 9 / 8a0 = m0 + m2 / 2 + m4 * 3 / 8 + m6 * 5 / 16 + m8 * 35 / 128 a2 = m2 / 2 + m4 / 2 + m6 * 15 / 32 + m8 * 7 / 16a4 = m4 / 8 + m6 * 3 / 16 + m8 * 7 / 32a6 = m6 / 32 + m8 / 16a8 = m8 / 128bf = x / a0For j = 1 To 10bf = (x + a2 / 2 * Sin(2 * bf) - a4 * Sin(4 * bf) / 4 + a6 * Sin(6 * bf) / 6 - a8 * Sin(8 * bf) / 8) / a0Next je1 = Sqr(a * a - b * b) / bVf = Sqr(1 + e1 * e1 * Cos(bf) * Cos(bf))Wf = Sqr(1 - e * e * Sin(bf) * Sin(bf))Nf = a / Wfc = a * a / bMf = c / Vf ^ 3tf = Tan(bf)e1 = Sqr((a * a - b * b) / (b * b))etaf = Sqr(e1 * e1 * Cos(bf) * Cos(bf))ik = Int(y_ / 1000000)y = y_ - ik * 1000000 - 500000b_ = bf - tf * y * y / (2 * Mf * Nf) + tf * (5 + 3 * tf *tf + etaf * etaf - 9 * etaf * etaf * tf * tf) * y * y * y * y / (24 * Mf * Nf * Nf * Nf) + tf * (61 + 90 * tf * tf + 45 * tf * tf * tf * tf) * y * y * y * y * y * y / (720 * Mf * Nf * Nf * Nf * Nf * Nf)l = y / (Nf * Cos(bf)) - (1 + 2 * tf * tf + etaf * etaf) * y * y * y / (6 * Nf * Nf * Nf * Cos(bf)) + (5 + 28 * tf * tf + 24 * tf * tf * tf * tf + 6 * etaf * etaf + 8 * etaf * etaf * tf * tf) * y * y * y * y * y / (120 * Nf * Nf * Nf * Nf * Nf * Cos(bf))dh = ValIf dh = 6 Thenl0 = 6 * ik - 3ElseIf dh = 3 Thenl0 = 3 * ikElseMsgBox "error", 48, "error": Exit Sub End IfEnd Ifl_ = l0 * / 180 + lsec1 = l_ * 206265deg1 = Int(sec1 / 3600)min1 = Int((sec1 - deg1 * 3600) / 60) sec1 = sec1 - deg1 * 3600 - min1 * 60 sec2 = b_ * 206265deg2 = Int(sec2 / 3600)min2 = Int((sec2 - deg2 * 3600) / 60) sec2 = sec2 - deg2 * 3600 - min2 * 60 Text1 = deg1Text2 = min1Text3 = Round(sec1, 5)Text4 = deg2Text5 = min2Text6 = Round(sec2, 5) End SubPrivate Sub Command3_Click() EndEnd SubPrivate Sub Command4_Click() = ""= ""= ""= ""= ""= ""= ""= ""= ""= ""End SubPrivate Sub Option1_Click() a = 6378245b =End SubPrivate Sub Option2_Click()a = 6378140b = 6356755.End SubPrivate Sub Option3_Click()a = 6378137b =End SubPrivate Sub Option4_Click()a = 6378137b =End SubPrivate Sub Option5_Click()a = Val(InputBox("a", "plsase input"))b = Val(InputBox("b", "please input")) End Sub界面截图如下：使用方法：1.如下所示，高斯投影坐标正算时，在经度纬度相应的文本框里输入经纬度，再在度号对应的文本框里输入是六度带还是三度带，最后按坐标正算按钮即可，答案会显示在X，Y相应的文本框里；2.如果要进行坐标反算，则输入X，Y与度号，最后按坐标反算按钮即可得到所需的大地经纬度；3.注意选择需要的椭球。