大地坐标及直角空间坐标转换计算公式Word文件下载.docx
《大地坐标及直角空间坐标转换计算公式Word文件下载.docx》由会员分享，可在线阅读，更多相关《大地坐标及直角空间坐标转换计算公式Word文件下载.docx（10页珍藏版）》请在冰豆网上搜索。
二高斯投影及高斯直角坐标系
1、高斯投影概述
高斯-克吕格投影的条件：
1.是正形投影；
2.中央子午线不变形
高斯投影的性质：
1.投影后角度不变；
2.长度比与点位有关.与方向无关；
3.离中央子午线越远变形越大
为控制投影后的长度变形.采用分带投影的方法。
常用3度带或6度带分带.城市或工程控制网坐标可采用不按3度带中央子午线的任意带。
2、高斯投影正算公式：
3、高斯投影反算公式：
1
坐标转换简介
坐标系统之间的坐标转换既包括不同的参心坐标之间的转换.或者不同的地心坐标系之间的转换.也包括参心坐标系与地心坐标系之间的转换以及相同坐标系的直角坐标与大地坐标之间的坐标转换.还有大地坐标与高斯平面坐标之间的转换。
在两个空间角直坐标系中.假设其分别为O--XYZ和O--XYZ.如果两个坐标系的原点相同.通过三次旋转.就可以使两个坐标系重合；
如果两个直角坐标系的原点不在同一个位置.通过坐标轴的平移和旋转可以取得一致；
如果两个坐标系的尺度也不尽一致.就需要再增加一个尺度变化参数；
而对于大地坐标和高斯投影平面坐标之间的转换.则需要通过高斯投影正算和高斯投影反算.通过使用中央子午线的经度和不同的参考椭球以及不同的投影面的选择来实现坐标的转换。
如何使用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——（X.Y.Z）84.空间大地坐标到空间直角坐标的转换。
2）（X.Y.Z）84——（X.Y.Z）54.坐标基准的转换.即Datum转换。
通常有三种转换方法：
七参数、简化三参数、Molodensky。
3）（X.Y.Z）54——（B.L）54.空间直角坐标到空间大地坐标的转换。
4）（B.L）54——（x.y）54.高斯投影正算。
从以上步骤不难看出.转换的关键是第二步.转换的参数。
鉴于我国曾使用不同的坐标基准（BJ54、State80、Correct54）.各地的重力值又有很大差异.所以很难确定一套适合全国且精度较好的转换参数。
在WGS-84坐标和54坐标之间是不存在一套转换参数可以全国通用的.在每个地方会不一样。
必须了解.在不同的椭球之间的转换是不严密的。
那么.两个椭球间的坐标转换应该是怎样的呢？
一般而言比较严密的是用七参数法.即3个平移因子（X平移.Y平移.Z平移）.3个旋转因子（X旋转.Y旋转.Z旋转）.一个比例因子（也叫尺度变化K）。
国内参数来源的途径不多.一般当地测绘部门会有。
通行的做法是：
在工作区内找三个以上的已知点.利用已知点的BJ54坐标和所测WGS84坐标.通过一定的数学模型.求解七参数。
若多选几个已知点.通过平差的方法可以获得较好的精度。
如果区域范围不大.最远点间的距离不大于30Km（经验值）.这可以用三参数.即只考虑3个平移因子（X平移.Y平移.Z平移）.而将旋转因子及比例因子（X旋转.Y旋转.Z旋转.尺度变化K）都视为0.所以三参数只是七参数的一种特例。
54和XX80也是两种不同的大地基准面.不同的参考椭球体.他们之间的转换也是同理。
在ArcGIS中提供了三参数、七参数转换法。
而在同一个椭球里的转换都是严密的.在同一个椭球的不同坐标系中转换需要用到四参数转换.举个例子.在XX既有54坐标又有XX坐标.在这两种坐标之间转换就用到四参数.计算四参数需要两个已知点
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]孔祥元、郭际明、X宗泉.《大地测量学基础》.XX大学.第一版.2001年9月
[2]李征航、黄劲松.《GPS测量与数据处理》.XX大学.第一版.2005年3月
[3]MAPGIS使用教程
大地坐标空间直角坐标转换
（2009-10-2221:
12:
41）
▼
标签：
杂谈
程序计算大地坐标与空间直角坐标转换
#include<
stdio.h>
math.h>
#definePI3.1415926535897932384626433832795
doublea,b,c,B,L,N,e,X,Y,Z,W,H;
intchoice;
doubleB1=0.0,B2=0.0;
doubledelta=0.0;
intmain（）
{
printf（"
pleaseinsertlongra\n"
）;
scanf（"
%f"
&
a）;
pleaseinsertshortrb\n"
b）;
e=sqrt（a*a-b*b）/a;
c=a*a/b;
1kongzhitodadi\n2daditokongzhi\0exit"
%d"
choice）;
while（choice!
=0）
{if（choice==2）
jingduL"
L）;
weiduB"
B）;
gaoduH"
H）;
W=sqrt（1-e*e*sin（B）*sin（B））;
N=a/W;
X=（N+H）*cos（B）*cos（L）;
Y=（N+H）*cos（B）*sin（L）;
Z=（N*（1-e*e）+H）*sin（B）;
X=%f,Y=%f,Z=%f"
X,Y,Z）;
}
if（choice==1）
zuobiao:
\nX="
X）;
\nY="
Y）;
\nZ="
Z）;
L=atan（Y/X）;
//
doublePP2=Z;
B1=atan（Z/sqrt（X*X+Y*Y））;
delta=PI/（180*60*60*1000.0）;
while（（B1-B2）>
=delta）
{
B2=B1;
B1=atan（（Z+N*e*e*sin（B1））/sqrt（X*X+Y*Y））;
B=B1;
H=Z/sin（B）-N*（1-e*e）;
L=%f,B=%f,H=%f"
elseprintf（"
entererror,pleaseenteragain"
return0;
/
高斯投影正、反算
//高斯投影正、反算
//////6度带宽
54年坐标系
//高斯投影由经纬度（Unit:
DD）反算大地坐标（含带号.Unit:
Metres）
voidGaussProjCal（doublelongitude,doublelatitude,double*X,double*Y）
intProjNo=0;
intZoneWide;
////带宽
double
longitude1,latitude1,longitude0,latitude0,X0,Y0,xval,yval;
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年XX坐标系参数
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））;
xval=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）;
X0=1000000L*（ProjNo+1）+500000L;
Y0=0;
xval=xval+X0;
yval=yval+Y0;
*X=xval;
*Y=yval;
//高斯投影由大地坐标（Unit:
Metres）反算经纬度（Unit:
DD）
voidGaussProjInvCal（doubleX,doubleY,double*longitude,double*latitude）
intProjNo;
doublelongitude1,latitude1,longitude0,latitude0,X0,Y0,xval,yval;
doublee1,e2,f,a,ee,NN,T,C,M,D,R,u,fai,iPI;
a=6378245.0;
f=1.0/298.3;
ProjNo=（int）（X/1000000L）;
//查找带号
longitude0=（ProjNo-1）*ZoneWide+ZoneWide/2;
//中央经线
X0=ProjNo*1000000L+500000L;
Y0=0;
xval=X-X0;
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=xval/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.0174532925199433//π
长半轴a=6378245.0;
扁率f=1.0/298.3;
长半轴a=6378140.0;
扁率f=1/298.257;
//80年XX坐标系参数
长半轴a=6378137m；
扁率f=1:
298.257223563。
//WGS-84坐标系
升级会员 