8高斯投影正反算 C#代码

合集下载

高斯投影正反算c代码

高斯投影正反算c代码

高斯投影正反算程序设计一.程序设计流程本程序(de)设计思路如下:(1),程序采用VS08版本作为开发平台,并采用C语言作为开发语言,设计为WindowsForm窗体程序形式.(2),程序主要(de)算法来自于教材.但是本程序为了更加实用,添加了更多(de)解算基准,包括:WGS-84,国际椭球1975,克氏椭球,和2000国家大地坐标系.(3),程序为了更方便(de)读取数据和输出数据,故需要自己定义了固定(de)数据输入格式和数据输出格式或形式,请老师注意查看.二.代码using System;using Systusing SystemponentModel;using System.Data;using System.Drawing;using System.Text;namespace Gauss{public partial class Form1 : Form{//大地坐标//Geodetic Coordinatepublic struct CRDGEODETIC{public double dLongitude;public double dLatitude;public double dHeight;}//笛卡尔坐标//Cartesian Coordinatepublic struct CRDCARTESIAN{public double x;public double y;public double z;}public Form1(){InitializeComponent();}private void button1_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;try{}catch{MessageBox.Show("Gauss Inverse: Choose datum error");return;}if (ttpareTo("克氏椭球")==0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;ee =}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}const double pai = 3.1415926;double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b);CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong;//求纬度string[] temp;double[] tempradius = new double[3];for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLatitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;//求经度for (int i = 0; i < 3; i++){tempradius[i] = Convert.ToDouble(temp[i]);}pcrdGeo.dLongitude = tempradius[0] / 180.0 pai + tempradius[1] / 180.0 / 60.0 pai + tempradius[2] / 180 / 60.0 / 60 pai;int deglon = Convert.ToInt32(pcrdGeo.dLongitude180 / pai);//求中央经度int num; //带号midlong = 0; //默认值,需要制定分带try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){num = Convert.ToInt32(deglon / 6 + 1);midlong = (6 num - 3) / 180.0 pai;}if (ttpareTo("6度带") == 0){num = Convert.ToInt32((deglon + 1.5) / 3);midlong = num 3 pai / 180;}double lp=pcrdGeo.dLongitude - midlong;double N = c / Math.Sqrt(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude));double M = c / Math.Pow(1 + epp eppMath.Cos(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude), 1.5);double ita = epp Math.Cos(pcrdGeo.dLatitude);double t = Math.Tan(pcrdGeo.dLatitude);double Nscnb = N Math.Sin(pcrdGeo.dLatitude) Math.Cos(pcrdGeo.dLatitude);double Ncosb = N Math.Cos(pcrdGeo.dLatitude);double cosb = Math.Cos(pcrdGeo.dLatitude);double X;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double B = pcrdGeo.dLatitude;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);X = a0 B - a2 / 2.0 s2b + a4 s4b / 4.0 - a6 / 6.0 s6b;pcrdCar.x = Nscnb lp lp / 2.0 + Nscnb cosb cosb Math.Pow(lp, 4) (5 - t t + 9 ita ita + 4 Math.Pow(ita, 4)) / 24.0 + Nscnb Math.Pow(cosb, 4) Math.Pow(lp, 6) (61 - 58 t t + Math.Pow(t, 4)) / 720.0 + X;pcrdCar.y = Ncosb Math.Pow(lp, 1) + Ncosb cosb cosb (1 - t t + ita ita) / 6.0 Math.Pow(lp, 3) + NcosbMath.Pow(lp, 5) Math.Pow(cosb, 4) (5 - 18 t t+ Math.Pow(t, 4) + 14 ita ita - 58 ita ita t t) / 120.0 ;if (pcrdCar.y < 0)pcrdCar.y += 500000;richTextBox1.Text = "Results:\nX:\t" +Convert.ToString(pcrdCar.x) +"\nY:\t"+ Convert.ToString(pcrdCar.y);}private void button2_Click(object sender, EventArgs e) {double ee = 0;double a = 0;string tt;int num; //带号string ytext; //利用y值求带号和中央经线try{}catch{MessageBox.Show("Gauss Inverse: Choose datumerror");return;}if (ttpareTo("克氏椭球") == 0){a = 6378245.00;}if (ttpareTo("WGS-84") == 0){a = 6378.00;}if (ttpareTo("1975国际椭球") == 0){a = 6378140.00;}if (ttpareTo("2000国家大地坐标系") == 0){a = 6378.0;}double b = Math.Sqrt(a a (1 - ee ee));double c = a a / b;double epp = Math.Sqrt((a a - b b) / b / b); CRDGEODETIC pcrdGeo;CRDCARTESIAN pcrdCar;double midlong = 0;//求X,Y和带号pcrdCar.x = Convert.ToDouble(textBox4.Text); ytext = textBox5.Text;string temp = ytext.Substring(0, 2);num = Convert.ToInt32(temp);ytext = ytext.Remove(0, 2);pcrdCar.y = Convert.ToDouble(ytext) - 500000; try{}catch{MessageBox.Show("Choose 3/6 error");return;}if (ttpareTo("3度带") == 0){midlong = num 3 pai / 180;}if (ttpareTo("6度带") == 0){midlong = (6 num - 3) pai / 180;}b = Math.Sqrt(a a (1 - ee ee));c = a a / b;epp = Math.Sqrt(a a - b b) / b;double m0, m2, m4, m6, m8;double a0, a2, a4, a6, a8;m0 = a (1 - ee ee);m2 = 3.0 / 2.0 m0 ee ee;m4 = 5.0 / 4.0 ee ee m2;m6 = 7.0 / 6.0 ee ee m4;m8 = 9.0 / 8.0 ee ee m6;a0 = m0 + m2 / 2.0 + 3.0 / 8.0 m4 + 5.0 / 16.0 m6 + 35.0 / 128.0 m8;a2 = m2 / 2 + m4 / 2 + 15.0 / 32.0 m6 + 7.0 / 16.0 m8;a4 = m4 / 8.0 + 3.0 / 16.0 m6 + 7.0 / 32.0 m8;a6 = m6 / 32.0 + m8 / 16.0;a8 = m8 / 128.0;double Bf, B;Bf = pcrdCar.x / a0;B = 0.0;while (Math.Abs(Bf - B) > 1E-10){B = Bf;double sb = Math.Sin(B);double cb = Math.Cos(B);double s2b = sb cb 2;double s4b = s2b (1 - 2 sb sb) 2;double s6b = s2b Math.Sqrt(1 - s4b s4b) + s4b Math.Sqrt(1 - s2b s2b);Bf = (pcrdCar.x - (-a2 / 2.0 s2b + a4 / 4.0 s4b - a6 / 6.0 s6b)) / a0;}double itaf, tf, Vf, Nf;itaf = epp Math.Cos(Bf);tf = Math.Tan(Bf);Vf = Math.Sqrt(1 + epp epp Math.Cos(Bf)Math.Cos(Bf));Nf = c / Vf;double ynf = pcrdCar.y / Nf;pcrdGeo.dLatitude = Bf - 1.0 / 2.0 Vf Vf tf (ynf ynf - 1.0 / 12.0 Math.Pow(ynf, 4) (5 + 3 tf tf + itaf itaf - 9 Math.Pow(itaf tf, 2)) +1.0 / 360.0 (61 + 90 tf tf + 45 Math.Pow(tf, 4)) Math.Pow(ynf, 6));pcrdGeo.dLongitude = (ynf / Math.Cos(Bf) - (1 + 2 tf tf + itaf itaf) Math.Pow(ynf, 3) / 6.0 / Math.Cos(Bf) +(5 + 28 tf tf + 24 Math.Pow(tf, 4) + 6 itaf itaf + 8 Math.Pow(itaf tf, 2)) Math.Pow(ynf, 5) / 120.0 /Math.Cos(Bf));pcrdGeo.dLongitude = pcrdGeo.dLongitude + midlong;//pcrdGeo.dLatitude = pcrdGeo.dLatitude;richTextBox2.Text = "Results:\nLatitude: " + Convert.ToString(pcrdGeo.dLatitude) + "\nLongtitude: " +Convert.ToString(pcrdGeo.dLongitude);}private void label13_Click(object sender, EventArgs e) {}}}三.程序运行结果分析通过选取书上(de)具体实例进行测试,本程序(de)精度大体满足要求,一般正算(de)精度在0.01米和0.001米之间,反算(de)精度在0.0001秒左右.以下是程序运行(de)截图.。

高斯投影正反算

高斯投影正反算

高斯投影正反算学院:资源与环境工程工程学院专业:测绘工程 学号:X51414012:超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面,并与某一条子午线相切,椭圆柱的中心轴通过椭球体的中心,然后用一定投影方法,将中央子午线两侧各一定经差围的地区投影到椭圆柱面上,再将此柱面展开即成为投影面。

高斯投影由于是正形投影,故保证了投影的角度不变性,图形的相似性以及在某点各方向上长度比的同一性。

由于采用了同样法则的分带投影,这即限制了长度变形,又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算,并且带与带间的互相换算也能用相同的公式和方法进行。

高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。

二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件 1)中央子午线投影后为直线 2)中央子午线投影后长度不变 3)投影具有正形性质,即正形投影条件2.高斯正算公式推导1)由第一个条件可知,由于地球椭球体是一个旋转椭球体,所以高斯投影必然有这样一个性质,即中央子午线东西两侧的投影必然对称于中央子午线。

2)由于高斯投影是换带投影,在每带经差l是不大的,lρ是一个微小量,所以可以将X=X (l,q ),Y=Y (l ,q )展开为经差为l 的幂级数,它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数,他们都是纬度B 的函数。

3)由第三个条件:∂y ∂l =∂x ∂q 和∂x ∂l =-∂y∂q ,将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩经过计算可以得出232244524632235242225sin cos sin cos (594)224 sin cos (6158)720cos cos (1)6cos (5181458)120N N x X B B l B B t l NB B t t l Ny N B l B t l NB t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路:级数展开,应用高斯投影三个条件,待定系数法求解。

(整理)高斯投影坐标正反算公式及适合电算的高斯投影公式

(整理)高斯投影坐标正反算公式及适合电算的高斯投影公式

§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的;②如果是正形投影,除了满足正形投影的条件外(C-R 偏微分方程),还有它本身的特殊条件。

8.3.1高斯投影坐标正算公式: B,l ⇒ x,y高斯投影必须满足以下三个条件:①中央子午线投影后为直线;②中央子午线投影后长度不变;③投影具有正形性质,即正形投影条件。

由第一条件知中央子午线东西两侧的投影必然对称于中央子午线,即(8-10)式中,x 为l 的偶函数,y 为l 的奇函数;0330'≤l ,即20/1/≈''''ρl ,如展开为l 的级数,收敛。

+++=++++=553316644220l m l m l m y l m l m l m m x (8-33)式中 ,,10m m 是待定系数,它们都是纬度B 的函数。

由第三个条件知:qy l x l y q x ∂∂-=∂∂∂∂=∂∂, (8-33)式分别对l 和q 求偏导数并代入上式----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m (8-34) 上两式两边相等,其必要充分条件是同次幂l 前的系数应相等,即dq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)(8-35)是一种递推公式,只要确定了0m 就可依次确定其余各系数。

由第二条件知:位于中央子午线上的点,投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ,即(8-33)式第一式中,当0=l时有:0m X x == (8-36) 顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2==== 得:B V cB N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-37,38)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-= (8-39)依次求得6543,,,m m m m 并代入(8-33)式,得到高斯投影正算公式6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ (8-42) 8.3.2高斯投影坐标反算公式x,y⇒B,l投影方程:),(),(21y x l y x B ϕϕ== (8-43)满足以下三个条件:①x 坐标轴投影后为中央子午线是投影的对称轴;② x 坐标轴投影后长度不变;③投影具有正形性质,即正形投影条件。

高斯平面直角坐标与大地坐标的相互转换——高斯投影的正算与反算.

高斯平面直角坐标与大地坐标的相互转换——高斯投影的正算与反算.

昆明冶金高等专科学校测绘学院 (4)计算公式
3 2 2 2 4 ( 5 3 t 9 t ) y f f f f 2M f N f 2 4M f N 3 f tf 2 4 6 (6 1 9 0t f 4 5t f ) y 7 2 0M f N 5 f 1 1 2 2 3 l y (1 2t f f ) y 3 N f co s B f 6 N f co s B f 1 2 5 (5 2 8t 2 t4 2 2 f 24 f 6 f 8 f t f )y 5 1 2 0N f co s B f B Bf tf y2 tf
式中:

2 e 2 cos2 B
t 2 tan2 B l (L L0) X为B对应子午线弧长 N为卯酉圈曲率半径 20626 5
昆明冶金高等专科学校测绘学院
2
高斯投影坐标反算公式
(1)高斯投影反算:
已知某点 x, y ,求该点 L, B ,即 x, y ( L, B) 的坐标变换。 (2)投影变换必须满足的条件
昆明冶金高等专科学校测绘学院
二、高斯投影坐标正反算得实用公式及算例
1 高斯投影坐标正算公式 (1)高斯投影正算: 已知某点的 L, B ,求该点的 x, y ,即 (2)投影变换必须满足的条件: 中央子午线投影后为直线; 中央子午线投影后长度不变; 投影具有正形性质,即正形投影条件。 (3)投影过程 在椭球面上有对称于中央子午线的两点 P1 和 P2 ,它们的大地坐标 分别为 ( L1 , B1 )或(l1 , B1)及 (L2 , B2)或(l2 , B2 ) 式中 l 为椭球面上点的经 度与中央子午线 ( L0 ) 的经度差:l L L0 ,点在中央子午线之东, l 为正,在西则为负,则投影后的平面坐标一定为P1 ( x1 , y1 ) 和 P2 ( x 2 , y 2 ) 。

高斯投影正反算原理

高斯投影正反算原理

高斯投影正反算原理高斯投影是一种常用于地图制图的投影方式,也被广泛应用于其他领域的空间数据处理。

高斯投影正反算是对于已知的地球坐标系上的位置(经纬度),通过计算得到该点的平面坐标(东、北坐标),或者对于已知的平面坐标(东、北坐标),通过计算得到该点的地球坐标系上的位置(经纬度)的过程。

本文将详细介绍高斯投影正反算的原理。

一、高斯投影简介高斯投影是一种圆锥投影,其投影面在地球表面的某个经线上,也就是说,投影面是以该经线为轴的圆锥面。

经过对圆锥体的调整后,使其切于地球椭球面,在该经线上进行投影,同时保持沿该经线方向的比例尺一致,从而达到地图上各点在包括该经线的垂直面上映射的目的。

这种投影方式在某一特定区域内得到高精度的结果,因此广泛应用于地图制图。

二、高斯投影数学模型对于高斯投影正反算,需要先建立高斯投影坐标系与地球坐标系的转换模型。

1.高斯投影坐标系的建立高斯投影坐标系的建立需要确定圆锥面的基本参数,首先需要确定其所处的中央子午线,再确定该子午线上的经度为零点,并利用该经线上某一点的经度和该点的高度来确定该点所在的圆锥体。

圆锥体的底面包括所有与地球椭球面相切的圆面,通过对这些圆面进行调整,使得圆锥体转动后能够在中央子午线上进行投影。

在此基础上,可建立高斯投影坐标系,其中投影面为圆锥面,且中央子午线与投影面的交点称为该投影坐标系的中心,投影面的上端点和下端点分别对应正北方向和正南方向。

2.地球坐标系的建立地球坐标系是以地球椭球体为基础建立的,其坐标系原点确定为地球椭球体上的一个特定点。

在已知该点经纬度和高度的前提下,可确定以该点为中心的地球椭球体,并可根据它与地球坐标系之间的转换关系得到平面坐标系。

3.高斯投影坐标系与地球坐标系之间的转换关系由于高斯投影坐标系与地球坐标系存在不同的坐标体系和基准面,因此需要通过数学关系式来建立它们之间的转换关系。

(1)高斯投影坐标系转地球坐标系:已知高斯投影坐标系中任意一点的东北坐标(N,E),以及所属的中央子午线经度λ0、椭球参数a和e,则可通过以下公式求出该点的地球坐标系经纬度(φ,λ)和高度H:A0为以地球椭球体中心为原点,高斯投影坐标系中心投影坐标为(0,0)的点到椭球面的距离。

高斯投影正反算

高斯投影正反算

高斯投影正反算LT2340123423401234...........q m m y m y m y m y l n n y n y n y n y '''''=+++++'''''=+++++3.引入高斯投影条件之一:正形条件4.由于可得到00n '=,带入上式可得到 2402435135...........q m m y m y l n y n y n y '''=+++'''=+++5.引入高斯投影条件之三:中央子午线投影后长度不变001122223322442255sec sec 122sec (12)6sec (56)24sec (5286)120f f f f f f f f f f f f f f f f f f f m q B dm n dx N t B dn m dx N B n t N t B m t N B n t N ηηη'=⎧⎪'⎪'==⎪⎪⎪''=-=-⎪⎪⎪⎨'=-++⎪⎪⎪'=++⎪⎪⎪⎪'=++⎪⎩四、高斯投影的特点1.当l 等于常数时,随着B 的增加x 的值增大,y 的值减小,无论B 值为正或为负,y 值不变。

这就是说,椭球面上除中央子午线外,其它子午线投影后,均向中央子午线弯曲,并向两极收敛,同时还对称于中央子午线和赤道。

2.当B 等于常数时,随着l 的增加,x 值和y 值都增大。

所以在椭球面上对称于赤道的纬圈,投影后仍成为对称的曲线,同时与子午线的投影曲线互相垂直凹向两极。

3.据中央子午线越远的子午线,投影后弯曲越厉害,长度变形越大。

五、MATLAB编程实现坐标正反算1.编写main函数function maindisp('欢迎使用高斯投影正反算及相邻带的坐标换算程序');disp('1:高斯正算 2:高斯反算 3:换带计算');K=0;while (K<1||K>3)K=input('请根据上列选择计算类型 K=');switch Kcase 1GSZS;case 2GSFS;case 3HDJS;otherwisedisp('K 值无效(1-3)');enddisp('程序作者:亚里士多墩');disp('指导老师:亚里士多德');end2.编写高斯正算GSZS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3-0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效(1-3)');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24; xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4);r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End3.编写高斯反算公式GSFS函数function GSZS%GSZS 是将大地坐标换算为高斯坐标的子函数%此函数要调用 DHH 和 HHD 两个子函数%此函数包含子午线收敛角的计算disp('你选择的是高斯正算');B=input('输入大地坐标 B=');L=input('输入大地坐标 L=');L0=input('输入所用中央子午线 L0=');B=DHH(B);L=DHH(L);L0=DHH(L0);disp('1:克拉索夫斯基椭球 2:1975 年国际椭球 3:WGS-84 椭球');T=0;while (T<1||T>3)T=input('请根据上列选择椭球模型 T=');switch Tcase 1a=6378245.0000000000;b=6356863.0187730473;X=111134.861*(B*180/pi)-16036.480*sin(2*B)+16.828*sin(4*B)-0.022*sin (6*B);case 2a=6378140.0000000000;b=6356755.2881575287;X=111133.005*(B*180/pi)-16038.528*sin(2*B)+16.833*sin(4*B)-0.022*sin (6*B);case 3a=6378137.0000000000;f=1/298.257223563;b=a*(1-f);X=6367449.1458*B-32009.8185*cos(B)*sin(B)-133.9975*cos(B)*(sin(B))^3 -0.6975*cos(B)*(sin(B))^5;otherwisedisp('T 值无效(1-3)');endende=(sqrt(a^2-b^2))/a;e1=(sqrt(a^2-b^2))/b;V=sqrt(1+(e1^2)*(cos(B))^2);c=(a^2)/b;M=c/(V^3);N=c/V;t=tan(B);n=sqrt((e1^2)*(cos(B))^2);l=L-L0;xp1=X;xp2=(N*sin(B)*cos(B)*l^2)/2;xp3=(N*sin(B)*((cos(B))^3)*(5-t^2+9*n^2+4*n^4)*l^4)/24;xp4=(N*sin(B)*((cos(B))^5)*(61-58*t^2+t^4)*l^6)/720;x=xp1+xp2+xp3+xp4;yp1=N*cos(B)*l;yp2=N*(cos(B))^3*(1-t^2+n^2)*l^3/6;yp3=N*(cos(B))^5*(5-18*t^2+t^4)*l^5/120;y=yp1+yp2+yp3;r1=l*sin(B);r2=(1/3)*sin(B)*(cos(B))^2*(l^3)*(1+3*n^2+2*n^4); r3=(1/15)*sin(B)*(cos(B))^2*(l^5)*(2-t^2);r=r1+r2+r3;format long gxyR=HHD(r)End六、代码测试1.L=111°47'24〞.8974,B=31°04'41〞.6832,L0=111°2.L=111.47532575,B=31.23484275,L0=1113.L=114.20,B=30.30,L0=1114.L=118.54152206,B=32.24576522,L0=117。

高斯投影坐标反算c语言代码

高斯投影坐标反算c语言代码

高斯投影坐标反算c语言代码#include<stdio.h>#include<math.h>#include<conio.h>main(){printf("#####################################################\n");printf("# 角度输入说明:如26°12′45.2″输入为26,12,45.2 #\n");printf("#####################################################\n");double x,y;int j,L0;printf("请输入高斯投影坐标(自然坐标),中间用逗号隔开:\n");scanf("%lf,%lf",&x,&y);//自然坐标输入printf("请输入中央子午线L0:\n");scanf("%d,%d,%lf",&L0); //中央子午线输入printf("请选择参考椭球:1.北京1954参考椭球。

\n 2.西安1980参考椭球。

\n");printf("选择的参考椭球为:");scanf("%d",&j);//选择椭球参数if(j==1){long double Bf0=0.157046064172*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005051773759*sin(Bf0)-0.000029837302*pow(sin(Bf0),3)+0.00000023818 9*pow(sin(Bf0),5));long double t=tan(Bf);long double m=0.00673852541468*pow(cos(Bf),2);long double V=1+m;long double N=6378245.000/sqrt(1-0.00669342162297*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d3)*60);long double m3=((L-d3)*60-f3)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //北京1954参考椭球}if(j==2){long double Bf0=0.157048687473*pow(10,-6)*x;long double Bf=Bf0+cos(Bf0)*(0.005052505593*sin(Bf0)-0.000029847335*pow(sin(Bf0),3)+0.0000002416* pow(sin(Bf0),5)+0.0000000022*pow(sin(Bf0),7));long double t=tan(Bf);long double m=0.006739501819473*pow(cos(Bf),2);long double V=1+m;long double N=6378140.000/sqrt(1-0.006694384999588*pow(sin(Bf),2));long double B1=Bf-1.0/2*V*t*pow(y/N,2)+1.0/24*(5+3*pow(t,2)+m-9*m*pow(t,2))*V*t*pow(y/N,4)-1.0/72 0*(61+90*pow(t,2)+45*pow(t,4))*V*t*pow(y/N,6);long double l1=(1/cos(Bf))*(y/N)-1.0/6*(1+2*pow(t,2)+m)*(1/cos(Bf))*pow(y/N,3)+1.0/120*(5+28*pow(t,2)+24*pow(t,4)+6*m+8*m*pow(t,2))*(1/cos(Bf))*pow(y/N,5);long double B=B1*57.29577951;long double l=l1*57.29577951;long double L=L0+l;int d2=int(B);int f2=int((B-d2)*60);long double m2=((B-d2)*60-f2)*60;printf("B=%d %d %.12lf\n",d2,f2,m2);int d3=int(L);int f3=int((L-d2)*60);long double m3=((L-d2)*60-f2)*60;printf("L=%d %d %.12lf\n",d3,f3,m3); //西安1980参考椭球}getch();}。

测绘程序设计(VS2008)实验报告--高斯投影正反算

测绘程序设计(VS2008)实验报告--高斯投影正反算

《测绘程序设计()》上机实验报告(Visual C++.Net)班级:学号:姓名:序号:二零一一年五月实验7 常用测量程序设计1.实验目的:1.1 巩固类的创建与使用;1.2 掌握数组参数的传递;1.3 掌握常用测绘程序设计的技巧。

1.42.实验内容:编写高斯投影正反算程序。

3.设计思路:这次的实验目的是实现高斯正反算。

需要考虑投影方式即分带的方式,又要考虑椭球参数的类型,所以我添加了两个函数来完成此功能。

分别是int SetProjectType(int m)和void SetParameter(int m,double &a,double &b)。

4.界面设计:界面设计很简单,具体见运行结果。

5.主要代码:文件名:GaussProjectDlg.cpp代码:const double PI=4*atan(1.0);//获得分带方式返回中央子午线经度int CGaussProjectDlg::SetProjectType(int m){UpdateData(TRUE);int n; //记录分带带号double L; //经度L=iDegreeL+iMinL/60+dSecondL/3600;if(m==1) //6度带{n=int(L/6)+1;L0=6*n-3;}else if(m==2) //3度带{n=int((L+1.5)/3);L0=3*n;}else if(m==3) //自主分带L0=L0;return L0;}//获取椭球参数void CGaussProjectDlg::SetParameter(int m,double &a,double &b) {if(m==1) //克拉索夫斯基椭球{a=6378245.0;b=6356863.0187730473;//e=sqrt(0.006693421622966);}else if(m==2) //1975国际协议椭球{a=6378140.0;b=6356755.2881575287;//e=sqrt(0.006694384999588);}else if(m==3) //WGS-84椭球{a=6378137.0;b=6356752.3142;//e=sqrt(0.0066943799013);}}void CGaussProjectDlg::OnBnClickedButtonpositivecal(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double N;double t;double Eta;double X;double A0,A2,A4,A6,A8;double RadB;double Rou;Rou=180*3600/PI;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double l;L0=SetProjectType(iProjectType);double L;L=iDegreeL+double(iMinL)/60+dSecondL/3600;l=(L-L0)*3600;RadB=(iDegreeB+double(iMinB)/60+dSecondB/3600)*PI/180;N=a/sqrt(1-e1*e1*sin(RadB)*sin(RadB));t=tan(RadB);Eta=e2*cos(RadB);A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);A2=-1.0/2*(3.0/4*e1*e1+60.0/64*pow(e1,4)+525.0/512*pow(e1,6)+17640.0/16384*pow(e1,8));A4=1.0/4*(15.0/64*pow(e1,4)+210.0/512*pow(e1,6)+8820.0/16384*pow(e1,8));A6=-1.0/6*(35.0/512*pow(e1,6)+2520.0/16384*pow(e1,8));A8=1.0/8*(315.0/16384*pow(e1,8));X=a*(1-e1*e1)*(A0*RadB+A2*sin(2*RadB)+A4*sin(4*RadB)+A6*sin(6*RadB)+A8*sin(8*RadB));x=X+N/(2*Rou*Rou)*sin(RadB)*cos(RadB)*l*l+N/(24*pow(Rou,4))*sin(RadB)*pow(cos(RadB),3)*(5-t*t+9*Eta*Eta+4*pow(Eta,4))*pow(l,4)+ N/(720*pow(Rou,6))*sin(RadB)*pow(cos(RadB),5)*(61-58*t*t+pow(t,4))*pow(l,6);y=N/Rou*cos(RadB)*l+N/(6*pow(Rou,3))*pow(cos(RadB),3)*(1-t*t+Eta*Eta)*pow(l,3)+N/(120*pow(Rou,5))*pow(cos(RadB),5)*(5-18*t*t+pow(t,4)+14*Eta*Eta-58*Eta*Eta*t*t)*pow( l,5);UpdateData(FALSE);}void CGaussProjectDlg::OnBnClickedButtonantical(){// TODO: 在此添加控件通知处理程序代码UpdateData(TRUE);double t_f;double Eta_f;double B_f;double N_f;double M_f;double X=x;double B0;double K0,K2,K4,K6;double a,b,e1,e2; //椭球参数SetParameter(iParameterType,a,b);e1=sqrt(a*a-b*b)/a;e2=sqrt(a*a-b*b)/b;double A0;A0=1+3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8);B0=X/(a*(1-e1*e1)*A0);K0=1.0/2*(3.0/4*e1*e1+45.0/64*pow(e1,4)+350.0/512*pow(e1,6)+11025.0/16384*pow(e1,8));K2=-1.0/3*(63.0/64*pow(e1,4)+1108.0/512*pow(e1,6)+58239.0/16384*pow(e1,8));K4=1.0/3*(604.0/512*pow(e1,6)+68484.0/16384*pow(e1,8));K6=-1.0/3*(26328.0/16384*pow(e1,8));B_f=B0+sin(2*B0)*(K0+sin(B0)*sin(B0)*(K4+K6*sin(B0)*sin(B0)));t_f=tan(B_f);Eta_f=e2*cos(B_f);N_f=a/sqrt(1-e1*e1*sin(B_f)*sin(B_f));M_f=N_f/(1+e2*e2*cos(B_f)*cos(B_f));double B;B=B_f-t_f/(2*M_f*N_f)*y*y+t_f/(24*M_f*pow(N_f,3))*(5+3*t_f*t_f+Eta_f*Eta_f-9*Eta_f*Eta_f*t_f*t_f)*pow(y,4)- t_f/(720*M_f*pow(N_f,5))*(61+90*t_f*t_f+45*pow(t_f,4))*pow(y,6);double l;l=1.0/(N_f*cos(B_f))*y-1.0/(6*pow(N_f,3)*cos(B_f))*(1+2*t_f*t_f+Eta_f*Eta_f)*pow(y,3)+1.0/(120*pow(N_f,5)*cos(B_f))*(5+28*t_f*t_f+24*pow(t_f,4)+6*Eta_f*Eta_f+8*Eta_f*Eta_f* t_f*t_f)*pow(y,5);//将B转化为度分秒的形式double dDegB;dDegB=B*180/PI;iDegreeB=int(dDegB);iMinB=int((dDegB-iDegreeB)*60);dSecondB=((dDegB-iDegreeB)*60-iMinB)*60;double dDegL;dDegL=l*180/PI+L0;iDegreeL=int(dDegL);iMinL=int((dDegL-iDegreeL)*60);dSecondL=((dDegL-iDegreeL)*60-iMinL)*60;UpdateData(FALSE);}6.运行结果:实验的运行结果如下图所示:正算:反算:7.实验总结这次实验是实现高斯投影的正反算,方法很多,实现并不复杂,但是计算公式复杂,变量繁多,稍有不慎,就会造成计算错误。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
private void button1_Click(object sender, EventArgs e) {
double ee = 0; double a = 0; string tt;
try {
tt = boBox1.Items[boBox1.SelectedIndex].ToString(); } catch {
X = a0 * B - a2 / 2.0 * s2b + a4 * s4b / 4.0 - a6 / 6.0 * s6b;
pcrdCar.x = Nscnb * lp * lp / 2.0 + Nscnb * cosb * cosb * Math.Pow(lp, 4) * (5 - t * t + 9 * ita * ita + 4 * Math.Pow(ita, 4)) / 24.0
MessageBox.Show("Gauss Inverse: Choose datum error!"); return; }
if (pareTo("克氏椭球")==0)
{ a = 6378245.00; ee = Math.Sqrt(0.006693421622);
} if (pareTo("WGS-84") == 0) {
double ita = epp * Math.Cos(pcrdGeo.dLatitude); double t = Math.Tan(pcrdGeo.dLatitude); double Nscnb = N * Math.Sin(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude); double Ncosb = N * Math.Cos(pcrdGeo.dLatitude); double cosb = Math.Cos(pcrdGeo.dLatitude);
二.代码
using System; using System.Collections.Generic; using ponentModel; using System.Data; using System.Drawing; using System.Text; using System.Windows.Forms;
//求中央经度 int num; midlong = 0;
//带号 //默认值,需要制定分带
try {
tt = boBox3.Items[boBox3.SelectedIndex].ToString(); } catch {
MessageBox.Show("Choose 3/6 error!"); return; } if (pareTo("3度带") == 0) {
a = 6378137.00; ee = Math.Sqrt(0.00669437999013); } if (pareTo("1975国际椭球") == 0) { a = 6378140.00; ee = Math.Sqrt(0.006694384999588); } if (pareTo("2000国家大地坐标系") == 0) { a = 6378137.0; ee =Math.Sqrt(0.0066943802290); }
//求纬度 string[] temp; temp = textBox1.Text.Split(' '); double[] tempradius = new double[3];
for (int i = 0; i < 3; i++) {
tempradius[i] = Convert.ToDouble(temp[i]); }
num = Convert.ToInt32(deglon / 6 + 1); midlong = (6 * num - 3) / 180.0 * pai; } if (pareTo("6度带") == 0) { num = Convert.ToInt32((deglon + 1.5) / 3); midlong = num * 3 * pai / 180; } double lp=pcrdGeo.dLongitude - midlong;
a8 = m8 / 128.0;
s2b);
double B = pcrdGeo.dLatitude; double sb = Math.Sin(B); double cb = Math.Cos(B); double s2b = sb * cb * 2; double s4b = s2b * (1 - 2 * sb * sb) * 2; double s6b = s2b * Math.Sqrt(1 - s4b * s4b) + s4b * Math.Sutton2_Click(object sender, EventArgs e) {
double ee = 0; double a = 0; string tt; int num;
//带号
string ytext;
//利用y值求带号和中央经线
try {
tt = boBox2.Items[boBox2.SelectedIndex].ToString(); } catch {
+ Nscnb * Math.Pow(cosb, 4) * Math.Pow(lp, 6) * (61 - 58 * t * t + Math.Pow(t, 4)) / 720.0 + X;
pcrdCar.y = Ncosb * Math.Pow(lp, 1) + Ncosb * cosb * cosb * (1 - t * t + ita * ita) / 6.0 * Math.Pow(lp, 3) + Ncosb * Math.Pow(lp, 5)* Math.Pow(cosb, 4) * (5 - 18 *t*t
double N = c / Math.Sqrt(1 + epp * epp * Math.Cos(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude));
double M = c / Math.Pow(1 + epp * epp * Math.Cos(pcrdGeo.dLatitude) * Math.Cos(pcrdGeo.dLatitude), 1.5);
高斯投影正反算程序设计
一.程序设计流程
本程序的设计思路如下: (1),程序采用 VS08 版本作为开发平台,并采用 C#语言作为开发语言,设计为 WindowsForm 窗体程序形式。 (2),程序主要的算法来自于教材。但是本程序为了更加实用,添加了更多的解算基准, 包括:WGS-84,国际椭球 1975,克氏椭球,和 2000 国家大地坐标系。 (3),程序为了更方便的读取数据和输出数据,故需要自己定义了固定的数据输入格式和 数据输出格式或形式,请老师注意查看。
MessageBox.Show("Gauss Inverse: Choose datum error!"); return; }
if (pareTo("克氏椭球") == 0) {
a = 6378245.00; ee = Math.Sqrt(0.006693421622); } if (pareTo("WGS-84") == 0) { a = 6378137.00; ee = Math.Sqrt(0.00669437999013); } if (pareTo("1975国际椭球") == 0) { a = 6378140.00; ee = Math.Sqrt(0.006694384999588); } if (pareTo("2000国家大地坐标系") == 0) { a = 6378137.0; ee =Math.Sqrt(0.0066943802290); }
const double pai = 3.1415926; double b = Math.Sqrt(a * a * (1 - ee * ee)); double c = a * a / b; double epp = Math.Sqrt((a * a - b * b) / b / b);
CRDGEODETIC pcrdGeo; CRDCARTESIAN pcrdCar; double midlong;
tempradius[i] = Convert.ToDouble(temp[i]); } pcrdGeo.dLongitude = tempradius[0] / 180.0 * pai + tempradius[1] / 180.0 / 60.0 * pai + tempradius[2] / 180 / 60.0 / 60 * pai; int deglon = Convert.ToInt32(pcrdGeo.dLongitude * 180 / pai);
}
//笛卡尔坐标 //Cartesian Coordinate public struct CRDCARTESIAN {
public double x; public double y; public double z; }
public Form1() {
InitializeComponent(); }
const double pai = 3.1415926535898; double b = Math.Sqrt(a * a * (1 - ee * ee)); double c = a * a / b; double epp = Math.Sqrt((a * a - b * b) / b / b); CRDGEODETIC pcrdGeo; CRDCARTESIAN pcrdCar; double midlong = 0;
相关文档
最新文档