問題

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

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

利用七參數進行坐標轉換

得到空間直角坐標后就可以利用七參數進行坐標轉換了。具體的七參數求解這里不進行討論，有意向的網友可以自行百度。轉換公式如下圖所示：
　　
　　經 貓醉 提醒，上圖中a1為縮放因子，a2, a3, a4分別為xyz的旋轉量（感謝）。

/*七參數運算*/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）下的高斯投影，其他的坐標系下的投影需要修改下參數。

總結

至此，左右的轉換已經完成。用了大概一周吧（中間有很長一段時間糾結于誤差大的驚人，今天才知道甲方給的測試數據有問題），算是復習了下專業知識。
　　下載源代碼

鏈接：https://pan.baidu.com/s/1sdglmFpIOhjTZvvp4zX_SQ 提取碼：cmyz

使用方法

（補充自2018年4月19日）
程序是用java寫的，封裝成了CoordsTransformer 類。如果要用于C#、Python語言需要自行修改。

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

如果數據量比較多的話，后面這句可以放在一個循環里面。

總結

以上是生活随笔為你收集整理的利用七参数进行CGCS2000坐标系到西安80坐标系的转换的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。