問題

　　因為工作，需要把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};
    }

注：函式裡面的ａ，ｂ的值是ＣＧＣＳ２０００（和ＷＧＳ８４相同）下的。其他座標系下的大地座標向空間直角座標的轉換公式都是相同的，只是ａｂ的值需要改動。

利用七引數進行座標轉換

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

    /*七引數運算*/
    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);


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