1. 程式人生 > >利用七引數進行CGCS2000座標系到西安80座標系的轉換

利用七引數進行CGCS2000座標系到西安80座標系的轉換

問題

  因為工作,需要把CGCS2000座標系下的座標轉到西安80座標系下,中間由於用到了七引數,所以要進經過到空間直角座標系的轉換,然後再轉換到西安80大地座標下,最後再投影到西安80座標的某度帶。
  要求是輸入CGCS2000下的大地座標,最後輸出西安80下的平面座標。那麼這項工作可以分為以下幾個步驟:
  1.CGCS2000下的大地座標到CGCS2000下的空間直角座標的轉換;
  2.CGCS2000下的空間直角座標經過七引數轉換,得到西安80下的空間直角座標;
  3.西安80下的空間直角座標向西安80下的大地座標轉換;
  4.西安80下的大地座標進行高斯投影,得到平面座標。

大地座標到空間直角座標的轉換

  根據轉換公式:
圖片來源於網際網路

  private double[] dd2kjzj(double[] dd){
        double a=6378137,b=6356752.314;
        double ee=(a*a-b*b)/(a*a);
        double L=Math.toRadians(dd[0]),B=Math.toRadians(dd[1]),H=dd[2];
        double N=a/Math.sqrt(1-ee*Math.sin(B)*Math.sin(B));
        double X=(N+H)*Math.cos(B)*Math.cos
(L); double Y=(N+H)*Math.cos(B)*Math.sin(L); double Z=(N*(1-ee)+H)*Math.sin(B); return new double[]{X,Y,Z}; }

注:函式裡面的a,b的值是CGCS2000(和WGS84相同)下的。其他座標系下的大地座標向空間直角座標的轉換公式都是相同的,只是ab的值需要改動。

利用七引數進行座標轉換

  得到空間直角座標後就可以利用七引數進行座標轉換了。具體的七引數求解這裡不進行討論,有意向的網友可以自行百度。轉換公式如下圖所示:
  圖片來源於網際網路


  

    /*七引數運算*/
    private double[] qicanshu(double[] source){
        double tX,tY,tZ;

        tX=dx+source[0]*(1+k)+Oz*source[1]-Oy*source[2];
        tY=dy+source[1]*(1+k)-Oz*source[0]+Ox*source[2];
        tZ=dz+source[2]*(1+k)+Oy*source[0]-Ox*source[1];

        return new double[]{tX,tY,tZ};
    }

  方法中的source陣列是CGCS2000座標系下的空間直角座標
  dx,dy,dz是偏移量
  Ox,Oy,Oz是旋轉量
  k是縮放因子

空間直角座標到大地座標的轉換

  這個轉換就是前面所示的逆運算,如下圖
圖片來源於網際網路
  轉換比較複雜,需要進行迭代運算。

    private double[] kjzj2dd(double[] kjzj){
        double X=kjzj[0],Y=kjzj[1],Z=kjzj[2];
        double a=6378140;
        double f=1/298.257;
        double e2=2*f-f*f; //e^2;
        double L=Math.toDegrees(Math.atan(Y/X)+Math.PI);
        double B2=Math.atan(Z/Math.sqrt(X*X+Y*Y));
        double B1;
        double N;
        while (true){
            N=a/Math.sqrt(1-f*(2-f)*Math.sin(B2)*Math.sin(B2));
            B1=Math.atan((Z+N*f*(2-f)*Math.sin(B2))/Math.sqrt(X*X+Y*Y));
            if(Math.abs(B1-B2)<0.0000000001)
                break;
            B2=B1;
        }
        double H=Z/Math.sin(B2)-N*(1-e2);
        double B=Math.toDegrees(B2);

        return new double[]{L,B,H};
    }

高斯投影(正算)

  在得到西安80座標下的大地座標後,需要進行高斯投影,這樣才能夠得到平面座標。由於這個計算公式更加的複雜,公式什麼的就不再列了。
  

    //只適用於西安80下的座標投影,程式碼來源於網際網路,經測試還不錯
    private double[] dd2pm(double[] dd){
        double L=Math.toRadians(dd[0]),B=Math.toRadians(dd[1]);

        //輔助量
        double cosB = Math.cos(B);
        double sinB = Math.sin(B);
        double cosB_2 = cosB * cosB;
        double l = L - Math.toRadians(Lo);
        double ll = l * l;

        //計算係數
        double N = 6399596.652 - (21565.045 - (108.996 - 0.603 * cosB_2) * cosB_2) * cosB_2;
        double a0 = 32144.5189 - (135.3646 - (0.7034 - 0.0041 * cosB_2) * cosB_2) * cosB_2;
        double a4 = (0.25 + 0.00253 * cosB_2) * cosB_2 - 0.04167;
        double a6 = (0.166 * cosB_2 - 0.084) * cosB_2;
        double a3 = (0.3333333 + 0.001123 * cosB_2) * cosB_2 - 0.1666667;
        double a5 = 0.00878 - (0.1702 - 0.20382 * cosB_2) * cosB_2;

        //計算高斯平面座標值
        double x = 6367452.1328 * B - (a0 - (0.5 + (a4 + a6 * ll) * ll) * ll * N) * cosB * sinB;
        double y = (1 + (a3 + a5 * ll) * ll) * l * N * cosB + 500000;

        double[] xy = new double[2];
        xy[0] = x;
        xy[1] = y;
        return xy;
    }

  這段程式碼只適用於IAG75(即西安80)下的高斯投影,其他的座標系下的投影需要修改下引數。

總結

  至此,左右的轉換已經完成。用了大概一週吧(中間有很長一段時間糾結於誤差大的驚人,今天才知道甲方給的測試資料有問題),算是複習了下專業知識。
  下載原始碼

使用方法

(補充自2018年4月19日)
程式是用java寫的,封裝成了CoordsTransformer 類。如果要用於C#、Python語言需要自行修改。


// 首先獲取該類的例項
CoordsTransformer coordsTransformer=CoordsTransformer.getInstance();
// 然後獲取執行結果
double[] resutls=coordsTransformer.fromCgcs2Xian80(119.8774232,28.44428546,112.66);

如果資料量比較多的話,後面這句可以放在一個迴圈裡面。