经纬度转平面坐标xy并算距离代码

经纬度算直线距离通常有2种方法，比较常用的是第一种。

方法一：半正矢公式

第一种方法是通过半正矢公式推导
这种方法的优点是公式简单，计算快，缺点是地球并不是标准的球体，存在一定的误差

        public double dis_Haversine(double[] BL1, double[] BL2)
        {
            //B为纬度，L为精度
            double B1 = BL1[0] / 180 * Math.PI;
            double B2 = BL2[0] / 180 * Math.PI;
            double L1 = BL1[1] / 180 * Math.PI;
            double L2 = BL2[1] / 180 * Math.PI;
            double a = B1 - B2;
            double b = L1 - L2;
            double s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a / 2), 2) + Math.Cos(B1) * Math.Cos(B2) * Math.Pow(Math.Sin(b / 2), 2)));
            s = s * 6378137;
            s = Math.Round(s * 10000) / 10000;
            return s;
        }

方法二：先转成平面坐标再计算

先将经纬度坐标转为平面坐标（x，y），再根据勾股定理计算距离。
该种方法优点是计算误差较小，结果较为精确，缺点是计算公式较为复杂，原理是大地坐标公式，采用不同投影坐标系结果都会不同。

       public double dis_xy(double[] BL1, double[] BL2)//求两点距离函数
        {
            double[] r1 = BL2xy84(BL1, 121.464467);//121.464467是投影坐标系的中心纬度，不同城市不同
            double[] r2 = BL2xy84(BL2, 121.464467);
            double s = Math.Pow(Math.Pow((r1[0] - r2[0]), 2) + Math.Pow((r1[1] - r2[1]), 2), 0.5);
            return s;
        }
        public static double[] BL2xy84(double[] BL, double L0)//经纬度转换成平面坐标
        {
            double[] xy = new double[2];
            double X, B, L, a, b, e, e1, V, c, M, N, t, n, l, x, y;
            X = B2S(0, BL[0]);
            B = BL[0] / 180 * Math.PI;//纬度
            L = BL[1] / 180 * Math.PI;//精度
            L0 = L0 / 180 * Math.PI;
            a = 6378137.0; //WGS_84参考椭球参数
            // f=1/298.257223563;//椭球扁率
            b = 6356752.3142;//短轴
            e = (Math.Sqrt(a * a - b * b)) / a;
            e1 = (Math.Sqrt(a * a - b * b)) / b;
            V = Math.Sqrt(1 + (e1 * e1) * (Math.Cos(B)) * (Math.Cos(B)));
            c = (a * a) / b;
            M = c / (V * V * V);
            N = c / V;
            t = Math.Tan(B);
            n = Math.Sqrt((e1 * e1) * (Math.Cos(B)) * (Math.Cos(B)));
            l = L - L0;
            x = X + N * t * Math.Cos(B) * Math.Cos(B) * l * l * (0.5 + (1 / 24) * (5 - t * t + 9 * n * n + 4 * n * n * n * n) * Math.Cos(B) * Math.Cos(B) * l * l + 1 / 720 * (61 - 58 * t * t + t * t * t * t) * Math.Pow((Math.Cos(B)), 4) * l * l * l * l);
            y = N * Math.Cos(B) * l * (1 + 1 / 6 * (1 - t * t + n * n) * Math.Cos(B) * Math.Cos(B) * l * l + 1 / 120 * (5 - 18 * t * t + t * t * t * t + 14 * n * n - 58 * t * t * n * n) * Math.Pow((Math.Cos(B)), 4) * l * l * l * l);
            double[] d = { x, y };
            return d;

        }
        public static double B2S(double B1, double B2)//BL2xy84函数的内部一个函数
        {
            double a, b, e2, A, B, C, D, E, F, G, S;
            a = 6378137.0;//WGS_84参考椭球参数
            //f=1/298.257223563;%椭球扁率
            b = 6356752.3142;// 短轴
            e2 = (a * a - b * b) / (a * a);
            A = 1 + 3 * e2 / 4 + 45 * e2 * e2 / 64 + 175 * e2 * e2 * e2 / 256 + 11025 * Math.Pow(e2, 4) / 16384 + 43659 * Math.Pow(e2, 5) / 65536 + 693693 * Math.Pow(e2, 6) / 1048576;
            B = A - 1;
            C = 15 * e2 * e2 * e2 / 32 + 175 * e2 * e2 * e2 / 384 + 3675 * Math.Pow(e2, 4) / 8192 + 14553 * Math.Pow(e2, 5) / 32768 + 231231 * Math.Pow(e2, 6) / 524288;
            D = 35 * e2 * e2 * e2 / 96 + 735 * Math.Pow(e2, 4) / 2048 + 14553 * Math.Pow(e2, 5) / 40960 + 231231 * Math.Pow(e2, 6) / 655360;
            E = Math.Pow(e2, 4) * 315 / 1024 + Math.Pow(e2, 5) * 6237 / 20480 + Math.Pow(e2, 6) * 99099 / 327680;
            F = 693 * Math.Pow(e2, 5) / 2560 + 11011 * Math.Pow(e2, 6) / 40960;
            G = 1001 * Math.Pow(e2, 6) / 4096;
            B1 = B1 / 180 * Math.PI;
            B2 = B2 / 180 * Math.PI;
            S = a * (1 - e2) * (A * (B2 - B1) - B * (Math.Sin(B2) * Math.Cos(B2) - Math.Sin(B1) * Math.Cos(B1)) - C * (Math.Pow(Math.Sin(B2), 3) * Math.Cos(B2) - Math.Pow(Math.Sin(B1), 3) * Math.Cos(B1)) - D * (Math.Pow(Math.Sin(B2), 5) * Math.Cos(B2) - Math.Pow(Math.Sin(B1), 5) * Math.Cos(B1)) - E * (Math.Pow(Math.Sin(B2), 7) * Math.Cos(B2) - Math.Pow(Math.Sin(B1), 7) * Math.Cos(B1)) - F * (Math.Pow(Math.Sin(B2), 9) * Math.Cos(B2) - Math.Pow(Math.Sin(B1), 9) * Math.Cos(B1)) - G * (Math.Pow(Math.Sin(B2), 11) * Math.Cos(B2) - Math.Pow(Math.Sin(B1), 11) * Math.Cos(B1)));
            return S;
        }