大地坐标直角空间坐标转换计算公式.docx

1、大地坐标直角空间坐标转换计算公式. . .坐标与直角空间坐标转换计算公式一、参心坐标与参心空间直角坐标转换1名词解释：A：参心空间直角坐标系：a)以参心 0 为坐标原点；b)Z 轴与参考椭球的短轴（旋转轴）相重合；c)X 轴与起始子午面和赤道的交线重合；d) Y 轴在赤道面上与 X 轴垂直，构成右手直角坐标系 0-XYZ ；e)地面点 P 的点位用（ X ，Y ，Z ）表示；B：参心坐标系：a)以参考椭球的中心为坐标原点，椭球的短轴与参考椭球旋转轴重合；b) 纬度 B ：以过地面点的椭球法线与椭球赤道面的夹角为纬度 B；c) 经度 L ：以过地面点的椭球子午面与起始子午面之间的夹角为经度 L

2、；d) 高 H ：地面点沿椭球法线至椭球面的距离为高 H ；e)地面点的点位用（ B， L ， H ）表示。2参心坐标转换为参心空间直角坐标：X ( N H ) * cos B * cos LY (N H ) * cos B * sin LZ N * (1 e2 ) H * sin B公式中， N 为椭球面卯酉圈的曲率半径， e 为椭球的第一偏心率， a、 b 椭球的长短半径， f 椭球扁率， W 为第一辅助系数a2b22 * f 1e或 eaf22 BW(1e * sinaNW80椭球参数：长半轴 a=6378140 5（ m）. . .c. . .短半轴 b=6356755.2882m扁

3、率 =1/298.2573参心空间直角坐标转换参心坐标Larctan(Y )XB arctan(Z*(N H)( X 2)Y 2 ) * N * (1 e2 ) HX 2Y 2NHcosB二 高斯投影及高斯直角坐标系1、高斯投影概述高斯 - 克吕格投影的条件： 1. 是正形投影； 2. 中央子午线不变形高斯投影的性质： 1. 投影后角度不变； 2. 长度比与点位有关，与方向无关；3.离中央子午线越远变形越大为控制投影后的长度变形， 采用分带投影的方法。 常用 3 度带或 6 度带分带， 城市或工程控制网坐标可采用不按 3 度带中央子午线的任意带。2、高斯投影正算公式 ：xXN sin B co

4、s Bl 2N sin B cos3 B(5 t 29 24 4） l 4224N sin B cos5 B(6158t2t 4 )l 6720yN cos BlN cos3 B(1t 22 )l 36Ncos5 B(5 18t 2t 414258t 2 2 )l 51203、高斯投影反算公式：. . .c.1y1 (1 2t 2fy2l1f2 )cosB fN f6N f415 28t f224t f46 f28f2 t 2fy120N f2B B ft fy122922y2M fy15 3t ffft fN fN f12416190t 2f45t 4fy360N f1坐标转换简介坐标系统之

5、间的坐标转换既包括不同的参心坐标之间的转换， 或者不同的地心坐标系之间的转换，也包括参心坐标系与地心坐标系之间的转换以及相同坐标系的直角坐标与坐标之间的坐标转换， 还有坐标与高斯平面坐标之间的转换。 在两个空间角直坐标系中， 假设其分别为 O-XYZ 和 O-XYZ，如果两个坐标系的原点相同，通过三次旋转，就可以使两个坐标系重合；如果两个直角坐标系的原点不在同一个位置， 通过坐标轴的平移和旋转可以取得一致；如果两个坐标系的尺度也不尽一致， 就需要再增加一个尺度变化参数； 而对于坐标和高斯投影平面坐标之间的转换， 则需要通过高斯投影正算和高斯投影反算， 通过使用中央子午线的经度和不同的参考椭球以

6、及不同的投影面的选择来实现坐标的转换。如何使用 ArcGIS 实现 WGS84经纬度坐标到 BJ54 高斯投影坐标的转换？这是很多从事 GIS 工作或者测绘工作者普遍遇到的问题。本文目的在于帮助用户解决这个问题。我们通常说的 WGS-84坐标是指经纬度这种坐标表示方法， 54 坐标通常是指经过高斯投影的平面直角坐标这种坐标表示方法。 为什么要进行坐标转换？我们先来看两组参数， 如表1所示：表 1 BJ54 与 WGS84基准参数参考椭球体长半轴短半轴扁率BJ54 基准参数Krasovsky_194063782456356863.0188298.3WGS84 基准参数WGS 8463781376

7、356752.3142298.257224很显然， WGS84与 BJ54 是两种不同的基准面，不同的参考椭球体，因而两种地图下，同一个点的坐标是不同的，无论是三度带六度带坐标还是经纬度坐标都是不同的。当要把 GPS接收到的点（ WGS84坐标系统的）叠加到 BJ54 坐标系统的底图上，那就会发现这些 GPS点不能准确的在它该在的地方， 即“与实际地点发生了偏移” 。这就要求把这些 GPS点从 WGS84的坐标系统转换成 BJ54 的坐标系统了。有关 WGS84与 BJ54 的坐标转换问题， 实质是 WGS-84椭球体到 BJ54 椭球体的转换问题。如果我们是需要把 WGS84的经纬度坐标转换

8、成 BJ54 的高斯投影坐标，那就还会涉及到投影变换问题。因此，这个转换过程，一般的 GPS数据处理软件都是采用下述步骤进行的：1）（ B， L） 84（ X，Y， Z） 84，空间坐标到空间直角坐标的转换。2）（ X， Y， Z）84（ X， Y， Z） 54，坐标基准的转换，即 Datum 转换。通常有三种. . .c. . .转换方法：七参数、简化三参数、 Molodensky 。3）（ X， Y， Z）54（ B， L） 54，空间直角坐标到空间坐标的转换。4）（ B， L） 54（ x，y） 54， 高斯投影正算。从以上步骤不难看出， 转换的关键是第二步， 转换的参数。 鉴于我国曾使

9、用不同的坐标基准（ BJ54、State80 、Correct54 ），各地的重力值又有很大差异，所以很难确定一套适合全国且精度较好的转换参数。 在 WGS-84坐标和 54 坐标之间是不存在一套转换参数可以全国通用的，在每个地方会不一样。必须了解， 在不同的椭球之间的转换是不严密的。 那么， 两个椭球间的坐标转换应该是怎样的呢？一般而言比较严密的是用七参数法， 即 3 个平移因子 （ X 平移，Y 平移， Z 平移），3 个旋转因子（ X 旋转， Y 旋转， Z 旋转），一个比例因子（也叫尺度变化 K）。国参数来源的途径不多，一般当地测绘部门会有。通行的做法是：在工作区找三个以上的已知点，

10、利用已知点的 BJ54 坐标和所测 WGS84坐标，通过一定的数学模型，求解七参数。若多选几个已知点，通过平差的方法可以获得较好的精度。 如果区域围不大， 最远点间的距离不大于 30Km（经验值），这可以用三参数，即只考虑 3 个平移因子（ X 平移， Y 平移， Z 平移），而将旋转因子及比例因子（ X 旋转， Y 旋转， Z 旋转，尺度变化 K）都视为 0，所以三参数只是七参数的一种特例。 54 和 80 也是两种不同的基准面，不同的参考椭球体，他们之间的转换也是同理。 在 ArcGIS 中提供了三参数、 七参数转换法。 而在同一个椭球里的转换都是严密的，在同一个椭球的不同坐标系中转换需要

11、用到四参数转换，举个例子，在既有 54 坐标又有坐标，在这两种坐标之间转换就用到四参数，计算四参数需要两个已知点2 ArcGIS 坐标转换例子2.1 应注意问题使用 ArcGIS 如何实现 WGS84经纬度坐标到 BJ54 高斯投影坐标的转换呢？在 ArcGIS 中，这个坐标转换步骤简化了， 用户只需要两个步骤就能够直接从最初的 WGS84经纬度坐标转换到 BJ54 高斯投影坐标。这就是 ArcGIS 的强大之处。接下来， 我们做一个例子。 假设我们已经知道了 7 参数， 应该如何操作呢？在具体的操作前，请大家一定注意以下三点：WGS84的经纬度坐标值是用度来表示，而不能是度分秒表示七参数的平

12、移因子单位是米，旋转因子单位是秒，比例因子单位是百万。在 ArcGIS 中， 7 参数法的名字是 Coordinate_Frame 方法。有人在用 ArcGIS 进行不同椭球体间的坐标转换时，转换出来的结果不对，然后就写文章说变形如何如何，很可能是由于他们没有注意上面这三个关键的问题造成的。2.2 转换步骤a、定义 7 参数的地理转换（ Create Custom Geographic Transformation ）在 Arctool 中打开 Create Custom Geographic Transformation 工具，如图 1 所示：. . .c. . .在弹出的窗口中，输入一个转

13、换的名字，如 wgs84ToBJ54。在定义地理转换方法下面，在 Method 中选择合适的转换方法如 COORDINATE_FRAME，然后输入平移参数、旋转角度和比例因子，如图 2 所示：. . .c. . .b、投影变换打开工具箱下的 Projections and TransformationsFeatureProject ，在弹出的窗口中输入要转换的数据以及 Output Coordinate System ，然后输入第一步自定义的地理坐标系如 wgs84ToBJ54，开始投影变换，如图 3 所示：点击“确定”，完成坐标转换。3结束语我国现已启用新的坐标系统 2000 国家坐标系，

14、2000 国家坐标系与现行国家坐标系转换、衔接的过渡期仍需一段较长时期， 在实际工作、 工程中还遇到不同坐标系之间转换， 本文针对在生产中从事测绘工作遇到的坐标转换问题提供解决方法和经验，希望对同行有所参考。【参考文献】1孔祥元、郭际明、宗泉 . 测量学基础 . 大学，第一版， 2001 年 9 月2征航、黄劲松 . GPS测量与数据处理 . 大学，第一版， 2005 年 3 月3MAPGIS 使用教程坐标空间直角坐标转换(2009-10-22 21:12:41)标签：杂谈程序计算坐标与空间直角坐标转换#include #include . . .c. . .#define PI 3.1415

15、926535897932384626433832795double a,b,c,B,L,N,e,X,Y,Z,W,H;int choice;double B1=0.0,B2=0.0;double delta=0.0;int main()printf(please insert long r an);scanf(%f,&a);printf(please insert short r bn );scanf(%f,&b);e=sqrt(a*a-b*b)/a;c=a*a/b;printf(1kongzhitodadin2daditokongzhi0exit);scanf(%d,&choice);whil

16、e(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;X=(N+H)*cos(B)*cos(L);Y=(N+H)*cos(B)*sin(L);Z=(N*(1-e*e)+H)*sin(B);printf(X=%f,Y=%f,Z=%f,X,Y,Z);if(choice=1)printf(zuobiao:nX=);scanf(%f,&X);printf(zuobiao:n

17、Y=);scanf(%f,&Y);printf(zuobiao:nZ=);scanf(%f,&Z);L=atan(Y/X);. . .c. . ./double PP2=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);printf(L=%f,B=%f,H=%f);else printf(enter error,please enter again);pri

18、ntf(1kongzhitodadin2daditokongzhi0exit);scanf(%d,&choice);return 0;/高斯投影正、反算/高斯投影正、反算/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, X

19、0,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 = longitud

20、e0 * iPI ;. . .c. . .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/

21、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*l atitude1);xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/12 0);yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A

22、*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)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude)int ProjNo; int ZoneWide; /带宽double longit

23、ude1,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 = (P

24、rojNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;/中央经线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);. . .c. . .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)+(

25、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;/计算经度 (Longitu

26、de) 纬度 (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; /54 年坐标系参数长半轴 a=6378140.0; 扁率 f=1/298.257; /80 年坐标系参数长半轴 a=6378137m；扁率 f=1:298.257223563。/WGS-84 坐标系. . .c