高斯投影正反算公式
更新时间:2023-12-04 16:09:01 阅读量: 教育文库 文档下载
高斯投影坐标正反算
一、基本思想:
高斯投影正算公式就是由大地坐标(L,B)求解高斯平面坐标(x,y),而高斯投影反算公式则是由高斯平面坐标(x,y)求解大地坐标(L,B)。
二、计算模型:
基本椭球参数: 椭球长半轴a 椭球扁率f
椭球短半轴:b?a(1?f)
a2?b2椭球第一偏心率 :e?
aa2?b2椭球第二偏心率 :e?? b高斯投影正算公式:此公式换算的精度为0.001m
x?X??NN2??sinBcosB?l?simBcos3B(5?t2?9?2?4?4)l??4242???24???
N5246??sinBcosB(61?58t?t)l720???6y?NN3223??cosB?l???cosB(1?t??)l???6???3N5242225???cosB(5?18t?t?14??58?t)l120???5
其中:角度都为弧度
B为点的纬度,l???L?L0,L为点的经度,L0为中央子午线经度;
N为子午圈曲率半径,N?a(1?esinB); t?tanB;
22?12?2?e?2cos2B
????180??3600
其中X为子午线弧长:
1616??X?a0B?sinBcosB?(a2?a4?a6)?(2a4?a6)sin2B?a6sin4B?
33??a0,a2,a4,a6,a8为基本常量,按如下公式计算: m23535?a?m??m?m?m8046?02816128??a?m2?m4?15m?7m68?2223216?m437? ?m6?m8?a4?81632?m6m8?a??632?16??a?m88?128?m0,m2,m4,m6,m8为基本常量,按如下公式计算:
m0?a(1?e2);m2?3279em0;m4?5e2m2;m6?e2m4;m8?e2m6; 268
高斯投影反算公式:此公式换算的精度为0.0001’’.
B?Bf??tftf2MfNfy?2tf24MfN3f?5?3t2f224??2f?9?ftf?y720MfN5f46y?61?90t2f?45tf?yyy322l??1?2t???ff?NfcosBf6N3cosBffy5422?5?28t2?2?f?24tf?6?f?8ftf?5120NfcosBfL?l?L0
其中:
L0为中央子午线经度。
Bf为底点纬度,也就是当x?X时的子午线弧长所对应的纬度。按照子午线弧长公式:
X?a0B?aaa2asin2B?4sin4B?6sin6B?8sin8B,迭代进行计算; 24681初始开始时设:Bf?Xa0 以后每次迭代按下式计算:
Bif?1?X(X?F(Bif))a0 a6a8a2a4iiiiF(B)??sin2Bf?sin4Bf?sin6Bf?sin8Bf2468if重复迭代至Bif?1?Bif??为止。
Nf?a(1?esinBf); Mf?a(1?e)(1?esinBf)
222?3222?12tf?tanBf;
22??2?ecosBf f
海福特椭球(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}
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*cby=(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*y=y+500000;
ow(sin(B),4)/3.0);
*cb*cb*cb*cb*(61-58*t*t+t*t*t*t)*l*l*l*l*l*l/720); ita-58*ita*ita*t*t)*l*l*l*l*l/120);
/**************高斯反算函数***************
输入 : 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 {
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;
double Bf1,Bfi0,Bfi1,FBfi;
} 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(Nfl=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
,3))-tf*(61+90*tf*tf+45*pow(tf,4))*pow(y,6)/(720*Mf*pow(Nf,5));
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
正在阅读:
高斯投影正反算公式12-04
陕西省殡葬服务收费管理办法(新修订版-陕价行发154号)11-20
循环流化床秸秆锅炉项目04-06
国家基金申请代码大全(2013版)04-12
固定资产管理制度(标准)09-09
最新教科版小学科学五年级下册《给冷水加热》公开课教学设计103-22
N6210型柴油机说明书11-05
三国杀武将台词09-20
电信光猫与TPLINK路由器如何连接05-16
- exercise2
- 铅锌矿详查地质设计 - 图文
- 厨余垃圾、餐厨垃圾堆肥系统设计方案
- 陈明珠开题报告
- 化工原理精选例题
- 政府形象宣传册营销案例
- 小学一至三年级语文阅读专项练习题
- 2014.民诉 期末考试 复习题
- 巅峰智业 - 做好顶层设计对建设城市的重要意义
- (三起)冀教版三年级英语上册Unit4 Lesson24练习题及答案
- 2017年实心轮胎现状及发展趋势分析(目录)
- 基于GIS的农用地定级技术研究定稿
- 2017-2022年中国医疗保健市场调查与市场前景预测报告(目录) - 图文
- 作业
- OFDM技术仿真(MATLAB代码) - 图文
- Android工程师笔试题及答案
- 生命密码联合密码
- 空间地上权若干法律问题探究
- 江苏学业水平测试《机械基础》模拟试题
- 选课走班实施方案
- 高斯
- 正反
- 公式
- 投影
- 新SQL-SERVER实验练习
- BRC巧克力及其制品HACCP计划(中英文) - 图文
- 桥梁上部构造T梁预制分项开工报告
- 帕尔马斯岛仲裁案
- 面试考务准备工作有关要求
- 奥鹏北语14秋《心理学》作业4满分答案
- 公安局110接处警系统 - 图文
- 15秋福师《教育学》在线作业
- 外国文学史名词解释
- 工作报告之结题报告教师评语
- Java常用类习题(附部分答案)
- 滨州市2014年3+4高职本科(试点)、
- 中国社会主义改革开放30年(92分)
- 2015年上《作物育种学各论》考试A答案
- 2017浙大《天然药物化学》在线作业
- 某校用水、用火、用电、用气等设施设备安全管理制度
- 护理考核标准
- 农药生产定点核准
- 小学六一儿童节节目串词与小学六年级六一儿童节主持开场白汇编
- 高中议论文满分作文(共8篇)