高斯投影正算与反算的理论方法与实现代码

高斯投影正算与反算的理论方法与实现代码

高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系) ，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。

//6度带宽 54北京坐标系

//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)

void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) {

int ProjNo; int ZoneWide; ////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;

iPI = 0.[**************]3; ////3.[1**********]98/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 = (ProjNo-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));

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;

}

//高斯投影由经纬度(Unit:DD)正算平面坐标(含带号，Unit:Metres)

void GaussProjCal(double longitude, double latitude, double *X, double *Y) {

int ProjNo=0; int ZoneWide; ////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double a,f, e2,ee, NN, T,C,A, M, iPI;

iPI = 0.[**************]3; ////3.[1**********]98/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 = 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;

}

NN 卯酉圈曲率半径，测量学里面用N 表示

M 为子午线弧长，测量学里用大X 表示

fai 为底点纬度，由子午弧长反算公式得到，测量学里用Bf 表示

R 为底点所对的曲率半径，测量学里用Nf 表示

高斯投影正算与反算的理论方法与实现代码

高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系) ，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。

//6度带宽 54北京坐标系

//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)

void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) {

int ProjNo; int ZoneWide; ////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;

iPI = 0.[**************]3; ////3.[1**********]98/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 = (ProjNo-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));

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;

}

//高斯投影由经纬度(Unit:DD)正算平面坐标(含带号，Unit:Metres)

void GaussProjCal(double longitude, double latitude, double *X, double *Y) {

int ProjNo=0; int ZoneWide; ////带宽

double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval; double a,f, e2,ee, NN, T,C,A, M, iPI;

iPI = 0.[**************]3; ////3.[1**********]98/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 = 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;

}

NN 卯酉圈曲率半径，测量学里面用N 表示

M 为子午线弧长，测量学里用大X 表示

fai 为底点纬度，由子午弧长反算公式得到，测量学里用Bf 表示

R 为底点所对的曲率半径，测量学里用Nf 表示


相关文章

  • 图幅理论面积与图斑椭球面积计算公式及要求

    • 图幅理论面积与图斑椭球面积计算公式及要求 一. 图幅理论面积计算公式 (1) 式中: a-椭球长半轴(单位:米) ,α-椭球扁率,b-椭球短半轴(单位:米) . е2﹦(a2﹣b2)/a2. A ﹦1﹢(3/6)е2﹢(30/80)е4﹢( ...查看


  • 毕业设计论文_人脸识别系统

    • 摘要 人脸识别因其在安全验证系统.信用卡验证.医学.档案管理.视频会议.人机 交互.系统公安(罪犯识别等)等方面的巨大应用前景而越来越成为当前模式识别和人工智能领域的一个研究热点. 本文提出了基于24位彩色图像对人脸进行识别的方法,介绍的主 ...查看


  • 房地产测绘总结

    • 1. 测量学总结 第一章绪论 2. 测量学简介 3. 测量学是研究地球的形状.大小确定地球表面各种物体的形状.大小.空间位置的科学. 4. 地物和地貌 a) 地物--地面上天然或人工形成的物体, b) 包括湖泊.河流.海洋.房屋.道路与桥梁 ...查看


  • 测绘测量学基础知识

    • 测绘测量学基础知识 作者 : fffsky | 分类 : 测量资料 | 超过 1189人阅读 | 已有 0 人留下了评论 1. 什么是测绘 测绘是测量和地图制图的统称.测绘科学研究的对象主要是地球的形状.大小和地球表面的各种固定物体的几何形 ...查看


  • 2012年注册测绘师考试[测绘综合能力]真题及答案

    • 2012年注册测绘师考试<测绘综合能力>真题及答案 一.单项选择题(每题3 分,共10 题) 1.使用N 台(N>3)GPS 接收机进行同步观测所获得的GPS 边中,独立的GPS 边的数量是( ). A.N B.N-1 C ...查看


  • 控制测量学试题

    • 控制测量学试卷(A) 一.填空题(20分) 1.控制网一般分为:(1) 和(2) 两大类.(1)类解决控制点的坐标 ,(2)类解决控制点的高程. 2.导线网由于通视方向少,在布网时受障碍物限制较少,又由于边角同测,其横向位移小 于 ,纵向位 ...查看


  • 地籍测量学习题

    • 地籍测量学 1. 地籍测量的含义和内容是什么? 地籍测量是为获取和表达地籍信息所进行的测绘工作.其基本内容是测定土地及其附着物的权属.位置.数量.质量和利用状况等. 地籍测量的具体内容 (1) 地籍控制测量,测量地籍基本控制点和地籍图根控制 ...查看


  • 基于Matlab的人脸识别课程设计

    • 摘要 人脸识别因其在安全验证系统.信用卡验证.医学.档案管理.视频会议.人机交互.系统公安(罪犯识别等) 等方面的巨大应用前景而越来越成为当前模式识别和人工智能领域的一个研究热点. 本文提出了基于24位彩色图像对人脸进行识别的方法,介绍的主 ...查看


  • [大地测量学基础]试卷(2)

    • <大地测量学基础>课程试卷(2) 一.解释下列术语(每个2分,共10分) 1.大地水准面:假定海水面完全处于静止和平衡状态(没有风浪.潮汐及大气压变化的影响),把这个海水面伸延到大陆下面,形成一个封闭曲面,在这个面上都保持与重力 ...查看


热门内容