空间直角坐标系与空间大地坐标系的相互转换及其C++源程序
- 格式:doc
- 大小:148.00 KB
- 文档页数:10
大地坐标转空间直角坐标方法1.准备工作:在进行大地坐标转换之前,首先要明确所采用的基准椭球参数,并且将大地坐标系转换为所采用的基准椭球上的坐标。
通常采用的基准椭球有WGS84、北京54和CGCS2000等。
这些基准椭球都有自己的参数,如长半轴a、偏心率e等。
根据所采用的基准椭球的参数,可以计算出该基准椭球的第一偏心率的平方(e^2)和扁率(f)等重要参数。
2.大地坐标转换为大地球面坐标:大地坐标的表示方法通常为经纬度(经度、纬度和高程)。
将经度和纬度转换为弧度形式,通过正弦定理和余弦定理等基本几何关系,可以计算出大地坐标在基准椭球上的投影参数。
利用这些参数,可以将大地坐标转换为大地球面坐标。
3.大地球面坐标转换为空间直角坐标:大地球面坐标是指基于基准椭球的坐标系,它只考虑地球的曲率而不考虑地球的引力场。
为了将其转换为直角坐标系,需要引入地球的引力场因素。
一种常见的方法是采用摄动参数法。
摄动参数法是通过导引纬度和经度等参数,计算出地球的重力梯度和坐标变换矩阵,并利用这些参数将大地球面坐标转换为空间直角坐标。
4.空间直角坐标的后处理:在将大地坐标转换为空间直角坐标之后,还需要进行一些后处理工作,以满足具体应用的要求。
例如,需要确定一个局部坐标系的原点和方向,进行坐标轴旋转和缩放等操作。
这些后处理工作可以在计算中进行,也可以在实际应用中进行。
总结起来,大地坐标到空间直角坐标的转换过程包括准备工作、大地坐标转大地球面坐标和大地球面坐标转空间直角坐标三个步骤。
在每个步骤中,需要根据具体问题选择合适的算法和参数。
同时,还需要注意坐标系之间的转换精度和误差控制,以确保转换结果的准确性。
坐标转换源代码--GPS定位程序(C++)GPS数据处理中为了满足不同的需要,处理的数据要进行坐标转换,得到在不同坐标系统下的结果,下面是笛卡尔坐标系,大地坐标系,站心地平坐标系(线型和极坐标形式)之间的转换源代码:头文件:#ifndef _COORDCOVERT_H#define _COORDCOVERT_H#include "stdlib.h"//WGS-84椭球体参数const double a=6378137.0;//长半轴const double flattening=1/298.257223563;//扁率const double delta=0.0000001;typedef struct tagCRDCARTESIAN{double x;double y;double z;}CRDCARTESIAN;typedef CRDCARTESIAN *PCRDCARTESIAN;//笛卡尔坐标系typedef struct tagCRDGEODETIC{double longitude;double latitude;double height;}CRDGEODETIC;typedef CRDGEODETIC *PCRDGEODETIC;//大地坐标系typedef struct tagCRDTOPOCENTRIC{double northing;double easting;double upping;}CRDTOPOCENTRIC;typedef CRDTOPOCENTRIC *PCRDTOPOCENTRIC; //站心地平坐标系(线坐标形式)typedef struct tagCRDTOPOCENTRICPOLAR{ double range;double azimuth;double elevation;}CRDTOPOCENTRICPOLAR;typedef CRDTOPOCENTRICPOLAR *PCRDTOPOCENTRICPOLAR;//站心地平坐标系(极坐标形式)//由笛卡尔坐标转换为大地坐标void CartesianToGeodetic (PCRDGEODETIC pcg, PCRDCARTESIAN pcc, double dSemiMajorAxis, double dFlattening);//pcg:指向所转换出的大地坐标的指针;//pcc:指向待转换的笛卡尔坐标的指针;//dSemiMajorAxis:参考椭球的长半轴;//dFlattening:参考椭球的扁率。
7.5 常用坐标系之间的关系与转换一、大地坐标系和空间大地直角坐标系及其关系大地坐标系用大地纬度企丈地经度L 和丈地髙H 来表示点的位置°这种坐标系是经 典大地测量甬:両用座标紊7屜据地图投影的理论,大地坐标系可以通过一定的投影转 化为投影平面上的直角坐标系,为地形测图和工程测量提供控制基础。
同时,这种坐标系 还是研究地球形状和大小的 种有用坐标系°所以大地坐标系在大地测量中始终有着重要 的作用.空间大地直角坐标系是-种以地球质心为原点购亘墮®坐标系,一般用X 、化Z 表 示点BSSTSTT 逐碇SS 範菇飞両H 绕禎扭转冻其轨道平面随时通过 地球质心。
对它们的跟踪观测也以地球质心为坐标原点,所以空间大地直角坐标系是卫星 大地测量中一种常用的基本坐标系。
现今,利用卫星大地测量的手段*可以迅速地测定点的空间大地直角坐拯,广泛应用于导航定位等空间技术。
同时经过数学变换,还可求岀点 的大地坐标I 用以加强和扩展地面大地网,进行岛屿和洲际联测,使传统的大地测量方法 发生了深刻的变化,所以空间大地宜角坐标系对现今大地测量的发展’具有重要的意义。
、大地坐标系和空间大地直角坐标系的转换如图7- 23所示’尸点的位置用空间 大地直角坐标〔X, Y, Z)表示,其相应 的大地坐标为(E, L)a 将该图与图?一5上式表明了 2种基本坐标系之间的关系。
加以比较可见,图7-5中的子午椭圆平面 相当于图7-23中的OJVP 平面.其中 PPz=Z.相当于图7-5中的j7;OP 3相当 丫于图7-5中的仏两平面的经度乙可视为相同,等于"叽 于是可以直接写岀X=jrcQsi f Y=jrsinL, Z=y将式(7-21).式(7-20)分别代入上式, 井考虑式(7-26)得X=Ncos^cosZr ”Y =NcQsBsinL > (7—78)Z=N (1—护〉sin^ ;BB 7-231.由大地坐标求空间大地直角坐标当已知椭球面上任一点P 的大地坐标(B, L)时,可以按式(7-78)直接求该点的 空间大地直角坐标(X, Y, Z)。
空间直角坐标系与大地坐标系转换程序.doc本文将介绍一种实现空间直角坐标系与大地坐标系转换的程序实现方法。
在编写程序时,需要使用一些数学库和函数,比如C++标准库中的cmath和iostream 等。
首先,我们需要了解空间直角坐标系和大明坐标系之间的转换公式。
假设空间直角坐标系为(x, y, z),大地坐标系为(L, B, H),则它们之间的转换公式为:x = cosLcosBsinHy = cosLsinBsinHz = sinLsinH其中,L为经度,B为纬度,H为高程。
根据上述公式,我们可以编写一个C++程序来实现空间直角坐标系与大地坐标系之间的转换。
程序实现如下:#include <iostream>#include <cmath>using namespace std;void transform() {double x, y, z;double L, B, H;cout << "Enter x, y, and z coordinates: ";cin >> x >> y >> z;cout << "Enter L and B coordinates: ";cin >> L >> B;H = acos(z / sqrt(x * x + y * y + z * z));cout << "The converted coordinates are: " << x << " " << y << " " << H << endl;}int main() {transform();return 0;}在上述程序中,我们首先定义了变量x、y、z、L、B和H,分别代表空间直角坐标系和大明坐标系的坐标值。
空间直角坐标系与空间大地坐标系的相互转换1.空间直角坐标系/笛卡尔坐标系坐标轴相互正交的坐标系被称作笛卡尔坐标系。
三维笛卡尔坐标系也被称为空间直角坐标系。
在空间直角坐标系下,点的坐标可以用该点所对应的矢径在三个坐标轴上的投影长度来表示,只有确定了原地、三个坐标轴的指向和尺度,就定义了一个在三维空间描述点的位置的空间直角坐标系。
以椭球体中心O为原点,起始子午面与赤道面交线为X轴,在赤道面上与X轴正交的方向为Y轴,椭球体的旋转轴为Z轴构成右手坐标系O.XYZ,在该坐标系中,P点的位置用X,Y,Z表示。
在测量应用中,常将地球空间直角坐标系的坐标原点选在地球质心(地心坐标系)或参考椭球中心(参心坐标系),z轴指向地球北极,x轴指向起始子午面与地球赤道的交点,y轴垂直于XOZ面并构成右手坐标系。
空间直角坐标系2.空间大地坐标系由于空间直角坐标无法明确反映出点与地球之间的空间关系,为了解决这一问题,在测量中引入了大地基准,并据此定义了大地坐标系。
大地基准指的是用于定义地球参考椭球的一系列参数,包括如下常量:2.1椭球的大小和形状2.2椭球的短半轴的指向:通常与地球的平自转轴平息。
2.3椭球中心的位置:根据需要确定。
若为地心椭球,则其中心位于地球质心。
2.4本初子午线:通过固定平极和经度原点的天文子午线,通常为格林尼治子午线。
以大地基准为基础建立的坐标系被称为大地坐标系。
由于大地基准又以参考椭球为基准,因此,大地坐标系又被称为椭球坐标系。
大地坐标系是参心坐标系,其坐标原点位于参考椭球中心,以参考椭球面为基准面,用大地经度L、纬度B 和大地高H表示地面点位置。
过地面点P的子午面与起始子午面间的夹角叫P 点的大地经度。
由起始子午面起算,向东为正,叫东经(0°~180°),向西为负,叫西经(0°~-180°)。
过P点的椭球法线与赤道面的夹角叫P点的大地纬度。
由赤道面起算,向北为正,叫北纬(0°~90°),向南为负,叫南纬(0°~-90°)。
#include "stdio、h"#include "math、h"#include "stdlib、h"#include "iostream"#define PI 3、14323double a,b,c,e2,ep2;int main(){int m,n,t;double RAD(double d,double f,double m);void RBD(double hd);void BLH_XYZ();void XYZ_BLH();void B_ZS();void B_FS();void GUS_ZS();void GUS_FS();printf(" 大地测量学\n");sp1:printf("请选择功能:\n");printf("1、大地坐标系到大地空间直角坐标的转换\n");printf("2、大地空间直角坐标到大地坐标系的转换\n");printf("3、贝塞尔大地问题正算\n");printf("4、贝塞尔大地问题反算\n");printf("5、高斯投影正算\n");printf("6、高斯投影反算\n");printf("0、退出程序\n");scanf("%d",&m);if(m==0)exit(0);sp2:printf("请选择椭球参数(输入椭球序号):\n");printf("1、克拉索夫斯基椭球参数\n");printf("2、IUGG_1975椭球参数\n");printf("3、CGCS_2000椭球参数\n");printf("0、其她椭球参数(自行输入)\n");scanf("%d",&n);switch(n){case 1:a=6378245、0;b=6356863、0188;c=6399698、9018;e2=0、297;ep2=0、468;break;case 2:a=6378140、0;b=6356755、2882;c=6399596、6520;e2=0、959;ep2=0、947;break;case 3:a=6378137、0;b=6356752、3141;c=6399593、6259;e2=0、290;ep2=0、547;break;case 0:{printf("请输入椭球参数:\n");printf("长半径a=");scanf("%lf",&a);printf("短半径b=");scanf("%lf",&b);c=a*a/b;ep2=(a*a-b*b)/(b*b);e2=(a*a-b*b)/(a*a);break;}default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp2 ;}while(1){switch(m){case 1:BLH_XYZ();break;case 2:XYZ_BLH();break;case 3:B_ZS();break;case 4:B_FS();break;case 5:GUS_ZS();break;case 6:GUS_FS();break;default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp1 ;}printf("就是否继续进行此功能计算? \n\n");printf("( 若继续进行此功能计算,则输入1;\n 若选择其她功能进行计算,则输入2;\n 若退出, 则输入0、)\n");scanf("%d",&t);switch(t){case 1:break;case 2:goto sp1;case 0:exit(0);}}}double RAD(double d,double f,double m) {double e;double sign=(d<0、0)?-1、0:1、0;if(d==0){sign=(f<0、0)?-1、0:1、0;if(f==0){sign=(m<0、0)?-1、0:1、0;}}if(d<0)d=d*(-1、0);if(f<0)f=f*(-1、0);if(m<0)m=m*(-1、0);e=sign*(d*3600+f*60+m)*PI/(3600*180);return e;}void RBD(double hd){int t;int d,f;double m;double sign=(hd<0、0)?-1、0:1、0;if(hd<0)hd=fabs(hd);hd=hd*3600*180/PI;t=int(hd/3600);d=sign*t;hd=hd-t*3600;f=int(hd/60);m=hd-f*60;printf("%d'%d'%lf'\n",d,f,m);}void BLH_XYZ(){double B,L,H,N,W;double d,f,m;double X,Y,Z;printf(" 请输入大地坐标(输入格式为角度(例如:30'40'50')):\n");printf(" 大地经度L=");scanf("%lf'%lf'%lf'",&d,&f,&m);L=RAD(d,f,m);printf(" 大地纬度B=");scanf("%lf'%lf'%lf'",&d,&f,&m);B=RAD(d,f,m);printf(" 大地高H=");scanf("%lf",&H);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;X=(N+H)*cos(B)*cos(L);Y=(N+H)*cos(B)*sin(L);Z=(N*(1-e2)+H)*sin(B);printf("\n\n 转换后得到大地空间直角坐标为:\n\n");printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);}void XYZ_BLH(){double B,L,H,N,W;double X,Y,Z;double tgB0,tgB1;printf(" 请输入大地空间直角坐标:\n");printf(" X=");scanf("%lf",&X);printf(" Y=");scanf("%lf",&Y);printf(" Z=");scanf("%lf",&Z);printf("\n\n 转换后得到大地坐标为:\n\n");L=atan(Y/X);printf(" 大地经度为: L=");RBD(L);printf("\n");tgB0=Z/sqrt(X*X+Y*Y);tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));while(fabs(tgB0-tgB1)>5*pow(10,-10)){tgB0=tgB1;tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));}B=atan(tgB1);printf(" 大地纬度为:B=");RBD(B);printf("\n");W=sqrt(1-e2*sin(B)*sin(B));N=a/W;H=sqrt(X*X+Y*Y)/cos(B)-N;printf(" 大地高为:H=%lf\n\n",H);}void B_ZS(){double L1,B1,A1,s,d,f,mi;double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;printf("请输入已知点的大地坐标(输入格式为角度(例如:30'40'50'),下同):\nL1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);L1=RAD(d,f,mi);printf("\nB1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);B1=RAD(d,f,mi);printf("请输入大地方位角:\nA1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);A1=RAD(d,f,mi);printf("请输入该点至另一点的大地线长:\ns=");scanf("%lf",&s);u1=atan(sqrt(1-e2)*tan(B1));m=asin(cos(u1)*sin(A1));M=atan(tan(u1)/cos(A1));m=(m>0)?m:m+2*PI;M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;sgm0=alfa*s;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); while(fabs(sgm0-sgm1)>2、8*PI/180*pow(10,-7)){sgm0=sgm1;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0); }sgm0=sgm1;A2=atan(tan(m)/cos(M+sgm0));A2=(A2>0)?A2:A2+PI;A2=(A1>PI)?A2:A2+PI;u2=atan(-cos(A2)*tan(M+sgm0));lmd1=atan(sin(u1)*tan(A1));lmd1=(lmd1>0)?lmd1:lmd1+PI;lmd1=(m<PI)?lmd1:lmd1+PI;lmd2=atan(sin(u2)*tan(A2));lmd2=(lmd2>0)?lmd2:lmd2+PI;lmd2=(m<PI)?(((M+sgm0)<PI)?lmd2:lmd2+PI):(((M+sgm0)>PI)?lmd2:lmd2+PI);lmd=lmd2-lmd1;B2=atan(sqrt(1+ep2)*tan(u2));kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;l=lmd-sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)+rp*sin(2*sgm0)*cos(4*M+2*sg m0));L2=L1+l;printf("\n\n得到另一点的大地坐标与大地线在该点的大地方位角为:\n\n");printf("L2=");RBD(L2);printf("\n");printf("B2=");RBD(B2);printf("\n");printf("A2=");RBD(A2);printf("\n");}void B_FS(){double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;printf("请输入第一个点大地坐标(输入格式为角度(例如:30'40'50'),下同):\n大地经度L1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L1=RAD(du,f,mi);printf("大地纬度B1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B1=RAD(du,f,mi);printf("\n请输入第二个点大地坐标:\n大地经度:L2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L2=RAD(du,f,mi);printf("大地纬度:B2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B2=RAD(du,f,mi);l=L2-L1;u1=atan(sqrt(1-e2)*tan(B1));u2=atan(sqrt(1-e2)*tan(B2));sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));m0=asin(cos(u1)*cos(u2)*sin(l)/sin(sgm0));dit_lmd=0、003351831*sgm0*sin(m0);lmd0=l+dit_lmd;dit_sgm=sin(m0)*dit_lmd;sgm1=sgm0+dit_sgm;m=asin(cos(u1)*cos(u2)*sin(lmd0)/sin(sgm1));A1=atan(sin(lmd0)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd0)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;M=atan(sin(u1)*tan(A1)/sin(m));M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));while(fabs(sgm0-sgm1)>1*PI/180*pow(10,-8)){sgm0=sgm1;sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));}sgm=sgm1;lmd=l+sin(m)*(alfap*sgm+btp*sin(sgm)*cos(2*M+sgm));s=(sgm-bt*sin(sgm)*cos(2*M+sgm)-r*sin(2*sgm)*cos(4*M+2*sgm))/alfa;A1=atan(sin(lmd)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2)));A2=(A2>0)?A2:A2+PI;A2=(m<0)?A2:A2+PI;printf("\n\n得到两点间大地线长S与大地正反方位角A1、A2如下:\n\n");printf("s=%lf\n",s);printf("A1=");RBD(A1);printf("\n");printf("A2=");RBD(A2);printf("\n");}void GUS_ZS(){double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;int DH3,DH6;printf("请输入大地坐标(输入格式为角度(例如:30'40'50')):\n大地经度L=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L=RAD(du,f,mi);printf("\n大地纬度B=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B=RAD(du,f,mi);At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;X=a*(1-e2)*(At*B-Bt*sin(2*B)/2+Ct*sin(4*B)/4-Dt*sin(6*B)/6);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;n=sqrt(ep2)*cos(B);t=tan(B);DH3=(L-(1、5*PI/180))/(3*PI/180)+1;DH6=L/(6*PI/180)+1;L03=DH3*(3*PI/180);L06=DH6*(6*PI/180)-(3*PI/180);l3=L-L03;l6=L-L06;m3=cos(B)*l3;m6=cos(B)*l6;x3=X+N*t*(m3*m3/2+(5-t*t+9*n*n+4*n*n*n*n)*m3*m3*m3*m3/24+(61-58*t*t+t*t*t*t)* m3*m3*m3*m3*m3*m3/720);x6=X+N*t*(m6*m6/2+(5-t*t+9*n*n+4*n*n*n*n)*m6*m6*m6*m6/24+(61-58*t*t+t*t*t*t)* m6*m6*m6*m6*m6*m6/720);y3=N*(m3+(1-t*t+n*n)*m3*m3*m3/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m3*m3*m3 *m3*m3/120);y6=N*(m6+(1-t*t+n*n)*m6*m6*m6/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m6*m6*m6 *m6*m6/120);Y3=DH3*1000000+500000+y3;Y6=DH6*1000000+500000+y6;printf("\n\n 得到的高斯平面坐标为:\n\n");printf(" 对于3度带:\n 纵坐标x=%、3lf\n 横坐标y=%、3lf(通用坐标Y=%、3lf)\n\n",x3,y3,Y3);printf(" 对于6度带:\n 纵坐标x=%、3lf\n 横坐标y=%、3lf(通用坐标Y=%、3lf)\n\n",x6,y6,Y6);}void GUS_FS(){double x,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;long DH;printf(" 请输入高斯平面坐标:\n\n");printf(" 纵坐标X=");scanf("%lf",&x);printf("\n");printf(" 自然坐标y=");scanf("%lf",&y);printf("\n");printf(" 通用坐标Y=");scanf("%lf",&Y);printf("\n");At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;B0=x/(a*(1-e2)*At);B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);while(fabs(B1-B0)>1*pow(10,-8)){B0=B1;B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);}Bf=B1;nf=sqrt(ep2)*cos(Bf);tf=tan(Bf);Vf=sqrt(1+ep2*cos(Bf)*cos(Bf));Nf=c/Vf;B=Bf-Vf*Vf*tf/2*((y/Nf)*(y/Nf)-(5+3*tf*tf+nf*nf-9*nf*nf*tf*tf)*pow((y/Nf),4)/12+(61+90*tf *tf+45*tf*tf)*pow((y/Nf),6)/360);L=((y/Nf)-(1+2*tf*tf+nf*nf)*(y/Nf)*(y/Nf)*(y/Nf)/6+(5+28*tf*tf+24*pow(tf,4)+6*nf*nf+8*nf *nf*tf*tf)*pow((y/Nf),5)/120)/cos(Bf);DH=Y/1000000;L3=3*PI/180*double(DH)+L;L6=6*PI/180*double(DH)-3*PI/180+L;printf("\n\n 得到的大地坐标为:\n\n");printf(" 大地纬度B=");RBD(B);printf("\n");printf(" 若为6度带,大地经度L=");RBD(L6);printf("\n");printf(" 若为3度带,大地经度L=");RBD(L3);printf("\n"); }。
大地坐标与直角空间坐标转换计算公式一、参心大地坐标与参心空间直角坐标转换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 ; d) 大地高H :地面点沿椭球法线至椭球面的距离为大地高H ; e) 地面点的点位用(B ,L ,H )表示。
2 参心大地坐标转换为参心空间直角坐标:⎪⎭⎪⎬⎫+-=+=+=B H e N Z L B H N Y L B H N X sin *])1(*[sin *cos *)(cos *cos *)(2 公式中,N 为椭球面卯酉圈的曲率半径,e 为椭球的第一偏心率,a 、b 椭球的长短半径,f 椭球扁率,W 为第一辅助系数a b a e 22-=或 f f e 1*2-= WaN B W e =-=22sin *1(西安80椭球参数:长半轴a=6378140±5(m )短半轴b=6356755.2882m 扁 率α=1/298.2573 参心空间直角坐标转换参心大地坐标[]NBY X H He N Y X H N Z B X YL -+=+-++==cos ))1(**)()(*arctan()arctan(22222二 高斯投影及高斯直角坐标系1、高斯投影概述高斯-克吕格投影的条件:1. 是正形投影;2. 中央子午线不变形高斯投影的性质:1. 投影后角度不变;2. 长度比与点位有关,与方向无关; 3. 离中央子午线越远变形越大为控制投影后的长度变形,采用分带投影的方法。
大地坐标系与空间直角坐标系的相互转换Python在地理信息系统(GIS)中,常常需要将大地坐标系(地理坐标系)与空间直角坐标系(笛卡尔坐标系)进行相互转换。
大地坐标系使用经纬度来表示地球表面上的任意点,而空间直角坐标系使用直角坐标来表示点在三维空间中的位置。
Python提供了一些库和工具,可以方便地进行这种转换。
大地坐标系与空间直角坐标系的基本概念大地坐标系(地理坐标系)大地坐标系是一种用经纬度来表示地球表面上任意点的坐标系。
经度表示点相对于本初子午线的位置(东经为正、西经为负),纬度表示点相对于赤道的位置(北纬为正、南纬为负)。
空间直角坐标系(笛卡尔坐标系)空间直角坐标系是一种使用直角坐标来表示点在三维空间中的位置的坐标系。
在空间直角坐标系中,每个点的位置由其相对于三个互相垂直的坐标轴的坐标值确定。
大地坐标系与空间直角坐标系的转换大地坐标系与空间直角坐标系之间的转换涉及到各种地球椭球参数和数学公式。
幸运的是,Python的一些库和工具已经实现了这些转换,使得我们可以很方便地进行转换操作。
Geopy库Geopy是一个Python库,提供了许多地理坐标系之间相互转换的功能。
使用Geopy,我们可以方便地进行大地坐标系到空间直角坐标系的转换。
首先,我们需要安装Geopy库。
可以使用pip命令来进行安装:pip install geopy接着,我们可以使用以下代码将大地坐标系的经纬度转换为空间直角坐标系的三维坐标:```python from geopy import Point from geopy.distance import distance定义大地坐标系的经纬度latitude = 40.7128 longitude = -74.0060将经纬度转换为空间直角坐标系的三维坐标point = Point(latitude, longitude) x, y, z = point.to_cartesian() print(f。
7.5 常用坐标系之间的关系与转换一、大地坐标系和空间大地直角坐标系及其关系 大地坐标系用大地纬度企丈地经度L 和丈地髙H 来表示点的位置°这种坐标系是经 典大地测量甬:両用座标紊7屜据地图投影的理论,大地坐标系可以通过一定的投影转 化为投影平面上的直角坐标系,为地形测图和工程测量提供控制基础。
同时,这种坐标系 还是研究地球形状和大小的 种有用坐标系°所以大地坐标系在大地测量中始终有着重要 的作用.空间大地直角坐标系是-种以地球质心为原点购亘墮®坐标系,一般用X 、化Z 表 示点BSSTSTT 逐碇SS 範菇飞両H 绕禎扭转冻其轨道平面随时通过 地球质心。
对它们的跟踪观测也以地球质心为坐标原点,所以空间大地直角坐标系是卫星 大地测量中一种常用的基本坐标系。
现今,利用卫星大地测量的手段*可以迅速地测定点的空间大地直角坐拯,广泛应用于导航定位等空间技术。
同时经过数学变换,还可求岀点 的大地坐标I 用以加强和扩展地面大地网,进行岛屿和洲际联测,使传统的大地测量方法 发生了深刻的变化,所以空间大地宜角坐标系对现今大地测量的发展’具有重要的意义。
、大地坐标系和空间大地直角坐标系的转换如图7- 23所示’尸点的位置用空间 大地直角坐标〔X, Y, Z)表示,其相应 的大地坐标为(E, L)a 将该图与图?一5加以比较可见,图7-5中的子午椭圆平面 相当于图7-23中的OJVP 平面.其中 PPz=Z.相当于图7-5中的j7;OP 3相当 丫于图7-5中的仏两平面的经度乙可视为相同,等于"叽 于是可以直接写岀X=jrcQsi f Y=jrsinL, Z=y将式(7-21).式(7-20)分别代入上式, 井考虑式(7-26)得X=Ncos^cosZr ”Y =NcQsBsinL > (7—78)Z=N (1—护〉sin^ ;上式表明了 2种基本坐标系之间的关系。
BB 7-231.由大地坐标求空间大地直角坐标当已知椭球面上任一点P 的大地坐标(B, L)时,可以按式(7-78)直接求该点的 空间大地直角坐标(X, Y, Z)。
《程序设计语言(C)》大作业报告题目:坐标系的转换完成人:小组构成及分工:*******独自完成程序的书写及调试.问题定义: 大地坐标和空间直角坐标系以及其他坐标系之间转换在卫星大地测量中经常用到的坐标系有,空间直角坐标系和大地直角坐标。
为了实现测量数据的快速高效的在不同的坐标系的转换,方便在学习及应用的中。
需要编写一程序实现数据的转换,实现空间直角坐标系与大地直角坐标之间在同一个系统中转换。
开发工具:Visual C++ 6.0数据结构描述:用不同的变量表示不同的坐标,变量选择时根据使用的习惯方便使用者的识别。
X:表示大地直角坐标的纵坐标;Y:表示大地直角坐标的横坐标;Z表示大地直角坐标的竖坐标L:表示空间直角坐标的经度;B:表示空间直角坐标的纬度;H:表示空间直角坐标的高度;算法描述:通过编写一个主函数描述出整个程序的主体不分,然后通过调用函数实现坐标的转换。
程序调试情况:坐标由大地直角坐标系中的转换为空间直角坐标系的坐标:大地直角坐标转换后空间直角坐标:B=60; X=2055059.130122;L=50; Y=2449123.986892;H=100; Z=5500477.615329;坐标由空间直角坐标系中的坐标转换为大地直角坐标系中的坐标;空间直角坐标转换后大地直角坐标;X=100; B=-127.103127844;Y=100; L=45.000000000;Z=10000; H=-6391994.685276;参考文献或网站:1.《控制测量学》(下册)第三版孔祥元郭际明主编武汉大学出版社;2. 《数字测图原理与方法》第二版潘正风程效军成枢王腾军宋伟东邹进贵编著武汉大学出版社;3.《C 程序设计语言》魏东平朱连章于广斌编著;电子工业出版社心得体会:写这个大作业确实让我收获了许多!1.写这次计算机大作业,让我经历了一个难忘的过程。
自己的是必须得自己独立自主的想办法去解决,没人会为与自己没多大关系的事分很多神的!2.经历了过程,让我学到了些东西也在解决困难的过程中认识了些学长,他们也教会了我许多学习经验。
#include "stdio.h"#include "math.h"#include "stdlib.h"#include "iostream"#define PI 3.1415926535897323double a,b,c,e2,ep2;int main(){int m,n,t;double RAD(double d,double f,double m);void RBD(double hd);void BLH_XYZ();void XYZ_BLH();void B_ZS();void B_FS();void GUS_ZS();void GUS_FS();printf(" 大地测量学\n");sp1:printf("请选择功能:\n");printf("1.大地坐标系到大地空间直角坐标的转换\n");printf("2.大地空间直角坐标到大地坐标系的转换\n");printf("3.贝塞尔大地问题正算\n");printf("4.贝塞尔大地问题反算\n");printf("5.高斯投影正算\n");printf("6.高斯投影反算\n");printf("0.退出程序\n");scanf("%d",&m);if(m==0)exit(0);sp2:printf("请选择椭球参数(输入椭球序号):\n");printf("1.克拉索夫斯基椭球参数\n");printf("2.IUGG_1975椭球参数\n");printf("3.CGCS_2000椭球参数\n");printf("0.其他椭球参数(自行输入)\n");scanf("%d",&n);switch(n){case1:a=6378245.0;b=6356863.0188;c=6399698.9018;e2=0.00669342162297;ep2=0.0067385254146 8;break;case2:a=6378140.0;b=6356755.2882;c=6399596.6520;e2=0.00669438499959;ep2=0.0067395018194 7;break;case3:a=6378137.0;b=6356752.3141;c=6399593.6259;e2=0.00669438002290;ep2=0.0067394967754 7;break;case 0:{printf("请输入椭球参数:\n");printf("长半径a=");scanf("%lf",&a);printf("短半径b=");scanf("%lf",&b);c=a*a/b;ep2=(a*a-b*b)/(b*b);e2=(a*a-b*b)/(a*a);break;}default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp2 ;}while(1){switch(m){case 1:BLH_XYZ();break;case 2:XYZ_BLH();break;case 3:B_ZS();break;case 4:B_FS();break;case 5:GUS_ZS();break;case 6:GUS_FS();break;default:printf("\n\n输入错误!\n请重新输入!\n\n");goto sp1 ;}printf("是否继续进行此功能计算? \n\n");printf("(若继续进行此功能计算,则输入1;\n 若选择其他功能进行计算,则输入2;\n 若退出,则输入0. )\n");scanf("%d",&t);switch(t){case 1:break;case 2:goto sp1;case 0:exit(0);}}}double RAD(double d,double f,double m){double e;double sign=(d<0.0)?-1.0:1.0;if(d==0){sign=(f<0.0)?-1.0:1.0;if(f==0){sign=(m<0.0)?-1.0:1.0;}}if(d<0)d=d*(-1.0);if(f<0)f=f*(-1.0);if(m<0)m=m*(-1.0);e=sign*(d*3600+f*60+m)*PI/(3600*180);return e;}void RBD(double hd){int t;int d,f;double m;double sign=(hd<0.0)?-1.0:1.0;if(hd<0)hd=fabs(hd);hd=hd*3600*180/PI;t=int(hd/3600);d=sign*t;hd=hd-t*3600;f=int(hd/60);m=hd-f*60;printf("%d'%d'%lf'\n",d,f,m);}void BLH_XYZ(){double B,L,H,N,W;double d,f,m;double X,Y,Z;printf(" 请输入大地坐标(输入格式为角度(例如:30'40'50')):\n");printf(" 大地经度L=");scanf("%lf'%lf'%lf'",&d,&f,&m);L=RAD(d,f,m);printf(" 大地纬度B=");scanf("%lf'%lf'%lf'",&d,&f,&m);B=RAD(d,f,m);printf(" 大地高H=");scanf("%lf",&H);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;X=(N+H)*cos(B)*cos(L);Y=(N+H)*cos(B)*sin(L);Z=(N*(1-e2)+H)*sin(B);printf("\n\n 转换后得到大地空间直角坐标为:\n\n");printf("X=%lf\nY=%lf\nZ=%lf\n\n",X,Y,Z);}void XYZ_BLH(){double B,L,H,N,W;double X,Y,Z;double tgB0,tgB1;printf(" 请输入大地空间直角坐标:\n");printf(" X=");scanf("%lf",&X);printf(" Y=");scanf("%lf",&Y);printf(" Z=");scanf("%lf",&Z);printf("\n\n 转换后得到大地坐标为:\n\n");L=atan(Y/X);printf(" 大地经度为:L=");RBD(L);printf("\n");tgB0=Z/sqrt(X*X+Y*Y);tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));while(fabs(tgB0-tgB1)>5*pow(10,-10)){tgB0=tgB1;tgB1=(1/sqrt(X*X+Y*Y))*(Z+a*e2*tgB0/sqrt(1+tgB0*tgB0-e2*tgB0*tgB0));}B=atan(tgB1);printf(" 大地纬度为:B=");RBD(B);printf("\n");W=sqrt(1-e2*sin(B)*sin(B));N=a/W;H=sqrt(X*X+Y*Y)/cos(B)-N;printf(" 大地高为:H=%lf\n\n",H);}void B_ZS(){double L1,B1,A1,s,d,f,mi;double u1,u2,m,M,k2,alfa,bt,r,kp2,alfap,btp,rp;double sgm0,sgm1,lmd,lmd1,lmd2,A2,B2,l,L2;printf("请输入已知点的大地坐标(输入格式为角度(例如:30'40'50'),下同):\nL1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);L1=RAD(d,f,mi);printf("\nB1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);B1=RAD(d,f,mi);printf("请输入大地方位角:\nA1=");scanf("%lf'%lf'%lf'",&d,&f,&mi);A1=RAD(d,f,mi);printf("请输入该点至另一点的大地线长:\ns=");scanf("%lf",&s);u1=atan(sqrt(1-e2)*tan(B1));m=asin(cos(u1)*sin(A1));M=atan(tan(u1)/cos(A1));m=(m>0)?m:m+2*PI;M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;sgm0=alfa*s;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0);while(fabs(sgm0-sgm1)>2.8*PI/180*pow(10,-7)){sgm0=sgm1;sgm1=alfa*s+bt*sin(sgm0)*cos(2*M+sgm0)+r*sin(2*sgm0)*cos(4*M+2*sgm0);}sgm0=sgm1;A2=atan(tan(m)/cos(M+sgm0));A2=(A2>0)?A2:A2+PI;A2=(A1>PI)?A2:A2+PI;u2=atan(-cos(A2)*tan(M+sgm0));lmd1=atan(sin(u1)*tan(A1));lmd1=(lmd1>0)?lmd1:lmd1+PI;lmd1=(m<PI)?lmd1:lmd1+PI;lmd2=atan(sin(u2)*tan(A2));lmd2=(lmd2>0)?lmd2:lmd2+PI;lmd2=(m<PI)?(((M+sgm0)<PI)?lmd2:lmd2+PI):(((M+sgm0)>PI)?lmd2:lmd2+PI);lmd=lmd2-lmd1;B2=atan(sqrt(1+ep2)*tan(u2));kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;l=lmd-sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos(2*M+sgm0)+rp*sin(2*sgm0)*cos(4*M+2*sgm0));L2=L1+l;printf("\n\n得到另一点的大地坐标和大地线在该点的大地方位角为:\n\n");printf("L2=");RBD(L2);printf("\n");printf("B2=");RBD(B2);printf("\n");printf("A2=");RBD(A2);printf("\n");}void B_FS(){double L1,B1,L2,B2,s,A1,A2,du,f,mi,m0,m,M;double l,u1,u2,alfa,bt,r,lmd0,dit_lmd,lmd,sgm,dit_sgm,sgm0,sgm1,alfap,btp,rp,k2,kp2;printf("请输入第一个点大地坐标(输入格式为角度(例如:30'40'50'),下同):\n大地经度L1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L1=RAD(du,f,mi);printf("大地纬度B1=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B1=RAD(du,f,mi);printf("\n请输入第二个点大地坐标:\n大地经度:L2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L2=RAD(du,f,mi);printf("大地纬度:B2=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B2=RAD(du,f,mi);l=L2-L1;u1=atan(sqrt(1-e2)*tan(B1));u2=atan(sqrt(1-e2)*tan(B2));sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));m0=asin(cos(u1)*cos(u2)*sin(l)/sin(sgm0));dit_lmd=0.003351831*sgm0*sin(m0);lmd0=l+dit_lmd;dit_sgm=sin(m0)*dit_lmd;sgm1=sgm0+dit_sgm;m=asin(cos(u1)*cos(u2)*sin(lmd0)/sin(sgm1));A1=atan(sin(lmd0)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd0)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;M=atan(sin(u1)*tan(A1)/sin(m));M=(M>0)?M:M+PI;k2=ep2*cos(m)*cos(m);alfa=(1-k2/4+7*k2*k2/64-15*k2*k2*k2/256)/b;bt=k2/4-k2*k2/8+37*k2*k2*k2/512;r=k2*k2/128-k2*k2*k2/128;kp2=e2*cos(m)*cos(m);alfap=(e2/2+e2*e2/8+e2*e2*e2/16)-e2/16*(1+e2)*kp2+3*e2*kp2*kp2/128;btp=e2*(1+e2)*kp2/16-e2*kp2*kp2/32;rp=e2*kp2*kp2/256;sgm0=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l));sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));while(fabs(sgm0-sgm1)>1*PI/180*pow(10,-8)){sgm0=sgm1;sgm1=acos(sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(l+sin(m)*(alfap*sgm0+btp*sin(sgm0)*cos( 2*M+sgm0))));}sgm=sgm1;lmd=l+sin(m)*(alfap*sgm+btp*sin(sgm)*cos(2*M+sgm));s=(sgm-bt*sin(sgm)*cos(2*M+sgm)-r*sin(2*sgm)*cos(4*M+2*sgm))/alfa;A1=atan(sin(lmd)/(cos(u1)*tan(u2)-sin(u1)*cos(lmd)));A1=(A1>0)?A1:A1+PI;A1=(m>0)?A1:A1+PI;A2=atan(sin(lmd)/(sin(u2)*cos(lmd)-tan(u1)*cos(u2)));A2=(A2>0)?A2:A2+PI;A2=(m<0)?A2:A2+PI;printf("\n\n得到两点间大地线长S和大地正反方位角A1、A2如下:\n\n");printf("s=%lf\n",s);printf("A1=");RBD(A1);printf("\n");printf("A2=");RBD(A2);printf("\n");}void GUS_ZS(){double B,L,x3,x6,y3,y6,Y3,Y6,du,f,mi,X,N,n,t;double At,Bt,Ct,Dt,m3,m6,l3,l6,W,L03,L06;int DH3,DH6;printf("请输入大地坐标(输入格式为角度(例如:30'40'50')):\n大地经度L=");scanf("%lf'%lf'%lf'",&du,&f,&mi);L=RAD(du,f,mi);printf("\n大地纬度B=");scanf("%lf'%lf'%lf'",&du,&f,&mi);B=RAD(du,f,mi);At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;X=a*(1-e2)*(At*B-Bt*sin(2*B)/2+Ct*sin(4*B)/4-Dt*sin(6*B)/6);W=sqrt(1-e2*sin(B)*sin(B));N=a/W;n=sqrt(ep2)*cos(B);t=tan(B);DH3=(L-(1.5*PI/180))/(3*PI/180)+1;DH6=L/(6*PI/180)+1;L03=DH3*(3*PI/180);L06=DH6*(6*PI/180)-(3*PI/180);l3=L-L03;l6=L-L06;m3=cos(B)*l3;m6=cos(B)*l6;x3=X+N*t*(m3*m3/2+(5-t*t+9*n*n+4*n*n*n*n)*m3*m3*m3*m3/24+(61-58*t*t+t*t*t*t)*m3*m3*m3*m3*m3*m3/720);x6=X+N*t*(m6*m6/2+(5-t*t+9*n*n+4*n*n*n*n)*m6*m6*m6*m6/24+(61-58*t*t+t*t*t*t)*m6*m6*m6*m6*m6*m6/720);y3=N*(m3+(1-t*t+n*n)*m3*m3*m3/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m3*m3*m3*m3*m3/120);y6=N*(m6+(1-t*t+n*n)*m6*m6*m6/6+(5-18*t*t+t*t*t*t+14*n*n-58*n*n*t*t)*m6*m6*m6*m6*m6/120);Y3=DH3*1000000+500000+y3;Y6=DH6*1000000+500000+y6;printf("\n\n 得到的高斯平面坐标为:\n\n");printf(" 对于3度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x3,y3,Y3);printf(" 对于6度带:\n 纵坐标x=%.3lf\n 横坐标y=%.3lf(通用坐标Y=%.3lf)\n\n",x6,y6,Y6);}void GUS_FS(){double x,y,Y,B,B0,B1,Bf,Vf,tf,Nf,nf,L,At,Bt,Ct,Dt,L3,L6;long DH;printf(" 请输入高斯平面坐标:\n\n");printf(" 纵坐标X=");scanf("%lf",&x);printf("\n");printf(" 自然坐标y=");scanf("%lf",&y);printf("\n");printf(" 通用坐标Y=");scanf("%lf",&Y);printf("\n");At=1+3*e2/4+45*e2*e2/64+175*e2*e2*e2/256;Bt=3*e2/4+15*e2*e2/16+525*e2*e2*e2/512;Ct=15*e2*e2/64+105*e2*e2*e2/256;Dt=35*e2*e2*e2/512;B0=x/(a*(1-e2)*At);B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);while(fabs(B1-B0)>1*pow(10,-8)){B0=B1;B1=(x-a*(1-e2)*(-Bt*sin(2*B0)/2+Ct*sin(4*B0)/4-Dt*sin(6*B0)/6))/(a*(1-e2)*At);}Bf=B1;nf=sqrt(ep2)*cos(Bf);tf=tan(Bf);Vf=sqrt(1+ep2*cos(Bf)*cos(Bf));Nf=c/Vf;B=Bf-Vf*Vf*tf/2*((y/Nf)*(y/Nf)-(5+3*tf*tf+nf*nf-9*nf*nf*tf*tf)*pow((y/Nf),4)/12+(61+90*tf*tf+45*tf*tf)*pow((y/Nf),6)/360);L=((y/Nf)-(1+2*tf*tf+nf*nf)*(y/Nf)*(y/Nf)*(y/Nf)/6+(5+28*tf*tf+24*pow(tf,4)+6*nf*nf+8*nf*nf*tf*tf)*po w((y/Nf),5)/120)/cos(Bf);DH=Y/1000000;L3=3*PI/180*double(DH)+L;L6=6*PI/180*double(DH)-3*PI/180+L;printf("\n\n 得到的大地坐标为:\n\n");printf(" 大地纬度B=");RBD(B);printf("\n");printf(" 若为6度带,大地经度L=");RBD(L6);printf("\n");printf(" 若为3度带,大地经度L=");RBD(L3);printf("\n");}。
存档日期:存档编号:本科生毕业论文(设计)论文题目:空间直角坐标系统转换及程序设计姓名:学号:院系:专业:测绘工程随着空间技术的发展,全球一体化的形成,越来越多的要求全球测绘资料形成统一规范,尤其是坐标系统的统一。
由于各测量单位工作目的不同,所选择的椭球参考系也会有所不同,出现了许多不同形式的坐标系,例如WGS-84坐标系、国家80坐标系、北京54坐标系、独立地方坐标及各种坐标。
在同一坐标系下坐标的表示方式又有空间直角坐标、大地坐标、平面坐标。
根据不同的测绘需求,需要将不同的坐标系下的坐标进行相互转换,在这些坐标转换的过程中既会运用到同一坐标系下的坐标转换模型,又会用到不同参考系下各坐标系间的坐标转换模型。
首先本文介绍大地测量学坐标的相关知识,接着介绍了与坐标转换相关的知识以及坐标转换模型(包括同一坐标系下的坐标转换模型和不同参考系下各坐标系间的坐标转换模型),并利用vb语言实现坐标转换的过程。
关键词:地球椭球,坐标系,转换模型,坐标转换Along with the development of space technology, the formation of global integration, more and more requirements of surveying and mapping material form a unified global standard, especially the unity of the coordinate system. Because each measurement unit work purpose is different, choose the frame of reference ellipsoid would differ, the emergence of many different forms of coordinate system, such as WGS-84 coordinate system, the state 80 coordinate system, Beijing 54 coordinate system, independent local coordinate system and various kinds of urban construction coordinates. In the same coordinate system of representation and coordinate space right-angle coordinate, coordinate, coordinate the earth plane. According to the different needs of surveying and mapping, need different coordinate transformation coordinate system, in which the process of coordinate transformation can use to the same coordinate system coordinate transformation model, and will use different reference frame, the coordinate transformation between the coordinate system model.First of all this paper introduces the geodetic coordinates of relevant knowledge, then introduce the knowledge of coordinate transformation and coordinate transformation model(including the same coordinate system coordinate transformation model and different reference frame, the coordinate transformation between the coordinate system model), and the use of vb language realization of coordinate transformation process.Key words:the earth ellipsoid, coordinate system, transformation model, coordinate transformation摘要 (I)Abstract (II)1 绪论 (1)1.1研究的背景和意义 (1)1.2国内外研究现状 (1)1.3研究的主要内容 (2)2 相关理论和知识介绍 (3)2.1地球椭球 (3)2.2 基准 (4)2.3建立大地坐标系的基本原理 (5)2.3.1椭球定位、定向的概念 (5)2.3.2参考椭球定位与定向的实现方法 (5)2.3.3多点定位 (6)3 坐标系统简介 (7)3.1测量常用的坐标系 (7)3.1.1大地坐标系 (7)3.1.2空间直角坐标系 (8)3.1.3平面坐标系 (8)3.1.4地方独立坐标系 (9)3.2 我国常用的坐标系统 (9)3.2.1 1954年北京坐标系 (10)3.2.2 1980年国家大地坐标系 (11)3.2.3 1954年北京坐标系(整体平差转换值) (11)3.2.4 WGS-84世界大地坐标系 (12)3.2.5 2000国家大地坐标系 (12)4 坐标转换理论与程序设计 (14)4.1坐标转换的基本概念 (14)4.2坐标系转换的模型 (14)4.2.1同一参考椭球下大地坐标系与空间直角坐标系的相互转换 (14)4.3基准转换的模型 (21)4.3.1不同地球椭球坐标系的空间三参数或七参数转换 (21)5 全文总结 (26)参考文献 (27)附录 (28)致谢 (28)1 绪论1.1研究的背景和意义随着大地测量学,卫星大地测量学,摄影测量学的发展和电子计算机的普及,对各种坐标系的研究变得越来越重要了。
空间直角坐标系与空间大地坐标系的相互转换1.空间直角坐标系/笛卡尔坐标系坐标轴相互正交的坐标系被称作笛卡尔坐标系。
三维笛卡尔坐标系也被称为空间直角坐标系。
在空间直角坐标系下,点的坐标可以用该点所对应的矢径在三个坐标轴上的投影长度来表示,只有确定了原地、三个坐标轴的指向和尺度,就定义了一个在三维空间描述点的位置的空间直角坐标系。
以椭球体中心O为原点,起始子午面与赤道面交线为X轴,在赤道面上与X轴正交的方向为Y轴,椭球体的旋转轴为Z轴构成右手坐标系O.XYZ,在该坐标系中,P点的位置用X,Y,Z表示。
在测量应用中,常将地球空间直角坐标系的坐标原点选在地球质心(地心坐标系)或参考椭球中心(参心坐标系),z轴指向地球北极,x轴指向起始子午面与地球赤道的交点,y轴垂直于XOZ面并构成右手坐标系。
空间直角坐标系2.空间大地坐标系由于空间直角坐标无法明确反映出点与地球之间的空间关系,为了解决这一问题,在测量中引入了大地基准,并据此定义了大地坐标系。
大地基准指的是用于定义地球参考椭球的一系列参数,包括如下常量:2.1椭球的大小和形状2.2椭球的短半轴的指向:通常与地球的平自转轴平息。
2.3椭球中心的位置:根据需要确定。
若为地心椭球,则其中心位于地球质心。
2.4本初子午线:通过固定平极和经度原点的天文子午线,通常为格林尼治子午线。
以大地基准为基础建立的坐标系被称为大地坐标系。
由于大地基准又以参考椭球为基准,因此,大地坐标系又被称为椭球坐标系。
大地坐标系是参心坐标系,其坐标原点位于参考椭球中心,以参考椭球面为基准面,用大地经度L、纬度B 和大地高H表示地面点位置。
过地面点P的子午面与起始子午面间的夹角叫P 点的大地经度。
由起始子午面起算,向东为正,叫东经(0°~180°),向西为负,叫西经(0°~-180°)。
过P点的椭球法线与赤道面的夹角叫P点的大地纬度。
由赤道面起算,向北为正,叫北纬(0°~90°),向南为负,叫南纬(0°~-90°)。
空间直角坐标系与空间大地坐标系的相互转换1.空间直角坐标系/笛卡尔坐标系坐标轴相互正交的坐标系被称作笛卡尔坐标系。
三维笛卡尔坐标系也被称为空间直角坐标系。
在空间直角坐标系下,点的坐标可以用该点所对应的矢径在三个坐标轴上的投影长度来表示,只有确定了原地、三个坐标轴的指向和尺度,就定义了一个在三维空间描述点的位置的空间直角坐标系。
以椭球体中心O为原点,起始子午面与赤道面交线为X轴,在赤道面上与X轴正交的方向为Y轴,椭球体的旋转轴为Z轴构成右手坐标系O.XYZ,在该坐标系中,P点的位置用X,Y,Z表示。
在测量应用中,常将地球空间直角坐标系的坐标原点选在地球质心(地心坐标系)或参考椭球中心(参心坐标系),z轴指向地球北极,x轴指向起始子午面与地球赤道的交点,y轴垂直于XOZ面并构成右手坐标系。
空间直角坐标系2.空间大地坐标系由于空间直角坐标无法明确反映出点与地球之间的空间关系,为了解决这一问题,在测量中引入了大地基准,并据此定义了大地坐标系。
大地基准指的是用于定义地球参考椭球的一系列参数,包括如下常量:2.1椭球的大小和形状2.2椭球的短半轴的指向:通常与地球的平自转轴平息。
2.3椭球中心的位置:根据需要确定。
若为地心椭球,则其中心位于地球质心。
2.4本初子午线:通过固定平极和经度原点的天文子午线,通常为格林尼治子午线。
以大地基准为基础建立的坐标系被称为大地坐标系。
由于大地基准又以参考椭球为基准,因此,大地坐标系又被称为椭球坐标系。
大地坐标系是参心坐标系,其坐标原点位于参考椭球中心,以参考椭球面为基准面,用大地经度L、纬度B 和大地高H表示地面点位置。
过地面点P的子午面与起始子午面间的夹角叫P 点的大地经度。
由起始子午面起算,向东为正,叫东经(0°~180°),向西为负,叫西经(0°~-180°)。
过P点的椭球法线与赤道面的夹角叫P点的大地纬度。
由赤道面起算,向北为正,叫北纬(0°~90°),向南为负,叫南纬(0°~-90°)。
从地面点P沿椭球法线到椭球面的距离叫大地高。
大地坐标坐标系中,P点的位置用L,B表示。
如果点不在椭球面上,表示点的位置除L,B外,还要附加另一参数——大地高H。
空间大地坐标系3.空间直角坐标与大地坐标间的转换3.1大地坐标转换为空间直角坐标将同一坐标系下的大地坐标(B、L、H)转换成空间直角坐标(X、Y、Z)的转换公式为:式中N为卯酉圈的半径,a为参考椭球的长半轴;b为参考椭球的短半轴;e为参考椭球的第一偏心率;并且有若点在椭球面上,则大地高H=0,上式可简化为:3.2空间直角坐标转换为空间大地坐标将同一坐标系下的空间直角坐标(X、Y、Z)转换为空间大地坐标(B、L、H)的公式为:在使用上式进行空间直角坐标到大地坐标的转换过程中,由于计算大地纬度口时用到大地高Ⅳ,而计算大地高时又需要用到大地纬度口.因此不能直接由空间直角坐标计算出大地坐标,而需要采用迭代计算的方法。
具体计算时,可先根据下式求出大地纬度口的初值:然后利用该初值来求出H、N的初值,再利用所求出的H和N初值再次求出B 值.如此反复,直至求出的及日、Ⅳ收敛为止。
4.算例本文根据以上公式在Microsoft VC++6.0环境下编写了一段程序(见附录)。
算例中的坐标采用的是武汉大学信息学部友谊广场上的某点的大地坐标作为已知值,然后经过转换函数CRDGEODETICtoCRDCARTESEAN(pcg, pcc, dSemiMajorAxis,dFlatning)把大地坐标转换为空间直角坐标得到坐标X、Y、Z。
由得到的空间直角坐标X、Y、Z,经过转换函数CRDCARTESIANtoCRDGEODETIC(pcc,pcg,dSemiMajorAxis,dFlatning)把空间直角坐标还原成空间大地坐标,计算结果如下图所示:计算结果从上图可以看出结果比较满意,高程和精度基本能完全还原,而纬度还原后有较大的误差,在测量中这种误差不允许的,需要修改算法,完善结果。
可能引起的原因有可能是由于纬度计算公式并不完善,还有可能是由于计算机的截断误差引起的,还要找时间继续修改、完善。
5.心得体会这次编程自认为很简单,但真动手自己亲自编写,还是或多或少遇到了一些问题,并分析问题,最终解决问题。
虽然这次作业很简单,但经过自己这样一步一步的编写出来,还是有很多收获,加强了运用VC++编写程序的能力,也充分认识到了学习VC++的重要性,更找到了自己的一些缺点与不足。
6.附录(程序源代码)#include<iostream>#include<cmath>using namespace std;#define M_PI 3.1415926typedef struct tagCRDCARTESIAN{double x;double y;double z;}CRDCARTESIAN;typedef CRDCARTESIAN *PCRDCARTESIAN; typedef struct tagCRDGEODETIC{double longitude;double latitude;double height;}CRDGEODETIC;typedef CRDGEODETIC *PCRDGEODETIC;void DMS_RAD(double DMS,double *Rad){int Deg,Min;double Sec;Deg=(int)DMS;Min=(int)((DMS-Deg)*100);Sec=((DMS-Deg)*100-Min)*100;*Rad=(Deg+Min/60.0+Sec/3600.0)/180.0*M_PI;return;}void RAD_DMS(double Rad,double *DMS){int Deg,Min;double Sec;double AR,AM;AR=Rad;if (Rad<0)AR=-Rad;AR=AR+1.0e-10;AR=AR*180.0/M_PI;Deg=(int)AR;AM=(AR-Deg)*60.0;Min=(int)AM;Sec=(AM-Min)*60;*DMS=Deg+Min/100.0+Sec/10000.0;if(Rad<0)*DMS=-*DMS;return;}bool CRDCARTESIANtoCRDGEODETIC(PCRDCARTESIAN pcc,PCRDGEODETIC pcg,double dSemiMajorAxis,double dFlattening) {double B0,R,N;double B_,L_;double X=pcc->x;double Y=pcc->y;double Z=pcc->z;R=sqrt(X*X+Y*Y);B0=atan2(Z,R);while (1){N=dSemiMajorAxis/sqrt(1.0-dFlattening*(2-dFlattening)*sin(B0)*sin(B0));B_=atan2(Z+N*dFlattening*(2-dFlattening)*sin(B0),R);if(fabs(B_-B0)<1.0e-10)break;B0=B_;}L_=atan2(Y,X);pcg->height=R/cos(B_)-N;RAD_DMS(B_,&pcg->latitude);RAD_DMS(L_,&pcg->longitude);return true;}bool CRDGEODETICtoCRDCARTESEAN(PCRDGEODETIC pcg,PCRDCARTESIAN pcc,double dSemiMajorAxis,double dFlattening) {double N;double B_,L_;double B=pcg->latitude;double L=pcg->longitude;double H=pcg->height;DMS_RAD(B,&B_);DMS_RAD(L,&L_);N=dSemiMajorAxis/sqrt(1.0-dFlattening*(2-dFlattening)*sin(B_)*sin(B_));pcc->x=(N+H)*cos(B_)*cos(L_);pcc->y=(N+H)*cos(B_)*sin(L_);pcc->z=(N*(1.0-dFlattening*(2-dFlattening))+H)*sin(B_);return true;}void main(){PCRDCARTESIAN pcc=new CRDCARTESIAN;PCRDGEODETIC pcg=new CRDGEODETIC;//B=30.31.40.23 L=114.21.20.51 h=41double rad;rad=(30*3600+31*60+40.23)/3600;pcg->latitude=rad;rad=(114*3600+21*60+20.51)/3600;pcg->height=41;pcg->longitude=rad;double dSemiMajorAxis=6378137;double dFlatning=1/298.257223563;cout<<std::fixed;cout<<"转换前已知的大地坐标:"<<endl;cout<<"H(高度)="<<pcg->height<<" "<<"L(经度)="<<pcg->longitude<<" "<<"B(纬度)="<<pcg->latitude<<endl;CRDGEODETICtoCRDCARTESEAN(pcg, pcc, dSemiMajorAxis,dFlatning);cout<<endl;cout<<"转换后的空间直角坐标:";cout<<"x="<<pcc->x<<" "<<"y="<<pcc->y<<" "<<"z="<<pcc->z<<" "<<endl;CRDCARTESIANtoCRDGEODETIC(pcc,pcg,dSemiMajorAxis,dFlatning);cout<<endl;cout<<"由转换后的空间直角坐标x,y,z还原成大地坐标:"<<endl;cout<<"H(高度)="<<pcg->height<<" "<<"L(经度)="<<pcg->longitude<<" "<<"B(纬度)="<<pcg->latitude<<endl;cout<<endl;}。