大地坐标及直角空间坐标转换计算公式.docx
《大地坐标及直角空间坐标转换计算公式.docx》由会员分享，可在线阅读，更多相关《大地坐标及直角空间坐标转换计算公式.docx（9页珍藏版）》请在冰豆网上搜索。
大地坐标及直角空间坐标转换计算公式
大地坐标与直角空间坐标转换计算公式
一、参心大地坐标与参心空间直角坐标转换
1名词解释：
A：
参心空间直角坐标系：
a）以参心0为坐标原点；
b）Z轴与参考椭球的短轴〔旋转轴〕相重合；
c）*轴与起始子午面和赤道的交线重合；
d）Y轴在赤道面上与*轴垂直.构成右手直角坐标系0-*YZ；
e）地面点P的点位用〔*.Y.Z〕表示；
B：
参心大地坐标系：
a）以参考椭球的中心为坐标原点.椭球的短轴与参考椭球旋转轴重合；
b）大地纬度B：
以过地面点的椭球法线与椭球赤道面的夹角为大地纬度B；
c）大地经度L：
以过地面点的椭球子午面与起始子午面之间的夹角为大地经度L；
d）大地高H：
地面点沿椭球法线至椭球面的距离为大地高H；
e）地面点的点位用〔B.L.H〕表示。
2参心大地坐标转换为参心空间直角坐标：
公式中.N为椭球面卯酉圈的曲率半径.e为椭球的第一偏心率.a、b椭球的长短半径.f椭球扁率.W为第一辅助系数
或
80椭球参数：
长半轴a=6378140±5〔m〕
短半轴b=6356755.2882m
扁率α=1/298.257
3参心空间直角坐标转换参心大地坐标
二高斯投影及高斯直角坐标系
1、高斯投影概述
高斯-克吕格投影的条件：
1.是正形投影；2.中央子午线不变形
高斯投影的性质：
1.投影后角度不变；2.长度比与点位有关.与方向无关；3.离中央子午线越远变形越大
为控制投影后的长度变形.采用分带投影的方法。
常用3度带或6度带分带.城市或工程控制网坐标可采用不按3度带中央子午线的任意带。
2、高斯投影正算公式：
3、高斯投影反算公式：
1坐标转换简介
坐标系统之间的坐标转换既包括不同的参心坐标之间的转换.或者不同的地心坐标系之间的转换.也包括参心坐标系与地心坐标系之间的转换以及一样坐标系的直角坐标与大地坐标之间的坐标转换.还有大地坐标与高斯平面坐标之间的转换。
在两个空间角直坐标系中.假设其分别为O--*YZ和O--*YZ.如果两个坐标系的原点一样.通过三次旋转.就可以使两个坐标系重合；如果两个直角坐标系的原点不在同一个位置.通过坐标轴的平移和旋转可以取得一致；如果两个坐标系的尺度也不尽一致.就需要再增加一个尺度变化参数；而对于大地坐标和高斯投影平面坐标之间的转换.则需要通过高斯投影正算和高斯投影反算.通过使用中央子午线的经度和不同的参考椭球以及不同的投影面的选择来实现坐标的转换。
如何使用ArcGIS实现WGS84经纬度坐标到BJ54高斯投影坐标的转换.这是很多从事GIS工作或者测绘工作者普遍遇到的问题。
本文目的在于帮助用户解决这个问题。
我们通常说的WGS-84坐标是指经纬度这种坐标表示方法.54坐标通常是指经过高斯投影的平面直角坐标这种坐标表示方法。
为什么要进展坐标转换.我们先来看两组参数.如表1所示：
表1BJ54与WGS84基准参数
参考椭球体
长半轴
短半轴
扁率
BJ54基准参数
Krasovsky_1940
6378245
6356863.0188
298.3
WGS84基准参数
WGS84
6378137
6356752.3142
298.257224
很显然.WGS84与BJ54是两种不同的大地基准面.不同的参考椭球体.因而两种地图下.同一个点的坐标是不同的.无论是三度带六度带坐标还是经纬度坐标都是不同的。
当要把GPS接收到的点〔WGS84坐标系统的〕叠加到BJ54坐标系统的底图上.那就会发现这些GPS点不能准确的在它该在的地方.即“与实际地点发生了偏移〞。
这就要求把这些GPS点从WGS84的坐标系统转换成BJ54的坐标系统了。
有关WGS84与BJ54的坐标转换问题.实质是WGS-84椭球体到BJ54椭球体的转换问题。
如果我们是需要把WGS84的经纬度坐标转换成BJ54的高斯投影坐标.那就还会涉及到投影变换问题。
因此.这个转换过程.一般的GPS数据处理软件都是采用下述步骤进展的：
1〕〔B.L〕84——〔*.Y.Z〕84.空间大地坐标到空间直角坐标的转换。
2〕〔*.Y.Z〕84——〔*.Y.Z〕54.坐标基准的转换.即Datum转换。
通常有三种转换方法：
七参数、简化三参数、Molodensky。
3〕〔*.Y.Z〕54——〔B.L〕54.空间直角坐标到空间大地坐标的转换。
4〕〔B.L〕54——〔*.y〕54.高斯投影正算。
从以上步骤不难看出.转换的关键是第二步.转换的参数。
鉴于我国曾使用不同的坐标基准〔BJ54、State80、Correct54〕.各地的重力值又有很大差异.所以很难确定一套适合全国且精度较好的转换参数。
在WGS-84坐标和54坐标之间是不存在一套转换参数可以全国通用的.在每个地方会不一样。
必须了解.在不同的椭球之间的转换是不严密的。
则.两个椭球间的坐标转换应该是怎样的呢.一般而言比拟严密的是用七参数法.即3个平移因子〔*平移.Y平移.Z平移〕.3个旋转因子〔*旋转.Y旋转.Z旋转〕.一个比例因子〔也叫尺度变化K〕。
国参数来源的途径不多.一般当地测绘部门会有。
通行的做法是：
在工作区找三个以上的点.利用点的BJ54坐标和所测WGS84坐标.通过一定的数学模型.求解七参数。
假设多项选择几个点.通过平差的方法可以获得较好的精度。
如果区域围不大.最远点间的距离不大于30Km〔经历值〕.这可以用三参数.即只考虑3个平移因子〔*平移.Y平移.Z平移〕.而将旋转因子及比例因子〔*旋转.Y旋转.Z旋转.尺度变化K〕都视为0.所以三参数只是七参数的一种特例。
54和80也是两种不同的大地基准面.不同的参考椭球体.他们之间的转换也是同理。
在ArcGIS中提供了三参数、七参数转换法。
而在同一个椭球里的转换都是严密的.在同一个椭球的不同坐标系中转换需要用到四参数转换.举个例子.在既有54坐标又有坐标.在这两种坐标之间转换就用到四参数.计算四参数需要两个点
2ArcGIS坐标转换例子
2.1应注意问题
使用ArcGIS如何实现WGS84经纬度坐标到BJ54高斯投影坐标的转换呢.在ArcGIS中.这个坐标转换步骤简化了.用户只需要两个步骤就能够直接从最初的WGS84经纬度坐标转换到BJ54高斯投影坐标。
这就是ArcGIS的强大之处。
接下来.我们做一个例子。
假设我们已经知道了7参数.应该如何操作呢.在具体的操作前.请大家一定注意以下三点：
WGS84的经纬度坐标值是用度来表示.而不能是度分秒表示
七参数的平移因子单位是米.旋转因子单位是秒.比例因子单位是百万。
在ArcGIS中.7参数法的名字是Coordinate_Frame方法。
有人在用ArcGIS进展不同椭球体间的坐标转换时.转换出来的结果不对.然后就写文章说变形如何如何.很可能是由于他们没有注意上面这三个关键的问题造成的。
2.2转换步骤
a、定义7参数的地理转换〔CreateCustomGeographicTransformation〕
在Arctool中翻开CreateCustomGeographicTransformation工具.如图1所示：
在弹出的窗口中.输入一个转换的名字.如wgs84ToBJ54。
在定义地理转换方法下面.在Method中选择适宜的转换方法如COORDINATE_FRAME.然后输入平移参数、旋转角度和比例因子.如图2所示：
b、投影变换
翻开工具箱下的ProjectionsandTransformations>Feature>Project.在弹出的窗口中输入要转换的数据以及OutputCoordinateSystem.然后输入第一步自定义的地理坐标系如wgs84ToBJ54.开场投影变换.如图3所示：
c:
\iknow\docshare\data\cur_work\photo.blog.sina..\showpic.html-blogid=a9091a330101aeal&url=s2.sinaimg./orignal/a9091a33gd8bf857c9e21
点击“确定〞.完成坐标转换。
3完毕语
我国现已启用新的坐标系统2000大地坐标系.2000大地坐标系与现行大地坐标系转换、衔接的过渡期仍需一段较长时期.在实际工作、工程中还遇到不同坐标系之间转换.本文针对在生产中从事测绘工作遇到的坐标转换问题提供解决方法和经历.希望对同行有所参考。
【参考文献】[1]孔祥元、郭际明、宗泉."大地测量学根底".大学.第一版.2001年9月
[2]征航、黄劲松."GPS测量与数据处理".大学.第一版.2005年3月
[3]MAPGIS使用教程
大地坐标空间直角坐标转换
（2009-10-2221:
12:
41）
▼
标签：
杂谈
程序计算大地坐标与空间直角坐标转换
*include
*include
*definePI3.795
doublea,b,c,B,L,N,e,*,Y,Z,W,H;
intchoice;
doubleB1=0.0,B2=0.0;
doubledelta=0.0;
intmain（）
{
printf（"pleaseinsertlongra\n"）;
scanf（"%f",&a）;
printf（"pleaseinsertshortrb\n"）;
scanf（"%f",&b）;
e=sqrt（a*a-b*b）/a;
c=a*a/b;
printf（"1kongzhitodadi\n2daditokongzhi\0e*it"）;
scanf（"%d",&choice）;
while（choice!
=0）
{if（choice==2）
{
printf（"jingduL"）;
scanf（"%f",&L）;
printf（"weiduB"）;
scanf（"%f",&B）;
printf（"gaoduH"）;
scanf（"%f",&H）;
W=sqrt（1-e*e*sin（B）*sin（B））;
N=a/W;
*=（N+H）*cos（B）*cos（L）;
Y=（N+H）*cos（B）*sin（L）;
Z=（N*（1-e*e）+H）*sin（B）;
printf（"*=%f,Y=%f,Z=%f",*,Y,Z）;
}
if（choice==1）
{
printf（"zuobiao:
\n*="）;
scanf（"%f",&*）;
printf（"zuobiao:
\nY="）;
scanf（"%f",&Y）;
printf（"zuobiao:
\nZ="）;
scanf（"%f",&Z）;
L=atan（Y/*）;
//doublePP2=Z;
B1=atan（Z/sqrt（***+Y*Y））;
delta=PI/（*60*60*1000.0）;
while（（B1-B2）>=delta）
{B2=B1;
B1=atan（（Z+N*e*e*sin（B1））/sqrt（***+Y*Y））;
}
B=B1;
H=Z/sin（B）-N*（1-e*e）;
printf（"L=%f,B=%f,H=%f"）;
}
elseprintf（"entererror,pleaseenteragain"）;
printf（"1kongzhitodadi\n2daditokongzhi\0e*it"）;
scanf（"%d",&choice）;
}
return0;
}
/
高斯投影正、反算
//高斯投影正、反算
//////6度带宽 54年坐标系
//高斯投影由经纬度（Unit:
DD）反算大地坐标（含带号.Unit:
Metres）
voidGaussProjCal（doublelongitude,doublelatitude,double**,double*Y）
{
intProjNo=0; intZoneWide; ////带宽
double longitude1,latitude1,longitude0,latitude0,*0,Y0,*val,yval;
double a,f,e2,ee,NN,T,C,A,M,iPI;
iPI=0.99433; ////3.98/.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））;
*val=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）;
*0=1000000L*（ProjNo+1）+500000L;
Y0=0;
*val=*val+*0; yval=yval+Y0;
**=*val;
*Y=yval;
}
//高斯投影由大地坐标（Unit:
Metres）反算经纬度（Unit:
DD）
voidGaussProjInvCal（double*,doubleY,double*longitude,double*latitude）
{
intProjNo; intZoneWide; ////带宽
doublelongitude1,latitude1,longitude0,latitude0,*0,Y0,*val,yval;
doublee1,e2,f,a,ee,NN,T,C,M,D,R,u,fai,iPI;
iPI=0.99433; ////3.98/.0;
a=6378245.0; f=1.0/298.3; //54年坐标系参数
////a=6378140.0; f=1/298.257; //80年坐标系参数
ZoneWide=6; ////6度带宽
ProjNo=（int）（*/1000000L）; //查找带号
longitude0=（ProjNo-1）*ZoneWide+ZoneWide/2;
longitude0=longitude0*iPI; //中央经线
*0=ProjNo*1000000L+500000L;
Y0=0;
*val=*-*0; 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=*val/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;
}
如果有需要程序的.可以直接跟我联系.呵呵
附：
高斯正反算参数
pi=0.0174532925※※0.99433//π
长半轴a=6378245.0; 扁率f=1.0/298.3; //54年坐标系参数
长半轴a=6378140.0; 扁率f=1/298.257; //80年坐标系参数
长半轴a=6378137m；扁率f=1:
298.257223563。
//WGS-84坐标系
升级会员 