坐標轉換之高斯投影轉大地坐標（C#）

###以西安80坐標舉例###

1、大地原點位於陝西省涇陽縣永樂鎮，橢球參數采用IUG 1975年大會推薦的參數，基本參數為： 

 

2、解算過程包括：

高斯反算和子午線弧長反算

高斯反算需要用到底點緯度，由子午線弧長反算得到，中間涉及到迭代計算；

 

3、代碼

1）Coordin類

class Coordin
    {


        //西安80
        //高斯投影平面坐標值
        public double x_XA;
        public double y_XA;
        public double Hg_XA;
        
        //大地坐標值(角度°)
        public double B_XA_Angle;
        public double L_XA_Angle;
       
        //大地坐標值(弧長)
        public double B_XA;
        public double L_XA;
        public double H_XA;

        //空間直角坐標值
        public double X_XA;
        public double Y_XA;
        public double Z_XA;
    }

 

2）Ellips類

//常用橢球參數
public class Ellips
{
   //編號 為80時指西安80所用橢球，類推
　　public int num;

　　public double a;//長半軸a(m)
　　public double b;//長半軸b(m)
　　public double f;//扁率f
10 　　public double e1;//第一偏心率e1
　　public double e2;//第二偏心率e2
　　public double LevOrg85;//85水准原點(m)
　　public double L0;//中央子午線(度)

　　public class Ellips75_XA80
　　{
　　　　public const int num = 80;
　　　　//1975國際橢球_參心
　　　　public const double a = 6378140.000;//長半軸a(m)
　　　　public const double b = 6356755.288;//長半軸b(m)
　　　　public const double f = 1 / 298.257;//扁率f
　　　　public const double e1 = 0.006694385;//第一偏心率e1
　　　　public const double e2 = 0.006739502;//第二偏心率e2
　　　　public const double LevOrg85 = 72.2604;//85水准原點(m)

　　}
}

 

3）弧長反算

//子午線弧長反解
public double LengthInverse(double x_in, double a, double e1)
{
　　//基本常量m0,m2,m4,m6,m8
　　double[] mConst = new double[9];
　　//基本常量a0,a2,a4,a6,a8
　　double[] aConst = new double[9];

　　//計算常量m0,m2,m4,m6,m8
　　mConst[0] = a * (1 - e1 * e1);
　　mConst[2] = 3 * e1 * e1 * mConst[0] / 2.0;
　　mConst[4] = 5 * e1 * e1 * mConst[2] / 4.0;
　　mConst[6] = 7 * e1 * e1 * mConst[4] / 6.0;
　　mConst[8] = 9 * e1 * e1 * mConst[6] / 8.0;
　　//計算常量a0,a2,a4,a6,a8
　　aConst[0] = mConst[0] + mConst[2] / 2 + 3 * mConst[4] / 8 
　　　　　　　　　　  + 5 * mConst[6] / 16 + 35 * mConst[8] / 128;
　　aConst[2] = mConst[2] / 2 + mConst[4] / 2 + 15 * mConst[6] / 32  
　　　　　　　　　　  + 7 * mConst[8] / 16;
　　aConst[4] = mConst[4] / 8 + 3 * mConst[6] / 16 + 7 * mConst[8] / 32;
　　aConst[6] = mConst[6] / 32 + mConst[8] / 16;
　　aConst[8] = mConst[8] / 128;

　　//設置初始底點緯度B0
　　//進行迭代運算
　　double B0 = x_in / a;
　　bool Judge = true;
　　while (Judge)
　　{
　　　　double F = -aConst[2] * Math.Sin(2 * B0) / 2 
　　　　　　　　　　　　+ aConst[4] * Math.Sin(4 * B0) / 4 
　　　　　　　　　　　　- aConst[6] * Math.Sin(6 * B0) / 6 
　　　　　　　　　　　　+ aConst[8] * Math.Sin(8 * B0) / 8;
　　　　double B1 = (x_in - F) / a;
　　　　if (Math.Abs(B0 - B1) < 1.0E-6)
　　　　{
　　　　　　Judge = false;
　　　　}
　　　　else
　　　　{
　　　　　　B0 = B1;
　　　　}
　　}

　　return B0;
}

 

4）高斯反算

//高斯反算(xyHg->BLH)
public void GaussBacCalculation(Coordin coordin, Ellips ellips)
{
　　//初始化橢球參數
   int num = ellips.num;
   double a = ellips.a;
   double e1 = ellips.e1;
   double e2 = ellips.e2;
   double L0 = ellips.L0;
   double LevOrg85 = ellips.LevOrg85;
   //初始化為0
   double x = 0;
   double y = 0;
   double Hg = 0;
   //根據橢球判斷
   switch(num)
17　　 {
    　 case 54:
       //北京54高斯投影平面坐標
                        
       break;

       case 80:
       //西安80高斯投影平面坐標
       x = coordin.x_XA;
       y = coordin.y_XA;
       Hg = coordin.Hg_XA;
       break;
                    
30　　　}

                

34 　　//(弧長)解算底點緯度Bf
   double Bf = LengthInverse(x, a, e1);
   //中間變量Wf
   double Wf = Math.Sqrt(1 - Math.Pow(e1, 2) * Math.Sin(Bf));
   //子午圈曲率半徑Nf
   double Nf = a / Wf;
   //其余變量
   double Mf = a * (1 - Math.Pow(e1, 2)) / (Math.Pow(Wf, 3));
   double tf = Math.Tan(Bf);
   double pf = e2 * Math.Cos(Bf);
   //修正橫坐標y值
   y -= 500000;

   //(弧長)經度L
   double l = y / (Nf * Math.Cos(Bf))
              - Math.Pow(y, 3) * (
　　　　　　　　　　　　　　　　　　　　　　　　1 + 2 * Math.Pow(tf, 2) + Math.Pow(pf, 2)
　　　　　　　　　　　　　　　　　　　　　　) / (6 * Math.Pow(Nf, 3) * Math.Cos(Bf))
              + Math.Pow(y, 5) * (
　　　　　　　　　　　　　　　　　　　　　　　　5 + 28 * Math.Pow(tf, 2) + 24 * Math.Pow(tf, 4) 
　　　　　　　　　　　　　　　　　　　　　　　　+ 6 * Math.Pow(pf, 2) + 8 * Math.Pow(pf, 2) * Math.Pow(tf, 2)
　　　　　　　　　　　　　　　　　　　　　　) / (120 * Math.Pow(Nf, 5) * Math.Cos(Bf));
   double L = l + L0 * rPI;
   //(弧長)維度B
   double B = Bf - tf * Math.Pow(y, 2) / (2 * Mf * Nf) 
　　　　　　　　　　+ tf * Math.Pow(y, 4) * (
　　　　　　　　　　　　　　　　　　　　　　　　 　　5 + 3 * Math.Pow(tf, 2) + Math.Pow(pf, 2) 
　　　　　　　　　　　　　　　　　　　　　　　　 　　- 9 * Math.Pow(pf, 2) * Math.Pow(tf, 2)
　　　　　　　　　　　　　　　　　　　　　　　　 ) / (24 * Mf * Math.Pow(Nf, 3))
              - tf * Math.Pow(y, 6) * (
　　　　　　　　　　　　　　　　　　　　　　　　　　　61 + 90 * Math.Pow(tf, 2) + 45 * Math.Pow(tf, 4)
　　　　　　　　　　　　　　　　　　　　　　　　　) / (720 * Mf * Math.Pow(Nf, 5));
   //(不考慮高程異常)正常高轉大地高
   double H = Hg + LevOrg85;


   switch (num)
   {
　　　　case 54:
    　　//北京54大地坐標

　　　　break;

　　　　case 80:
　　　　//西安80大地坐標
68　　　　 　　coordin.L_XA = L;
          coordin.B_XA = B;
          coordin.H_XA = H;
       break;
　　}

                
}

 

 

 

參考博客及文檔

 

坐標轉換-大地轉高斯平面&平面坐標轉換

坐標高斯正反算

子午線弧長正反算

高斯投影正反算公式