高斯投影坐标正反算
一、基本思想：
高斯投影正算公式就是由大地坐标（L ，B ）求解高斯平面坐标（x ，y ），而高斯投影反算公式则是由高斯平面坐标（x ，y ）求解大地坐标（L ，B ）。
二、计算模型：
基本椭球参数：
椭球长半轴a
椭球扁率f
椭球短半轴：(1)b a f =-
椭球第一偏心率
：e a
= 椭球第二偏心率
：e b
'=高斯投影正算公式：此公式换算的精度为0.001m
64256
442234
22)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+''++-''+''?''+
=ρηηρρ 52224255
32233
)5814185(cos 120)1(cos 6cos l t t t B N l t B N l B N y ''-++-''+''+-''+''?''=ηηρηρρ
其中：角度都为弧度
B 为点的纬度，0l L L ''=-，L 为点的经度，0L 为中央子午线经度； N 为子午圈曲率半径，1222
(1sin )N a e B -=-；
tan t B =； 222cos e B η'=
1803600ρπ
''=*
其中X 为子午线弧长：
2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ??=--++-+????
02468,,,,a a a a a 为基本常量，按如下公式计算：
2004682426844686868
83535281612815722321637816323216128m a m m m m m m a m m m a m m m m a m a ?=++++???=+++???=++???=+???=??
02468,,,,m m m m m 为基本常量，按如下公式计算：
22222020426486379(1);;5;;268
m a e m e m m e m m e m m e m =-====；
高斯投影反算公式：此公式换算的精度为0.0001’’.
()()()()22222432465
3
2235
2422250
53922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f f
f f f f f f f
t t B B y t t y
M N M N t y t t y
M N y y l t N B N B y t t t N B L l L ηηηηη=-
+++--++=-+++++++=+
其中： 0L 为中央子午线经度。
f B 为底点纬度，也就是当x X =时的子午线弧长所对应的纬度。按照子午线弧长公式：68240sin 2sin 4sin 6sin82468
a a a a X a B B B B B =-+-+，迭代进行计算； 初始开始时设：10f B X a =
以后每次迭代按下式计算：
10
6824(())()sin 2sin 4sin 6sin82468i
i
f f i
i i i i f
f f f f B X X F B a a a a a F B B B B B +=-=-+-+
重复迭代至1i
i
f f B B ε+-<为止。
1222
(1sin )f f N a e B -=-；
32222(1)(1sin )f f M a e e B -=--
tan f f t B =；
222cos f f e B η'=
海福特椭球(1910) 我国52年以前基准椭球 a=6378388m b=6356911.9461279m α=0.33670033670
克拉索夫斯基椭球(1940 Krassovsky) 北京54坐标系基准椭球
a=6378245m b=6356863.018773m α=0.33523298692
1975年I.U.G.G 推荐椭球(国际大地测量协会1975) 西安80坐标系基准椭球
a=6378140m b=6356755.2881575m α=0.0033528131778
WGS-84椭球(GPS 全球定位系统椭球、17届国际大地测量协会) WGS-84 GPS 基准椭球
a=6378137m b=6356752.3142451m α=0.00335281006247
三、程序代码函数：
/************高斯投影正算函数***************
输入 : double a ,f 椭球参数，B,L 为大地坐标,L0为中央子午线的经度，单位为弧度，x,y 为高斯平面坐标,y 加上了500000常量
返回：none
******************************************/
void gaosiforward (double a ,double f ,double B ,double L ,double L0,double &x ,double &y ) {
double b , c ,e1, e2; //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率
double l , W ,N , M , daihao ;//W 为常用辅助函数，N 为子午圈曲率半径，M 为卯酉圈曲率半径
double X ;//子午线弧长，高斯投影的坐标
double ruo , ita , sb , cb ,t ;
double m [5],n [5];
//计算一些基本常量
{
b =a *(1-f );
e1=sqrt (a *a -b *b )/a ;
e2=sqrt (a *a -b *b )/b ;
c =a *a /b ;
m [0]=a *(1-e1*e1); m [1]=3*(e1*e1*m [0])/2.0;
m[2]=5*(e1*e1*m[1])/4.0;
m[3]=7*(e1*e1*m[2])/6.0;
m[4]=9*(e1*e1*m[3])/8.0;
n[0]=m[0]+m[1]/2+3*m[2]/8+5*m[3]/16+35*m[4]/128;
n[1]=m[1]/2+m[2]/2+15*m[3]/32+7*m[4]/16;
n[2]=m[2]/8+3*m[3]/16+7*m[4]/32;
n[3]=m[3]/32+m[4]/16;
n[4]=m[4]/128; /////by kjh 2014.5.22 把改成了
}
//由纬度计算子午线弧长
{
X=n[0]*B-sin(B)*cos(B)*((n[1]-n[2]+n[3])+(2*n[2]-(16*n[3]/3.0))*sin(B)*sin(B)+16*n[3]*p ow(sin(B),4)/3.0);
}
l=L-L0;//弧度
ita=e2*cos(B);
sb=sin(B);
cb=cos(B);
W=sqrt(1-e1*e1*sb*sb);
N=a/W;
t=tan(B);
ruo=(180/Pi)*3600;
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);
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);
y=y+500000;
}
/**************高斯反算函数***************
输入 : double a ,f 椭球参数， x,y为高斯平面坐标，L0为中央子午线的经度； B,L为大地坐标，单位为弧度
*返回：none
*****************************/
void gaosibackward(double a,double f,double x,double y,double L0,double &B,double &L)
{
double b, c,e1, e2; //短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率
double Bf,itaf,tf,Nf,Mf,Wf;
double l;
double m[5],n[5];
y=y-500000;
//计算一些基本常量
{
b=a*(1-f);
e1=sqrt(a*a-b*b)/a;
e2=sqrt(a*a-b*b)/b;
c=a*a/b;
m[0]=a*(1-e1*e1);
m[1]=3*(e1*e1*m[0])/2.0;
m[2]=5*(e1*e1*m[1])/4.0;
m[3]=7*(e1*e1*m[2])/6.0;
m[4]=9*(e1*e1*m[3])/8.0;
n[0]=m[0]+m[1]/2+3*m[2]/8+5*m[3]/16+35*m[4]/128;
n[1]=m[1]/2+m[2]/2+15*m[3]/32+7*m[4]/16;
n[2]=m[2]/8+3*m[3]/16+7*m[4]/32;
n[3]=m[3]/32+m[4]/16;
n[4]=m[4]/128;
}
//计算Bf
{
double Bf1,Bfi0,Bfi1,FBfi;
Bf1=x/n[0];
Bfi0=Bf1;
Bfi1=0;
FBfi=0;
int num=0;
do
{
num=0;
FBfi=0.0-n[1]*sin(2*Bfi0)/2.0+n[2]*sin(4*Bfi0)/4.0-n[3]*sin(6*Bfi0)/6.0;
Bfi1=(x-FBfi)/n[0];
if (fabs(Bfi1-Bfi0)>(Pi*pow(10.0,-8)/(36*18)))
{
num=1;
Bfi0=Bfi1;
}
} while (num==1);
Bf=Bfi1;
}
tf=tan(Bf);
Wf=sqrt(1-e1*e1*sin(Bf)*sin(Bf));
Nf=a/Wf;
Mf=a*(1-e1*e1)/(Wf*Wf*Wf);
itaf=e2*cos(Bf);
B=Bf-tf*y*y/(2*Mf*Nf)+tf*(5+3*tf*tf+itaf*itaf-9*itaf*itaf*tf*tf)*pow(y,4)/(24*Mf*pow(Nf ,3))-tf*(61+90*tf*tf+45*pow(tf,4))*pow(y,6)/(720*Mf*pow(Nf,5));
l=y/(Nf*cos(Bf))-(1+2*tf*tf+itaf*itaf)*pow(y,3)/(6*pow(Nf,3)*cos(Bf))+(5+28*tf*tf+24*po w(tf,4)+6*itaf*itaf+8*itaf*itaf*tf*tf)*pow(y,5)/(120*pow(Nf,5)*cos(Bf));
L=l+L0;
}
2014-5-22
' 输入: double a ,f 椭球参数，B,L为大地坐标,L0为中央子午线的经度，单位为弧度，x,y为高斯平面坐标,y加上了常量
Private Function gaosiforward(ByVal a As Double, ByVal f As Double, ByVal B As Double, ByVal L As Double, ByVal L0 As Double) As Double()
Dim x, y, xy(2) As Double
Dim bb, c, e1, e2 As Double'短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率
Dim ll, W, N, M, daihao As Double'W为常用辅助函数，N为子午圈曲率半径，M为卯酉圈曲率半径
Dim xx As Double'子午线弧长，高斯投影的坐标
Dim ruo, ita, sb, cb, t As Double
Dim mm(5), nn(5) As Double
bb = a * (1 - f)
e1 = Math.Sqrt(a * a - bb * bb) / a
e2 = Math.Sqrt(a * a - bb * bb) / bb
c = a * a / bb
mm(0) = a * (1 - e1 * e1)
mm(1) = 3 * (e1 * e1 * mm(0)) / 2.0
mm(2) = 5 * (e1 * e1 * mm(1)) / 4.0
mm(3) = 7 * (e1 * e1 * mm(2)) / 6.0
mm(4) = 9 * (e1 * e1 * mm(3)) / 8.0
nn(0) = mm(0) + mm(1) / 2 + 3 * mm(2) / 8 + 5 * mm(3) / 16 + 35 * mm(4) / 128
nn(1) = mm(1) / 2 + mm(2) / 2 + 15 * mm(3) / 32 + 7 * mm(4) / 16
nn(2) = mm(2) / 8 + 3 * mm(3) / 16 + 7 * mm(4) / 32
nn(3) = mm(3) / 32 + mm(4) / 16
nn(4) = mm(4) / 128
xx = nn(0) * B - Sin(B) * Cos(B) * ((nn(1) - nn(2) + nn(3)) + (2 * nn(2) - (16 * nn(3) / 3.0)) * Sin(B) * Sin(B) + 16 * nn(3) * Pow(Sin(B), 4) / 3.0)
ll = L - L0 '弧度
ita = e2 * Cos(B)
sb = Sin(B)
cb = Cos(B)
W = Sqrt(1 - e1 * e1 * sb * sb)
N = a / W
t = Tan(B)
ruo = (180 / PI) * 3600
x = (xx + N * sb * cb * ll * ll / 2 + N * sb * cb * cb * cb * (5 - t * t + 9 * ita * ita + 4 * ita * ita * ita * ita) * ll * ll * ll * ll / 24 + N * sb * cb * cb * cb * cb * cb * (61 - 58 * t * t + t * t * t * t) * ll * ll * ll * ll * ll * ll / 720)
y = (N * cb * ll + N * cb * cb * cb * (1 - t * t + ita * ita) * ll * ll * ll / 6 + N *
cb * cb * cb * cb * cb * (5 - 18 * t * t + t * t * t * t + 14 * ita * ita - 58 * ita * ita * t * t) * ll * ll * ll * ll * ll / 120)
y = y + 500000
xy(0) = x
xy(1) = y
Return xy
End Function
Private Function gaosibackward(ByVal a As Double, ByVal f As Double, ByVal x As Double, ByVal y As Double, ByVal L0 As Double) As Double()
Dim b, l, bl(2) As Double
Dim bb, c, e1, e2 As Double'短半轴，极点处的子午线曲率半径，第一偏心率，第二偏心率Dim Bf, itaf, tf, Nf, Mf, Wf As Double
Dim ll As Double
Dim m(5), n(5) As Double
y = y - 500000
bb = a * (1 - f)
e1 = Sqrt(a * a - bb * bb) / a
e2 = Sqrt(a * a - bb * bb) / bb
c = a * a / bb
m(0) = a * (1 - e1 * e1)
m(1) = 3 * (e1 * e1 * m(0)) / 2.0
m(2) = 5 * (e1 * e1 * m(1)) / 4.0
m(3) = 7 * (e1 * e1 * m(2)) / 6.0
m(4) = 9 * (e1 * e1 * m(3)) / 8.0
n(0) = m(0) + m(1) / 2 + 3 * m(2) / 8 + 5 * m(3) / 16 + 35 * m(4) / 128
n(1) = m(1) / 2 + m(2) / 2 + 15 * m(3) / 32 + 7 * m(4) / 16
n(2) = m(2) / 8 + 3 * m(3) / 16 + 7 * m(4) / 32
n(3) = m(3) / 32 + m(4) / 16
n(4) = m(4) / 128
'计算BF
Dim Bf1, Bfi0, Bfi1, FBfi As Double
Bf1 = x / n(0)
Bfi0 = Bf1
Bfi1 = 0
FBfi = 0
Dim num As Integer = 0
Do
num = 0
FBfi = 0.0 - n(1) * Sin(2 * Bfi0) / 2.0 + n(2) * Sin(4 * Bfi0) / 4.0 - n(3) * Sin(6 * Bfi0) / 6.0
Bfi1 = (x - FBfi) / n(0)
If (Abs(Bfi1 - Bfi0) > (PI * Pow(10.0, -8) / (36 * 18))) Then
num = 1
Bfi0 = Bfi1
End If
Loop While num = 1
Bf = Bfi1
tf = Tan(Bf)
Wf = Sqrt(1 - e1 * e1 * Sin(Bf) * Sin(Bf))
Nf = a / Wf
Mf = a * (1 - e1 * e1) / (Wf * Wf * Wf)
itaf = e2 * Cos(Bf)
b = Bf - tf * y * y / (2 * Mf * Nf) + tf * (5 + 3 * tf * tf + itaf * itaf - 9 * itaf * itaf * tf * tf) * Pow(y, 4) / (24 * Mf * Pow(Nf, 3)) - tf * (61 + 90 * tf * tf + 45 * Pow(tf, 4)) * Pow(y, 6) / (720 * Mf * Pow(Nf, 5))
ll = y / (Nf * Cos(Bf)) - (1 + 2 * tf * tf + itaf * itaf) * Pow(y, 3) / (6 * Pow(Nf, 3) * Cos(Bf)) + (5 + 28 * tf * tf + 24 * Pow(tf, 4) + 6 * itaf * itaf + 8 * itaf * itaf * tf * tf) * Pow(y, 5) / (120 * Pow(Nf, 5) * Cos(Bf))
l = ll + L0
bl(0) = hdtod(b)
bl(1) = hdtod(l)
Return bl
End Function
'度转换为弧度
Private Function dtohd(ByVal d As Double) As Double
Dim hd As Double
hd = d * PI / 180
Return hd
End Function
'弧度转换为度
Private Function hdtod(ByVal hd As Double) As Double
Dim d As Double
d = hd * 180 / PI
Return d
End Function
Private Function dfmtod(ByVal dfm As Double) As Double
Dim dfm2, dfm3(), du, fen, miao, miao1, miao2 As String
Dim duf As Double
dfm2 = dfm.ToString
dfm3 = dfm2.Split(".")
du = dfm3(0)
fen = dfm3(1).Substring(0, 2)
miao1 = dfm3(1).Substring(2, 2)
miao2 = dfm3(1).Substring(4)
miao = miao1 & "." & miao2
duf = Convert.ToInt32(du) + Convert.ToInt32(fen) / 60.0 + Convert.ToDouble(miao) / 3600 Return duf
End Function
'国家坐标系椭球参数2000,84,54,80 a2 = 6378137
f2 = 1 / 298.257222101 a84 = a2
f84 = 1 / 298.257223563 a54 = 6378245
f54 = 0.33523298692
a80 = 6378140
f80 = 0.0033528131778
高斯投影坐标正反算 V B程序 文件编码（GHTU-UITID-GGBKT-POIU-WUUI-8968）
高斯投影坐标正反算 学院： 班级： 学号： 姓名： 课程名称: 指导老师：
实验目的： 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 Double Private Sub Command1_Click() Dim x_ As Double, t#, eta#, N#, W#, k1#, k2#, ik1%, ik2%, dh% deg1 = Val min1 = Val sec1 = Val deg2 = Val min2 = Val sec2 = Val l_ = (deg1 * 3600 + min1 * 60 + sec1) / 206265 b_ = (deg2 * 3600 + min2 * 60 + sec2) / 206265 dh = Val k1 = ((l_ * 180 / + 3) / 6) k2 = (l_ * 180 / / 3) ik1 = Round(k1, 0) ik2 = Round(k2, 0) If dh = 6 Then l0 = 6 * ik1 - 3 Else
//84的椭球 final double a = 6378137; final double Alfa = 1.0 / 298.257223563; double centreL, x, y, b, e1, ee; double a0, a2, a4, a6, a8, Bf0; double[] Coeficient_a0 = new double[5]; double sinBf, cosBf; double FBf, Bf1, dB, bf; double c, v, Nf, Mf, tf; double itaf, dietaB, dietaL; double B, L; double dmsB, dmsL, dmsCentreL1; double radlat, radlon, radl0, l; double sb, cb, t, ita, X, N; public String MakeProject(double L, double B, double CentreLon) //高斯正算{ /*输入已知数据:经度\纬度\ 中央子午线*/ dmsB = B; dmsL = L; dmsCentreL1 = CentreLon; radlat = DMSTORAD(dmsB); radlon = DMSTORAD(dmsL); radl0 = DMSTORAD(dmsCentreL1); 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; a0a2a4a6a8(a, ee, Coeficient_a0); a0 = Coeficient_a0[0]; a2 = Coeficient_a0[1];
高斯投影正反算 班级：测绘九班C 姓名：塔娜 学号：2009301610323 指导老师：苏新洲 2011-11-02
高斯投影正反算编程 一、题目： 已知部分数据，根据高斯投影正反算思想进行编程，并采用克氏椭球，按3°或6°带投影。 正算：已知大地坐标B 、L， 二、已知数据： 正算： B=51.38439023 L=111.02131360 反算： x=5724004.723 y=19502559.920 三、计算结果： 正算结果： x=5724004.723 y=19502559.920 反算结果： B=51.38439023 L=111.02131360
（不予画出） 五、源代码： #include"gaosi.h" #include"math.h" #include"stdio.h" #include"tchar.h" #include"stdlib.h" #define pi 3.141592653589793 #define rho 206265 void Calculateellipse2plane(double B,double L); void Calculateplane2ellipse(double x,double y); double Dms2Rad(double Dms); double D2Dms(double D); double Dms2D(double Dms); int main() { double B=0,L=0; double x=0,y=0; int i=0; printf("如使用高斯投影坐标正算，请输入1；反算，请输入2\n"); scanf_s("%d",&i); if(i==1)
高斯投影坐标正反算编程报告 1. 编程思想 进行高斯投影坐标正反算的编程需要牵涉到大量的公式,为了使程序条理更清楚,各块的数据复用性更强,这里采取了结构化的编程思想。 程序由四大块组成。 GeodesyHomework 、cpp 文件用于存放main()函数,就是整个程序的入口。通过结构化的编程尽力使main()函数变得简单。 MyFunction 、h 与MyFunction 、cpp 用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。 Zhengsuan 、h 与Zhengsuan 、cpp 用于存放Zhengsuan 类,在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量,在类的构造函数中进行成员变量的初始化及正算计算。通过get 函数获得相应的正算结果。 Fansuan 、h 与Fansuan 、cpp 用于存放Fansuan 类,类似于Zhengsuan 类,Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量,在类的构造函数中进行成员变量的初始化及反算计算。通过get 函数获得相应的反算结果。 2. 计算模型 高斯投影正算公式 6 4256 4 42234 22)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+ ''++-''+''?''+=ρηηρρ 5 2224255 3 2233 )5814185(cos 120)1(cos 6cos l t t t B N l t B N l B N y ''-++-''+ ''+-''+''?''=ηηρηρρ 高斯投影反算公式 () () ()( ) 2 22425 52 23 36 4254 222232 8624285cos 12021cos 6cos 459061720935242f f f f f f f f f f f f f f f f f f f f f f f f f f f f f t t t B N y t B N y B N y l y t t y N M t y t t N M t y N M t B B ηηηηη+++++++-=++- -+++ -= 3. 程序框图
#include "stdafx.h" #include "iostream.h" #include "math.h" #include "stdio.h" #define P 206264.806247096355 #define PI 3.141592653589793 void GaosZ_fun() { printf("高斯投影的正算\n"); double l,L,B,n2,x,y,N,t,V,c,e2; double i,j,k,n,h,a0,a4,a6,a3,a5,cB2; int m; e2=0.006738525414683; c=6399698.901782711; B=17.33557339*3600/P; L=119.15521159*3600/P; l=L-111*3600/P // l=((m%6)*3600+n*60+h)/P; t=tan(B); n2=e2*cos(B)*cos(B); V=sqrt(1+n2); cB2=pow(cos(B),2); N=6399698.902-(21562.267-(108.973-0.612*cB2)*cB2)*cB2; // N=c/V; a0=32140.404-(135.3302-(0.7092-0.004*cB2)*cB2)*cB2; a4=(0.25+0.00252*cB2)*cB2-0.04166; a6=(0.166*cB2-0.084)*cB2; a3=(0.3333333+0.001123*cB2)*cB2-0.1666667; a5=0.0083-(0.1667-(0.1968+0.0040*cB2)*cB2)*cB2; // x=X+N*sin(B)*cos(B)*l*l/2+N*sin(B)*pow(cos(B),3)*(5-t*t+9*n2+4*n2*n2)*pow(l, 4)/24+N*sin(B)*pow(cos(B),5)*(61-58*t*t+pow(t,4))*pow(l,6)/720; // y=N*cos(B)*l+N*pow(cos(B),3)*(1-t*t+n2)*pow(l,5)/6+N*pow(cos(B),5)*(5-18*t*t +pow(t,4)+14*n2-58*n2*t*t)*pow(l,5)/120; x=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("x=%f\ny=%f\n",x,y); } void GaosF_fun() { printf("高斯投影的反算\n"); double B,Bf,Nf,b,b2,b3,b4,b5,Z,x,y,L0,l;
class Gauss { #region 高斯投影正反算 ///
高斯投影正反算编程一．高斯投影正反算基本公式（1）高斯正算基本公式 （2）高斯反算基本公式
以上主要通过大地测量学基础课程得到，这不进行详细的推导，只是列出基本公式指导编程的进行。 二．编程的基本方法和流程图 （1）编程的基本方法 高斯投影正反算基本上运用了所有的编程基本语句，本文中是利用C++语言进行基本的设计。高斯正算中对椭球参数和带宽的选择主要运用了选择语句。而高斯反算中除了选择语句的应用，在利用迭代算法求底点纬度还应用了循环语句。编程中还应特别注意相关的度分秒和弧度之间的相互转换，这是极其重要的。 （2）相关流程图 1）正算
选择带宽 3/6 度带 计算带号 输入大地坐标 B ，L 和经差 L0 6 度带 3 度带 选择椭球参数 计算带号 计算弧长 计算平面坐标 x,y 打印 x,y 开始 计算平面坐标 x,y 计算弧长 打印 x,y
开始 输入自然值坐标x,y 和经差L0 选择椭球参数 利用迭代算法 求解底点纬度 利用公式计算B 和L 打印B 和L 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;
高斯投影坐标正反算 一、相关概念 大地坐标系由大地基准面和地图投影确定，由地图投影到特定椭圆柱面后在南北两极剪开展开而成，是对地球表面的逼近，各国或地区有各自的大地基准面，我国目前主要采用的基准面为：基准面，为GPS基准面，17届国际大地测量协会上推荐，椭圆柱长半轴a=6378137m，短半轴b=； 2.西安80坐标系，1975年国际大地测量协会上推荐，椭圆柱长半轴a=6378140m，短半轴b=; 3.北京54坐标系，参照前苏联克拉索夫斯基椭球体建立,椭圆柱长半轴a=6378245m, 短半轴b=; 通常所说的高斯投影有三种，即投影后： a)角度不变（正角投影），投影后经线和纬线仍然垂直； b)长度不变； c)面积不变； 大地坐标一般采用高斯正角投影，即在地球球心放一点光源，地图投影到过与中央经线相切的椭圆柱面上而成；可分带投影，按中央经线经度值分带，有每6度一带或每3度一带两种(起始带中央经线经度为均为3度，即：6度带1带位置0-6度，3度带1带位置度)，即所谓的高斯-克吕格投影。
图表11高斯投影和分带 地球某点经度（L）为过该点和地球自转轴的半圆与子午线所在半圆夹角，东半球为东经，西半球为西经；地球某点纬度（B）为所在水平面法线与赤道圆面的线面角。 正算是已知大地坐标（L，B），求解高斯平面坐标（X，Y）,为确保Y值为正，Y增加500公里；反算则是由高斯平面坐标（X，Y）求解大地坐标（L，B）。 二、计算模型： 地球椭球面由椭圆绕地球自转轴旋转180度而成。 图表 1 椭圆 椭圆长半轴a，椭圆短半轴b, 椭圆方程：
(1) 图表2椭球面 椭球面方程： y2 a2+ x2 b2 + z2 a2 =1 /*************************************** 与网上充斥的将函数关系先展开为泰勒级数，再依据投影规则确定各参数不同，本文直接依据空间立体三角函数关系得出结果。 *****/ （一）正算 由图表1，
高斯投影坐标正反算编程报告 10021班 张鑫 学号：2010302590040 1. 编程思想 进行高斯投影坐标正反算的编程需要牵涉到大量的公式，为了使程序条理更清楚，各块的数据复用性更强，这里采取了结构化的编程思想。 程序由四大块组成。 GeodesyHomework.cpp 文件用于存放main()函数，是整个程序的入口。通过结构化的编程尽力使main()函数变得简单。 MyFunction.h 和MyFunction.cpp 用于存放计算过程中进行角度弧度换算时所要用到的一些自定的转换函数。 Zhengsuan.h 和Zhengsuan.cpp 用于存放Zhengsuan 类，在Zhengsuan 类中声明了高斯投影坐标正算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及正算计算。通过get 函数获得相应的正算结果。 Fansuan.h 和Fansuan.cpp 用于存放Fansuan 类，类似于Zhengsuan 类，Fansuan 类中声明了高斯投影坐标反算所要用到的所有变量，在类的构造函数中进行成员变量的初始化及反算计算。通过get 函数获得相应的反算结果。 2. 计算模型 高斯投影正算公式 6 4256 4 42234 22)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+ ''++-''+''?''+=ρηηρρ 5 2224255 3 2233)5814185(cos 120)1(cos 6cos l t t t B N l t B N l B N y ''-++-''+''+-''+''?''=ηηρηρρ 高斯投影反算公式 () () ()( ) 222425 52 2336 4254 2222 32 8624285cos 12021cos 6cos 459061720935242f f f f f f f f f f f f f f f f f f f f f f f f f f f f f t t t B N y t B N y B N y l y t t y N M t y t t N M t y N M t B B ηηηηη+++++++-=++- -+++ -=
高斯投影正反算 学院：资源与环境工程工程学院 专业：测绘工程 学号：X51414012 ：超 一、高斯投影概述 想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差围的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。 二、高斯投影坐标正算公式 1.高斯投影必须满足以下三个条件 1）中央子午线投影后为直线 2）中央子午线投影后长度不变 3）投影具有正形性质，即正形投影条件 2.高斯正算公式推导 1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。 2）由于高斯投影是换带投影，在每带经差ｌ是不大的，ｌ ρ是一个微小量，所以可 以将 X=X （ｌ，q ），Y=Y （l ，q ） 展开为经差为l 的幂级数，它可写成如下的形式 X=m 0+m 2l 2+m 4l 4 +…
Y=m 1l+m 3l 2+m 5l 5 +… 式中m 0,m1,m2,…是待定系数，他们都是纬度B 的函数。 3）由第三个条件：?y ?l =?x ?q 和?x ?l =-?y ?q ，将上式分别对l 和q 求偏导 234012342 3 4 01234........... x m m l m l m l m l y n n l n l n l n l =+++++=+++++ 可得到下式 03121234031212 34111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ?====?? ? ?=-=-=-=-? ? 经过计算可以得出 2322445246 32235242225 sin cos sin cos (594)224 sin cos (6158)720 cos cos (1) 6 cos (5181458)120N N x X B B l B B t l N B B t t l N y N B l B t l N B t t t l ηηηηη=+ ?+-+++-+=?+-++-++-三、高斯投影坐标反算公式推导 1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。 2.投影公式在底点处展开 12(,) (,) q f x y l f x y '='= 展开为
高斯投影坐标正反算 学院： 班级： 学号： 姓名： 课程名称: 指导老师：
实验目的： 1.了解高斯投影坐标正反算的基本思想； 2.学会编写高斯正反算程序，加深了解。 实验原理： 高斯投影正算公式中应满足的三个条件： 1. 中央子午线投影后为直线； 2. 中央子午线投影后长度不变； 3. 投影具有正形性质，即正形投影条件。 高斯投影反算公式中应满足的三个条件： 1. x坐标轴投影成中央子午线，是投影的对称轴； 2. x轴上的长度投影保持不变； 3. 正形投影条件，即高斯面上的角度投影到椭球面上后角度没 有变形，仍然相等。 操作工具： 计算机中的VB6.0 代码： 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 Double Private Sub Command1_Click() Dim x_ As Double, t#, eta#, N#, W#, k1#, k2#, ik1%, ik2%, dh% deg1 = Val(Text1.Text) min1 = Val(Text2.Text) sec1 = Val(Text3.Text) deg2 = Val(Text4.Text) min2 = Val(Text5.Text) sec2 = Val(Text6.Text) l_ = (deg1 * 3600 + min1 * 60 + sec1) / 206265 b_ = (deg2 * 3600 + min2 * 60 + sec2) / 206265 dh = Val(Text9.Text) k1 = ((l_ * 180 / 3.14159 + 3) / 6) k2 = (l_ * 180 / 3.14159 / 3) ik1 = Round(k1, 0) ik2 = Round(k2, 0) If dh = 6 Then l0 = 6 * ik1 - 3 Else If dh = 3 Then
高斯投影坐标正算公式 高斯投影坐标正反算公式 2.2.2. 1高斯投影坐标正算公式： B, x,y 高斯投影必须满足以下三个条件： ⑴中央子午线投影后为直线；⑵中央子午线投影后长度不变；⑶投影具有正形性质，即正形投影条件。 由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即 式中，x为的偶函数，y为的奇函数；，即， 如展开为的级数，收敛。 （2-10） 式中是待定系数，它们都是纬度B的函数。 由第三个条件知： 分别对和q求偏导数并代入上式 (2-11) 上两式两边相等，其必要充分条件是同次幂前的系数应相等，即
(2-12) (2-12)是一种递推公式，只要确定了就可依次确定其余各系数。 由第二条件知:位于中央子午线上的点，投影后的纵坐标x应等于投影前从赤道量至该点的子午线弧长X，即(2-10)式第一式中，当时有： (2-13) 顾及(对于中央子午线) 得： (2-14,15) (2-16) 依次求得并代入(2-10)式，得到高斯投影正算公式
(2-17) 2.2.2. 2高斯投影坐标反算公式 x,y B, 投影方程： (2-18) 高斯投影坐标反算公式推导要复杂些。 ⑴由x求底点纬度(垂足纬度),对应的有底点处的等量纬度，求x,y与 的关系式，仿照式有， 由于y和椭球半径相比较小(1/16.37)，可将展开为y的幂级数；又由于是对称投影，q必是y的偶函数，必是y的奇函数。 (2-19) 是待定系数，它们都是x的函数. 由第三条件知： ，
， (2-20) (2-19)式分别对x和y求偏导数并代入上式 上式相等必要充分条件，是同次幂y前的系数相等， 第二条件，当y=0时，点在中央子午线上，即x=X，对应的点称为底点，其纬度为底点纬度，也就是x=X时的子午线弧长所对应的纬度，设所对应的等量纬度为。也就是在底点展开为y的幂级数。 由(2-19)1式 依次求得其它各系数 (2-21) (2-21)1 ………… 将代入(2-19)1式得
//servey.h // #ifndef SERVEY_H // #define SERVEY_H #include
高斯投影正反算程序设计 一．程序设计流程 本程序的设计思路如下： （1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm 窗体程序形式。 （2），程序主要的算法来自于教材。但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。 （3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。 二．代码 using System; using Syst using https://www.doczj.com/doc/c6964139.html,ponentModel; using System.Data; using System.Drawing; using System.Text; namespace Gauss { public partial class Form1 : Form { //大地坐标 //Geodetic Coordinate public struct CRDGEODETIC { public double dLongitude; public double dLatitude; public double dHeight; } //笛卡尔坐标
//Cartesian Coordinate public struct CRDCARTESIAN { public double x; public double y; public double z; } public Form1() { InitializeComponent(); } private void button1_Click(object sender, EventArgs e) { double ee = 0; double a = 0; string tt; try { } catch { MessageBox.Show("Gauss Inverse: Choose datum error!"); return; } if (https://www.doczj.com/doc/c6964139.html,pareTo("克氏椭球")==0) { a = 6378245.00; } if (https://www.doczj.com/doc/c6964139.html,pareTo("WGS-84") == 0) { a = 6378137.00; } if (https://www.doczj.com/doc/c6964139.html,pareTo("1975国际椭球") == 0)
高斯投影正反算c代码 Coca-cola standardization office【ZZ5AB-ZZSYT-ZZ2C-ZZ682T-ZZT18】
高斯投影正反算程序设计 一．程序设计流程 本程序的设计思路如下： （1），程序采用VS08版本作为开发平台，并采用C#语言作为开发语言，设计为WindowsForm窗体程序形式。 （2），程序主要的算法来自于教材。但是本程序为了更加实用，添加了更多的解算基准，包括：WGS-84，国际椭球1975，克氏椭球，和2000国家大地坐标系。 （3），程序为了更方便的读取数据和输出数据，故需要自己定义了固定的数据输入格式和数据输出格式或形式，请老师注意查看。 二．代码 using System; using ; using ; using ; using ; using Gauss { public partial class Form1 : Form { double b = (a * a * (1 - ee * ee)); double c = a * a / b; double epp = ((a * a - b * b) / b / b); CRDGEODETIC pcrdGeo; CRDCARTESIAN pcrdCar; double midlong = 0;
//求X，Y和带号 = ; ytext = ; string temp = (0, 2); num = (temp); ytext = (0, 2); = (ytext) - 500000; try { tt = } catch { ("Choose 3/6 error!"); return; } if ("3度带") == 0) { midlong = num * 3 * pai / 180; } if ("6度带") == 0) { midlong = (6 * num - 3) * pai / 180; } b = (a * a * (1 - ee * ee)); c = a * a / b; epp = (a * a - b * b) / b; double m0, m2, m4, m6, m8; double a0, a2, a4, a6, a8; m0 = a * (1 - ee * ee);
高斯投影正反算编程 一、题目： 已知部分数据，根据高斯投影正反算思想进行编程，并采用克氏椭球，按3°或6°带投影。 正算：已知大地坐标B 、L， 二、已知数据： 正算： B=51.38439023 L=111.02131360 反算： x=5724004.723 y=19502559.920 三、计算结果： 正算结果： x=5724004.723 y=19502559.920 反算结果： B=51.38439023 L=111.02131360
五、源代码： #include"gaosi.h" #include"math.h" #include"stdio.h" #include"tchar.h" #include"stdlib.h" #define pi 3.141592653589793 #define rho 206265 void Calculateellipse2plane(double B,double L); void Calculateplane2ellipse(double x,double y); double Dms2Rad(double Dms); double D2Dms(double D); double Dms2D(double Dms);
int main() { double B=0,L=0; double x=0,y=0; int i=0; printf("如使用高斯投影坐标正算，请输入1；反算，请输入2\n"); scanf_s("%d",&i); if(i==1) { printf("请输入大地坐标纬度B(度分秒):\n"); scanf_s("%lf",&B); printf("请输入大地坐标经度L(度分秒):\n"); scanf_s("%lf",&L); Calculateellipse2plane(B,L); } if(i==2) { printf("请输入国家统一坐标x(m):\n"); scanf_s("%lf",&x); printf("请输入国家统一坐标y(m):\n"); scanf_s("%lf",&y); Calculateplane2ellipse(x,y); } return 0; } void Calculateellipse2plane(double B,double L) //高斯投影正算主体{ double l=0,Lo=0,a0=0,a3=0,a4=0,a5=0,a6=0,n=0,c=0; double x=0,y=0; double m=0,p=0,q=0; int N=0,i=0; //带号 printf("如使用6°带请输入1，使用3°带请输入2\n"); scanf_s("%d",&i); if(i==1) //已知a点在6°带的带号和中央子午线经度 { N=int(Dms2D(L)/6); Lo=6*N-3; } if(i==2) //已知a点在3°带的带号和中央子午线经度 { N=int((Dms2D(L)+1.5)/3);
高斯投影正、反算代码 //高斯投影正、反算 //////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)
§8.3高斯投影坐标正反算公式 任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。 8.3.1高斯投影坐标正算公式： B,l ? x,y 高斯投影必须满足以下三个条件： ①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。 由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为l 的偶函数，y 为l 的奇函数；0330'≤l ，即20/1/≈''''ρl ，如展开为l 的级数，收敛。 +++=++++=553316644220l m l m l m y l m l m l m m x （8-33） 式中 ,,10m m 是待定系数，它们都是纬度B 的函数。 由第三个条件知： q y l x l y q x ??-=????=??, (8-33)式分别对l 和q 求偏导数并代入上式 ----=++++++=+++553315 63424 42204 52 3164253l dq dm l dq dm l dq dm l m l m l m l dq dm l dq dm dq dm l m l m m (8-34) 上两式两边相等，其必要充分条件是同次幂l 前的系数应相等，即
dq dm m dq dm m dq dm m 231 20 13121? =? -== (8-35) (8-35)是一种递推公式，只要确定了 0m 就可依次确定其余各系数。 由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当0=l 时有： 0m X x == (8-36) 顾及(对于中央子午线) B V M r M B N dq dB M dB dX cos cos 2 ==== 得： B V c B N r dq dB dB dX dq dX dq dm m cos cos 01===?===(8-37,38) B B N dq dB dB dm dq dm m cos sin 2 2121112=?-=?-= (8-39) 依次求得6543,,,m m m m 并代入(8-33)式，得到高斯投影正算公式
高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系)，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。 //6度带宽54北京坐标系 //高斯投影由大地坐标(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; //带内大地坐标