这里的是没有给定参数所以没有用七参法,先将投影坐标系转了地理坐标系,然后做了横纵和比例纠偏
public static double[] Transform(PointF point) { double xParam = 94.362163134086399; double yParam = -310.26525523306055; double xMultiple = 1.19862910076924; double yMultiple = 1; //旋转中心点 double centerX = 114.00092403; double centerY = 36.14333070; //旋转角度 double Angle = 0.064894377180536; double X = Math.Round(point.X, 7) + 4000000; double Y = Math.Round(point.Y, 7) + 38500000; //double X = 60139 + 4000000; //double Y = 34944 + 38500000; // 由高斯投影坐标反算成经纬度 int ProjNo; int ZoneWide; ////带宽 double[] output = new double[2]; double longitude1, latitude1, longitude0, X0, Y0, xval, yval;//latitude0, double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI; iPI = 0.0174532925199433; ////3.1415926535898/180.0; a = 6378245.0; f = 1.0 / 298.3; //54年北京坐标系参数 //a = 6378140.0; f = 1 / 298.257; //80年西安坐标系参数 ZoneWide = 6; ////6度带宽 ProjNo = (int)(X / 1000000L); //查找带号 longitude0 = (ProjNo - 1) * ZoneWide + ZoneWide / 2; longitude0 = longitude0 * iPI; //中央经线 X0 = ProjNo * 1000000L + 500000L; Y0 = 0; xval = X - X0; 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) 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); 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); // 现状图 output[0] = longitude1 / iPI * xMultiple + xParam; output[1] = latitude1 / iPI * yMultiple + yParam; return output; }
横向有20m左右偏差,纵向有3m左右偏差。因为cad没有旋转,所以没有做旋转的纠偏。
这里的代码仅作参考,有更好建议可以提出来大家一起研究