C#——GPS坐标与地方坐标相互转换

　　自决定总结之前做的模块以后，过了三天了，现在在公司，也不忙就简单的写一下。代码回去后附上。

　　这是老师布置的第一份任务，一开始也没理解这到底是怎么回事，后来就按照错误的理解去做，做了、错了；重新理解、又做、又错，又重新理解，大约花了三周时间才把问题搞明白。当时大二，做事效率是真不行，唉……

　　现在说正题吧，首先介绍一下任务：每个矿区有自己独立的坐标系，现在已知矿区中的三个点的两套坐标，一个是GPS坐标（经纬度坐标），另一个是矿区的独立坐标系中的坐标（X,Y坐标），根据这三个点的已知坐标，计算出用于坐标转换的参数。也就是给出一个点的GPS坐标，能够计算出对应的矿区坐标。

　　当初接到这个问题，也没有查资料，单凭自己的理解去做。最初最笨最傻的理解是把地球近似为圆的，然后呢矿区认为是平的。基于这个假设，根据简单的几何关系去计算，计算的几何方法就不说了……没用，白瞎忙活了好久。

　　后面就查了资料，这个问题前人已经总结了好多遍，现将坐标转换的具体步骤详细介绍:

    首先定义点结构体：

/*---------------------------------------------------------
 * 　 // Copyright (C) 2011 东北大学
　　　// 版权所有。
　　　//
　　　// 文件名：点的结构体
　　　// 文件功能描述：定义一个点的结构体
　　　//
　　　// 创建标识：  潘腾 2013.09.12
　　　//
　　　// 修改标识：
　　　// 修改描述：
 * 　　　------------------------------------------------------------*/
namespace OMCC.Calculate
{
    /// <summary>
    /// 定义一个点的结构体
    /// </summary>
    public struct APoint
    {
        public double b;      //经度坐标
        public double l;      //纬度坐标
        public double x1;     //高斯平面坐标系下的坐标
        public double y1;     //高斯平面坐标系下的坐标
        public double x2;     //地方坐标系下的坐标
        public double y2;     //地方坐标系下的坐标

        /// <summary>
        /// 点的结构体
        /// </summary>
        /// <param name="b"></param> 
        /// <param name="l"></param>
        /// <param name="x1"></param>
        /// <param name="y1"></param>
        /// <param name="x2"></param>
        /// <param name="y2"></param>
        public APoint(double b, double l, double x1, double y1,double x2,double y2)
        {
            this.b = b;
            this.l = l;
            this.x1 = x1;
            this.y1 = y1;
            this.x2 = x2;
            this.y2 = y2;
        }
    }
}

经纬度坐标到平面坐标之间的转换

/*---------------------------------------------------------
 * 　 // Copyright (C) 2011 东北大学
　　　// 版权所有。
　　　//
　　　// 文件名：高斯转换
　　　// 文件功能描述：由大地坐标（经纬度坐标）到平面直角坐标
 *                     （X Y坐标）之间的转换
　　　//
　　　// 创建标识：潘腾 2013.09.12
　　　//
　　　// 修改标识：
　　　// 修改描述：
 * 　　　------------------------------------------------------------*/
using System;

namespace OMCC.Calculate
{
    /*----------------------------------------------------
     * 输入：输入某个点的大地坐标
     * 输出：该点的高斯坐标
     *------------------------------------------------*/
    class BL2XY
    {
        #region  类属性
        private APoint PointA;

        public APoint A1
        {
            get { return PointA; }
            set { PointA = value; }
        }
        #endregion

        #region 高斯转换

        /// <summary>
        /// 
        /// </summary>
        /// <param name="PointA"></param>
        /// <returns></returns>
        public double[] BL2GS(APoint PointA)
        {
            //由弧度制到角度制在化为以秒为单位的转化
            PointA.b = PointA.b * 180 / Math.PI;
            PointA.b = PointA.b * 3600;
            PointA.l = PointA.l * 180 / Math.PI;
            PointA.l = PointA.l * 3600;

            double B = PointA.b;
            double L = PointA.l;
            double x, y;
            double L0;
            double l1;
            double p = 206264.80625;
            double b = B / p;
            double N;


            //N = (int)(L / (3 * 3600) + 0.5);
            //L0 =3 * N * 3600;
            //l1 = (L - L0) / p;

            N = (int)(L / (6 * 3600) + 1);
            L0 = (6 * N - 3) * 3600;
            l1 = (L - L0) / p;

            double c = Math.Pow(Math.Cos(b), 2);
            double c1 = Math.Sin(b) * Math.Cos(b);
            double c2 = Math.Cos(b);
            double l2 = Math.Pow(l1, 2);
            double n = 6399698.902 - (21562.267 - (108.973 - 0.612 * c) * c) * c;
            double a0 = 32140.404 - (135.3302 - (0.7092 - 0.0040 * c) * c) * c;
            double a4 = (0.25 + 0.00252 * c) * c - 0.04166;
            double a6 = (0.166 * c - 0.084) * c;
            double a3 = (0.3333333 + 0.001123 * c) * c - 0.1666667;
            double a5 = 0.0083 - (0.1667 - (0.1968 + 0.0040 * c) * c) * c;


            x = 6367558.4969 * b - (a0 - (0.5 + (a4 + a6 * l2) * l2) * l2 * n) * c1;
            y = (1 + (a3 + a5 * l2) * l2) * l1 * n * c2;

            double y1 = y + 500000.0;
            double abc = new double();
            for (int i = 1; y1 / i > 1; i = i * 10)
            {
                abc = (N * i * 10 + y1);
            }
            return new double[] { x, abc };
        }
        #endregion
    }
}

 两个平面坐标之间的转换

/*---------------------------------------------------------
 * 　 // Copyright (C) 2011 东北大学
　　　// 版权所有。
　　　//
　　　// 文件名：平面坐标转换参数计算
　　　// 文件功能描述：返回两个平面直角坐标变换之间的参数
　　　//
　　　// 创建标识：潘腾 2013.09.12
　　　//
　　　// 修改标识：
　　　// 修改描述：
 * 　　　------------------------------------------------------------*/
using System.Collections.Generic;

namespace OMCC.Calculate
{
    /*----------------------------------------------------
     * 输入：三个点A B C在两个坐标系中的坐标
     * 输出：两个坐标系的转换参数
     *------------------------------------------------*/
    class XY2XY
    {
        #region  类属性

        private APoint pointA;
        /// <summary>
        /// 第一个点
        /// </summary>
        public APoint PointA
        {
            get { return pointA; }
            set { pointA = value; }
        }

        private APoint pointB;
        /// <summary>
        /// 第二个点
        /// </summary>
        public APoint PointB
        {
            get { return pointB; }
            set { pointB = value; }
        }
        /// <summary>
        /// 第三个点
        /// </summary>
        private APoint pointC;

        public APoint PointC
        {
            get { return pointC; }
            set { pointC = value; }
        }

        #endregion

        #region 返回两个平面坐标转换参数
        public List<double> GetXY2XY()
        {
            double ratioa3 = (((pointB.x1) - (pointA.x1)) * ((pointC.x2) - (pointA.x2)) - ((pointB.x2) 
                - (pointA.x2)) * ((pointC.x1) - (pointA.x1))) /(((pointC.y1) - (pointA.y1)) * ((pointB.x1)
                - (pointA.x1)) - ((pointC.x1) - (pointA.x1)) * ((pointB.y1) - (pointA.y1)));
            double ratioa2 = ((pointC.x2) - (pointA.x2) - ratioa3 * ((pointC.y1) - (pointA.y1))) / ((pointC.x1) - (pointA.x1));
            double ratioa1 = (pointA.x2) - ratioa2 * (pointA.x1) - ratioa3 * (pointA.y1);
            double ratiob3 = (((pointB.x1) - (pointA.x1)) * ((pointC.y2) - (pointA.y2)) - ((pointB.y2) - (pointA.y2)) * 
                ((pointC.x1) - (pointA.x1))) /(((pointC.y1) - (pointA.y1)) * ((pointB.x1) - (pointA.x1))
                - ((pointC.x1) - (pointA.x1)) * ((pointB.y1) - (pointA.y1)));
            double ratiob2 = ((pointC.y2) - (pointA.y2) - ratiob3 * ((pointC.y1) - (pointA.y1))) / ((pointC.x1) - (pointA.x1));
            double ratiob1 = (pointA.y2) - ratiob2 * (pointA.x1) - ratiob3 * (pointA.y1);
            List<double> result = new List<double>();
            result.Add(ratioa1);
            result.Add(ratioa2);
            result.Add(ratioa3);
            result.Add(ratiob1);
            result.Add(ratiob2);
            result.Add(ratiob3);
            return result;
        }
        #endregion
    }
}

OK了，就写到这儿吧，时间过太久了，也不想在花太多时间去回顾了。写的也许不够详细………………