利用七参数进行CGCS2000坐标系到西安80坐标系的转换

Bertha 。 2022-08-19 15:11 49阅读 0赞

## 问题 ##

因为工作，需要把CGCS2000坐标系下的坐标转到西安80坐标系下，中间由于用到了七参数，所以要进经过到空间直角坐标系的转换，然后再转换到西安80大地坐标下，最后再投影到西安80坐标的某度带。  
　　要求是输入CGCS2000下的大地坐标，最后输出西安80下的平面坐标。那么这项工作可以分为以下几个步骤：  
　　1.CGCS2000下的大地坐标到CGCS2000下的空间直角坐标的转换；  
　　2.CGCS2000下的空间直角坐标经过七参数转换，得到西安80下的空间直角坐标；  
　　3.西安80下的空间直角坐标向西安80下的大地坐标转换；  
　　4.西安80下的大地坐标进行高斯投影，得到平面坐标。

## 大地坐标到空间直角坐标的转换 ##

根据转换公式：  
![图片来源于互联网][format_png]

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};
        }

注：函数里面的ａ，ｂ的值是ＣＧＣＳ２０００（和ＷＧＳ８４相同）下的。其他坐标系下的大地坐标向空间直角坐标的转换公式都是相同的，只是ａｂ的值需要改动。

## 利用七参数进行坐标转换 ##

得到空间直角坐标后就可以利用七参数进行坐标转换了。具体的七参数求解这里不进行讨论，有意向的网友可以自行百度。转换公式如下图所示：  
　　![图片来源于互联网][format_png 1]  
　　经 [猫醉][Link 1] 提醒，上图中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是缩放因子

## 空间直角坐标到大地坐标的转换 ##

这个转换就是前面所示的逆运算，如下图  
![图片来源于互联网][format_png 2]  
　　转换比较复杂，需要进行迭代运算。

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）下的高斯投影，其他的坐标系下的投影需要修改下参数。

## 总结 ##

至此，左右的转换已经完成。用了大概一周吧（中间有很长一段时间纠结于误差大的惊人，今天才知道甲方给的测试数据有问题），算是复习了下专业知识。  
　　[下载源代码][Link 2]

链接：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);

如果数据量比较多的话，后面这句可以放在一个循环里面。

[format_png]: /images/20220731/38653a73055c43cbb895f11a3fa0ffeb.png
[format_png 1]: /images/20220731/f55ad4a637e0473881265d3df2220d74.png
[Link 1]: https://me.csdn.net/qq_38876657
[format_png 2]: /images/20220731/34222a086b3b43b4a411cfad26225715.png
[Link 2]: http://download.csdn.net/detail/dahongdahong/9632754