(完整word版)高斯投影正反算公式_新
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
高斯投影坐标正反算
一、相关概念
大地坐标系由大地基准面和地图投影确定,由地图投影到特定椭圆柱面后在南北两极剪开展开而成,是对地球表面的逼近,各国或地区有各自的大地基准面,我国目前主要采用的基准面为:
1.WGS84基准面,为GPS基准面,17届国际大地测量协会上推荐,椭圆柱长半轴a=6378137m,短半轴
b=6356752.3142451m;
2.西安80坐标系,1975年国际大地测量协会上推荐,椭圆柱长半轴a=6378140m,短半轴b=6356755.2881575m;
3.北京54坐标系,参照前苏联克拉索夫斯基椭球体建立,椭圆柱长半轴a=6378245m, 短半轴b=6356863.018773m;
通常所说的高斯投影有三种,即投影后:
a)角度不变(正角投影),投影后经线和纬线仍然垂直;
b)长度不变;
c)面积不变;
大地坐标一般采用高斯正角投影,即在地球球心放一点光源,地图投影到过与中央经线相切的椭圆柱面上而成;可分带投影,按中央经线经度值分带,有每6度一带或每3度一带两种(起始带中央经线
经度为均为3度,即:6度带1带位置0-6度,3度带1带位置1.5-4.5 度),即所谓的高斯-克吕格投影。
图表11高斯投影和分带
地球某点经度(L)为过该点和地球自转轴的半圆与子午线所在半圆夹角,东半球为东经,西半球为西经;地球某点纬度(B)为所在水平面法线与赤道圆面的线面角。
正算是已知大地坐标(L,B),求解高斯平面坐标(X,Y),为确保Y值为正,Y增加500公里;反算则是由高斯平面坐标(X,Y)求解大地坐标(L,B)。
二、计算模型:
地球椭球面由椭圆绕地球自转轴旋转180度而成。
图表 1 椭圆
椭圆长半轴a,椭圆短半轴b, 椭圆方程:
(1)
图表2椭球面
椭球面方程:
/***************************************
与网上充斥的将函数关系先展开为泰勒级数,再依据投影规则确定各参数不同,本文直接依据空间立体三角函数关系得出结果。
*****/
(一)正算
由图表1,
由方程式(1),
令,可得
在图表2中,
,则
由椭圆方程,令
可知:
正算依据公式(4)、(5)、(6)、(7)得到结果,其中a:地球椭球长半轴;
b:地球椭球短半轴;
B:该点纬度;
L:该点经度减去中央经线L0后的值;
X:大地x坐标值;
Y: 大地y坐标值。
(6)式积分按积分原理由计算机求积分。
(二)反算
由式(4)可得,
三、程序代码函数:
/************高斯投影正算函数***************
输入: double a ,double b,(m_B,m_L)为大地坐标, L0为带号(6度带),(x,y)为高斯平面坐标,y加上了500000常量
返回:none
******************************************/
void gaosiforward(double a,double b,double m_B,double m_L,double L0,double &x,double &y)
{
double B=m_B*3.14159265/180;//换算成弧度
double L=(m_L-6.0*L0+3.0)*3.14159265/180; //换算成弧度
double xita=atan(b*b*tan(B)/a/a/cos(L));
double dxita=0.000001;
double xi=dxita;
x=0.0;
double c=a*a/b/b;
while(xi { x+=dxita/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi)); xi+=dxita; } x*=a; y=a*cos(xita)*tan(L)/sqrt(c*sin(xita)*sin(xita)+cos(xita)*c os(xita)); y+=500000.0; } /**************高斯反算函数*************** 输入: double a ,double b,(B,L)为大地坐标, L0为带号(6度带),(x,y)为高斯平面坐标,y加上了500000常量 *返回:none *****************************/ void gaosibackward(double a,double b,double x,double y,double L0,double &B,double &L) { double dxi=0.000001; double xi=dxi; double X=0.0; double c= a* a/ b/b; while(X< x/ a) { X+=dxi/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi)); xi+=dxi; } double r= a/sqrt(c*sin(xi)*sin(xi)+cos(xi)*cos(xi)); double Y=y-500000.0; L=atan(Y/r/cos(xi)); L=L*180/3.14159265-3.0+6.0*L0; B=atan(m_a*m_a*tan(xi)*cos(L)/m_b/m_b); B=B*180/3.14159265; }