大地坐标转公里网坐标

using System;

namespace CoordinateTransportation
{
    static class Program
    {
        static void Main()
        {
            do
            {
                Console.WriteLine("Please input your Ground Coordinate:");

                Console.Write("X:");
                var x = Console.ReadLine().Trim();
                Console.Write("Y:");
                var y = Console.ReadLine().Trim();
                Console.Write("Belt number:");
                var beltNumber = Console.ReadLine().Trim();

                CoordinateTransportation(double.Parse(x), double.Parse(y), double.Parse(beltNumber), out var gpsX, out var gpsY);
                
                var gcj02Coordinate = gps84_To_Gcj02(gpsY, gpsX);
                Console.WriteLine("The GPS Coordinate result of Ground Coordinate {0}, {1} is {2}, {3}", x, y, gcj02Coordinate[0], gcj02Coordinate[1]);
            } while (true);
        }

        private static double a = 6378245.0;
        private static double ee = 0.00669342162296594323;

        /// <summary>
        /// 坐标转换
        /// </summary>
        /// <param name="x">公里网横坐标 X</param>
        /// <param name="y">公里网纵坐标 Y</param>
        /// <param name="beltNumber">带号</param>
        /// <param name="longitude">大地坐标经度</param>
        /// <param name="latitude">大地坐标纬度</param>
        private static void CoordinateTransportation(double x, double y, double beltNumber, out double longitude, out double latitude)
        {
            int projNo; ////带宽
            double x0, y0;
            double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
            iPI = Math.PI / 180.0;
            //a = 6378245.0; f = 1.0 / 298.3; //BJ54 参数
            a = 6378140.0; f = 1 / 298.257; //XA80 参数
            a = 6378137.0; f = 1 / 298.257223563; // WGS84 参数
            projNo = (int)(x / 1000000L); //查找带号
            var longitude0 = beltNumber * 3;
            longitude0 = longitude0 * iPI; //中央经线
            x0 = projNo * 1000000L + 500000L;
            y0 = 0;
            var xval = x - x0; var yval = y - y0; //带内大地坐标
            e2 = 2 * f - f * f;
            e1 = (1.0 - Math.Sqrt(1 - e2)) / (1.0 + Math.Sqrt(1 - e2));
            ee = e2 / (1 - e2);
            M = yval;
            u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
            fai = u + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * u) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(
            4 * u)
            + (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * u) + (1097 * e1 * e1 * e1 * e1 / 512) * Math.Sin(8 * u);
            C = ee * Math.Cos(fai) * Math.Cos(fai);
            T = Math.Tan(fai) * Math.Tan(fai);
            NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(fai) * Math.Sin(fai));
            R = a * (1 - e2) / Math.Sqrt((1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin
            (fai) * Math.Sin(fai)));
            D = xval / NN;
            //计算经度(Longitude) 纬度(Latitude)
            var longitude1 = longitude0 + (D - (1 + 2 * T + C) * D * D * D / 6 + (5 - 2 * C + 28 * T - 3 * C * C + 8 * ee + 24 * T * T) * D
                                                                                                                                        * D * D * D * D / 120) / Math.Cos(fai);
            var latitude1 = fai - (NN * Math.Tan(fai) / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
                                                              + (61 + 90 * T + 298 * C + 45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
            //转换为度 DD
            longitude = longitude1 / iPI;
            latitude = latitude1 / iPI;
        }

        private static double[] gps84_To_Gcj02(double lat, double lon)
        {
            if (OutOfChina(lat, lon))
            {
                return new[] { lat, lon };
            }
            var dLat = TransformLat(lon - 105.0, lat - 35.0);
            var dLon = TransformLon(lon - 105.0, lat - 35.0);
            var radLat = lat / 180.0 * Math.PI;
            var magic = Math.Sin(radLat);
            magic = 1 - ee * magic * magic;
            var sqrtMagic = Math.Sqrt(magic);
            dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * Math.PI);
            dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * Math.PI);
            var mgLat = lat + dLat;
            var mgLon = lon + dLon;
            return new[] { mgLat, mgLon };
        }


        private static bool OutOfChina(double lat, double lon)
        {
            if (lon is < 72.004 or > 137.8347)
                return true;
            if (lat is < 0.8293 or > 55.8271)
                return true;
            return false;
        }

        private static double TransformLat(double x, double y)
        {
            var ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y
            + 0.2 * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * Math.PI) + 20.0 * Math.Sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(y * Math.PI) + 40.0 * Math.Sin(y / 3.0 * Math.PI)) * 2.0 / 3.0;
            ret += (160.0 * Math.Sin(y / 12.0 * Math.PI) + 320 * Math.Sin(y * Math.PI / 30.0)) * 2.0 / 3.0;
            return ret;
        }

        private static double TransformLon(double x, double y)
        {
            var ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1
            * Math.Sqrt(Math.Abs(x));
            ret += (20.0 * Math.Sin(6.0 * x * Math.PI) + 20.0 * Math.Sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
            ret += (20.0 * Math.Sin(x * Math.PI) + 40.0 * Math.Sin(x / 3.0 * Math.PI)) * 2.0 / 3.0;
            ret += (150.0 * Math.Sin(x / 12.0 * Math.PI) + 300.0 * Math.Sin(x / 30.0 * Math.PI)) * 2.0 / 3.0;
            return ret;
        }
    }
}