using ESRI.ArcGIS.Client.Geometry; using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.IO; using System.Linq; using System.Runtime.Serialization.Json; using System.Text; namespace Kingo.Plugin.MakeTaskPackage.Common { public class CoordinateTrans { #region 四参数集合操作 /// /// 初始化四参数 /// /// public static void InitParams() { try { string selectSql = "SELECT * FROM TB_CSSZ"; DataTable dt = null;// Kingo.SanDiaoCheck.Common.CommonAPI.GetSystemTable(selectSql); if (dt != null && dt.Rows.Count > 0) { ParamsEntity hcryEntity = null; foreach (DataRow item in dt.Rows) { hcryEntity = new ParamsEntity(); hcryEntity.BSM = ObjectConvertOperate.ObjectToInt(item["BSM"]); hcryEntity.XZQDM = ObjectConvertOperate.ObjectToString(item["XZQDM"]); hcryEntity.X80OFFSET = ObjectConvertOperate.ObjectToDouble(item["X80OFFSET"]); hcryEntity.X80ROTATE = ObjectConvertOperate.ObjectToDouble(item["X80ROTATE"]); hcryEntity.X84OFFSET = ObjectConvertOperate.ObjectToDouble(item["X84OFFSET"]); hcryEntity.X84ROTATE = ObjectConvertOperate.ObjectToDouble(item["X84ROTATE"]); hcryEntity.Y80OFFSET = ObjectConvertOperate.ObjectToDouble(item["Y80OFFSET"]); hcryEntity.Y80ROTATE = ObjectConvertOperate.ObjectToDouble(item["Y80ROTATE"]); hcryEntity.Y84OFFSET = ObjectConvertOperate.ObjectToDouble(item["Y84OFFSET"]); hcryEntity.Y84ROTATE = ObjectConvertOperate.ObjectToDouble(item["Y84ROTATE"]); hcryEntity.DH = ObjectConvertOperate.ObjectToString(item["DH"]); AddOrUpdateParams(hcryEntity); } } } catch (Exception ex) { } } /// /// 四参数集合 /// private static Dictionary paramsEntityList = new Dictionary(); /// /// 增加或修改四参数 /// /// /// public static bool AddOrUpdateParams(ParamsEntity paramEntity) { lock (paramsEntityList) { if (paramsEntityList.ContainsKey(paramEntity.XZQDM)) paramsEntityList[paramEntity.XZQDM] = paramEntity; else paramsEntityList.Add(paramEntity.XZQDM, paramEntity); } return true; } /// /// 获取四参数根据行政区 /// /// /// public static ParamsEntity GetParamsByXzq(string xzq) { if (string.IsNullOrWhiteSpace(xzq)) return null; if (paramsEntityList.ContainsKey(xzq)) return paramsEntityList[xzq]; return null; } public static Param4 Get80To84(string xzq) { ParamsEntity paramsEntity = GetParamsByXzq(xzq); Param4 param4 = null; if (paramsEntity != null) { if (paramsEntity.X80OFFSET == 0 && paramsEntity.X80ROTATE == 0 && paramsEntity.Y80OFFSET == 0 && paramsEntity.Y80ROTATE == 0) return null; param4 = new Param4(); param4.PX = paramsEntity.X80OFFSET; param4.PY = paramsEntity.Y80OFFSET; param4.SX = paramsEntity.X80ROTATE; param4.SY = paramsEntity.Y80ROTATE; } return param4; } public static Param4 Get84To80(string xzq) { ParamsEntity paramsEntity = GetParamsByXzq(xzq); Param4 param4 = null; if (paramsEntity != null) { if (paramsEntity.X84OFFSET == 0 && paramsEntity.X84ROTATE == 0 && paramsEntity.Y84OFFSET == 0 && paramsEntity.Y84ROTATE == 0) return null; param4 = new Param4(); param4.PX = paramsEntity.X84OFFSET; param4.PY = paramsEntity.Y84OFFSET; param4.SX = paramsEntity.X84ROTATE; param4.SY = paramsEntity.Y84ROTATE; } return param4; } #endregion /// /// 计算四参数 /// /// 重合点个数 /// 重合点在旧坐标系中的坐标 /// 重合点的在新坐标系中的坐标 /// private static double[][] Calculate4Param(int N, double[][] oldp, double[][] newp1) { int i, j; double[][] coneB; //系数矩阵B coneB = new double[2 * N][]; for (i = 0; i < 2 * N; i++) coneB[i] = new double[4]; double[][] TconeB; //系数矩阵B的转置 TconeB = new double[4][]; for (i = 0; i < 4; i++) TconeB[i] = new double[2 * N]; double[][] consL; //常数项矩阵L consL = new double[2 * N][]; for (i = 0; i < 2 * N; i++) consL[i] = new double[1]; double[][] Nbb; //法方程系数矩阵Nbb Nbb = new double[4][]; for (i = 0; i < 4; i++) Nbb[i] = new double[4]; double[][] DNbb; //法方程系数矩阵Nbb的逆矩阵 DNbb = new double[4][]; for (i = 0; i < 4; i++) DNbb[i] = new double[4]; double[][] ENbb; //法方程系数矩阵Nbb的扩展矩阵 ENbb = new double[4][]; for (i = 0; i < 4; i++) ENbb[i] = new double[8]; double[][] W; //用于表示Nbb*TconeB W = new double[4][]; for (i = 0; i < 4; i++) W[i] = new double[2 * N]; double[][] parameter; //平差参数 parameter = new double[4][]; for (i = 0; i < 4; i++) parameter[i] = new double[1]; double[][] V; //平差改正数 V = new double[2 * N][]; for (i = 0; i < 2 * N; i++) V[i] = new double[1]; double[][] result2; //用于存储系数B与平差参数parameter的乘积 result2 = new double[2 * N][]; for (i = 0; i < 2 * N; i++) result2[i] = new double[1]; for (i = 0; i < 2 * N; i++) //写出方程系数矩阵B { for (j = 0; j < 4; j++) { coneB[i][0] = (0.5 + Math.Pow(-1, i) * 0.5); coneB[i][1] = (0.5 - Math.Pow(-1, i) * 0.5); coneB[i][2] = oldp[i / 2][0] * (0.5 + Math.Pow(-1, i) * 0.5) + oldp[i / 2][1] * (0.5 - Math.Pow(-1, i) * 0.5); coneB[i][3] = oldp[i / 2][0] * (0.5 - Math.Pow(-1, i) * 0.5) + oldp[i / 2][1] * (-0.5 - Math.Pow(-1, i) * 0.5); } } for (i = 0; i < 2 * N; i++) //写出常数项矩阵L for (j = 0; j < 1; j++) consL[i][j] = newp1[i / 2][0] * (0.5 + Math.Pow(-1, i) * 0.5) + newp1[i / 2][1] * (0.5 - Math.Pow(-1, i) * 0.5); transform(coneB, TconeB, 2 * N, 4); //调用转置函数 multiplication(TconeB, coneB, Nbb, 4, 2 * N, 4); //调用乘法函数计算法方程系数 det(Nbb, ENbb, DNbb, 4); //调用求逆函数,计算法方程系数的逆阵 multiplication(TconeB, consL, W, 4, 2 * N, 1);//调用乘法函数,计算常数项W multiplication(DNbb, W, parameter, 4, 4, 1); //调用乘法函数,解算转换系数a,b,c,d return parameter; } /// /// 四参数转换 /// /// /// /// /// private static double[][] Tras4Param(double[][] oldpoint, double[][] parameter, int M) { int i; double[][] newpoint; //用于存储旧坐标经转换后的新坐标 newpoint = new double[M][]; for (i = 0; i < M; i++) { newpoint[i] = new double[2]; } for (i = 0; i < M; i++) { newpoint[i][0] = parameter[0][0] + oldpoint[i][0] * parameter[2][0] - oldpoint[i][1] * parameter[3][0]; newpoint[i][1] = parameter[1][0] + oldpoint[i][1] * parameter[2][0] + oldpoint[i][0] * parameter[3][0]; } return newpoint; } //求逆函数 private static void det(double[][] p, double[][] b, double[][] c, int n) { int i, j, k; double m; //扩展为增广矩阵 for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { b[i][j] = p[i][j]; } for (j = n; j < 2 * n; j++) { if (i + n == j) b[i][j] = 1; else b[i][j] = 0; } } //将左半部分化为单位矩阵 for (i = 0; i < n; i++) { m = b[i][i]; for (j = i; j < n * 2; j++) b[i][j] /= m; for (j = i + 1; j < n; j++) { m = b[j][i]; for (k = i; k < n * 2; k++) b[j][k] = b[j][k] - m * b[i][k]; } } for (i = n - 1; i >= 0; i--) { m = b[i][i]; for (j = i - 1; j >= 0; j--) { m = b[j][i]; for (k = i; k < n * 2; k++) b[j][k] = b[j][k] - m * b[i][k]; } } for (i = 0; i < n; i++) { for (j = n; j < 2 * n; j++) { c[i][j - n] = b[i][j]; } } } //转置函数 private static void transform(double[][] a, double[][] b, int m, int n) { int i, j; for (i = 0; i < m; i++) for (j = 0; j < n; j++) { b[j][i] = a[i][j]; } } //乘法函数 public static void multiplication(double[][] a, double[][] b, double[][] c, int m, int s, int n) { int i, j, k; double sum; for (i = 0; i < m; i++) for (j = 0; j < n; j++) { sum = 0.0; for (k = 0; k < s; k++) sum += a[i][k] * b[k][j]; c[i][j] = sum; } } //矩阵加法函数 private static void add(double[][] a, double[][] b, double[][] ab, int m, int n, int sign) { for (int i = 0; i < m; i++) for (int j = 0; j < n; j++) ab[i][j] = a[i][j] + Math.Pow(-1, sign) * b[i][j]; //sign为偶数是表示矩阵相加,奇数表示相减 } /// /// 根据同位点计算平面转换四参数 /// /// /// /// public static Param4 Calculate4Param(MapPoint[] sourcePoints, MapPoint[] targetPoints) { if (sourcePoints == null || targetPoints == null || sourcePoints.Length != targetPoints.Length || sourcePoints.Length == 0) throw new ArgumentException("参数为空或者个数不一致", "sourcePoints,targetPoints"); int n = sourcePoints.Length; double[][] oldp1 = new double[n][]; double[][] newp1 = new double[n][]; for (int i = 0; i < n; i++) { oldp1[i] = new double[] { sourcePoints[i].X, sourcePoints[i].Y }; newp1[i] = new double[] { targetPoints[i].X, targetPoints[i].Y }; } double[][] param4 = Calculate4Param(n, oldp1, newp1); return new Param4() { PX = param4[0][0], PY = param4[1][0], SX = param4[2][0], SY = param4[3][0] }; } public static MapPoint[] Tras4Param(MapPoint[] sourcePoints, Param4 param4) { if (sourcePoints == null || sourcePoints.Length == 0) throw new ArgumentNullException("sourcePoints"); if (param4 == null) throw new ArgumentNullException("param4"); double[][] parameter = new double[4][]; parameter[0] = new double[] { param4.PX }; parameter[1] = new double[] { param4.PY }; parameter[2] = new double[] { param4.SX }; parameter[3] = new double[] { param4.SY }; int m = sourcePoints.Length; double[][] oldp = new double[m][]; for (int i = 0; i < m; i++) { oldp[i] = new double[] { sourcePoints[i].X, sourcePoints[i].Y }; } double[][] newp = Tras4Param(oldp, parameter, m); MapPoint[] targetPoints = new MapPoint[m]; for (int i = 0; i < m; i++) { targetPoints[i] = new MapPoint(newp[i][0], newp[i][1]); } return targetPoints; } /// /// 平面四参数转换方法 /// /// /// /// public static MapPoint Tras4Param(MapPoint sourcePoint, Param4 param4) { return Tras4Param(new MapPoint[] { sourcePoint }, param4)[0]; } //高斯投影正、反算 //////3度带宽 //高斯投影由经纬度(Unit:DD)反算大地坐标(含带号,Unit:Metres) private static void GaussProjCal(double longitude, double latitude, double a, double f, out double X, out double Y, int ProjNo = 0) { //int ProjNo = 0; int ZoneWide; ////带宽 double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval; double e2, ee, NN, T, C, A, M, iPI; iPI = Math.PI / 180; ////3.1415926535898/180.0; ZoneWide = 3; ////3度带宽 //a = 6378245.0; f = 1.0 / 298.3; //54年北京坐标系参数 //a=6378140.0; f=1/298.257; //80年西安坐标系参数 //a=6378137m;f=1/298.257223563;//WGS-84坐标系 if (ProjNo == 0) { ProjNo = (int)Math.Round(longitude / ZoneWide); } longitude0 = ProjNo * ZoneWide; longitude0 = longitude0 * iPI; latitude0 = 0; longitude1 = longitude * iPI; //经度转换为弧度 latitude1 = latitude * iPI; //纬度转换为弧度 e2 = 2 * f - f * f; ee = e2 * (1.0 - e2); NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(latitude1) * Math.Sin(latitude1)); T = Math.Tan(latitude1) * Math.Tan(latitude1); C = ee * Math.Cos(latitude1) * Math.Cos(latitude1); A = (longitude1 - longitude0) * Math.Cos(latitude1); M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.Sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.Sin(6 * latitude1)); xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120); yval = M + NN * Math.Tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720); X0 = 1000000L * (ProjNo) + 500000L; Y0 = 0; xval = xval + X0; yval = yval + Y0; X = xval; Y = yval; } /// /// 高斯投影正算 /// /// /// /// public static MapPoint GaussProjCal(MapPoint sourcePoint, EarthParam param, int ProjNo = 0) { double x, y; GaussProjCal(sourcePoint.X, sourcePoint.Y, param.A, param.F, out x, out y, ProjNo); return new MapPoint(x, y); } //高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD) private static void GaussProjInvCal(double X, double Y, double a, double f, out double longitude, out double latitude, int ProjNo = 0) { //int ProjNo; int ZoneWide; ////带宽 double longitude1, latitude1, longitude0, latitude0, X0, Y0, xval, yval; double e1, e2, ee, NN, T, C, M, D, R, u, fai, iPI; iPI = Math.PI / 180; ////3.1415926535898/180.0; //a = 6378245.0; f = 1.0 / 298.3; //54年北京坐标系参数 //a=6378140.0; f=1/298.257; //80年西安坐标系参数 //a=6378137m;f=1/298.257223563;//WGS-84坐标系 ZoneWide = 3; ////3度带宽 if (ProjNo == 0) { ProjNo = (int)(X / 1000000L); //查找带号 } longitude0 = ProjNo * ZoneWide; 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); //转换为度 DD longitude = longitude1 / iPI; latitude = latitude1 / iPI; } /// /// 高斯投影反算 /// /// /// /// public static MapPoint GaussProjInvCal(MapPoint sourcePoint, EarthParam param, int ProjNo = 0) { double x, y; GaussProjInvCal(sourcePoint.X, sourcePoint.Y, param.A, param.F, out x, out y, ProjNo); return new MapPoint(x, y); } /// /// 西安80平面坐标转84经纬度坐标 /// /// /// /// public static MapPoint Xian80ToWgs84(MapPoint sourcePoint, Param4 param = null, int ProjNo = 0) { if (sourcePoint == null) throw new ArgumentNullException("sourcePoint"); if (param != null && param.SX != 0 && param.PX != 0) { MapPoint targetPoint = Tras4Param(sourcePoint, param); return GaussProjInvCal(targetPoint, EarthParam.WGS84, ProjNo); } else { return GaussProjInvCal(sourcePoint, EarthParam.WGS84, ProjNo); } } /// /// 84经纬度坐标转西安80平面坐标 /// /// /// /// public static MapPoint Wgs84ToXian80(MapPoint sourcePoint, Param4 param = null, int ProjNo = 0) { if (sourcePoint == null) { throw new ArgumentNullException("sourcePoint"); } MapPoint targetPoint = GaussProjCal(sourcePoint, EarthParam.WGS84, ProjNo); if (param != null && param.PX != 0 && param.SX != 0) { return Tras4Param(targetPoint, param); } else { return targetPoint; } } public static MapPoint Wgs84ToXian80(MapPoint sourcePoint, string xzqdm, int ProjNo = 0) { return Wgs84ToXian80(sourcePoint, Get84To80(xzqdm), ProjNo); } public static Envelope Xian80ToWgs84(Envelope source, Param4 param = null, int ProjNo = 0) { if (source == null) { throw new ArgumentNullException("source"); } if (ProjNo == 0) { MapPoint pt = source.GetCenter(); ProjNo = (int)(pt.X / 1000000L); } MapPoint ptLT = Xian80ToWgs84(new MapPoint(source.XMin, source.YMax), param, ProjNo); MapPoint ptRB = Xian80ToWgs84(new MapPoint(source.XMax, source.YMin), param, ProjNo); return new Envelope(ptLT.X, ptRB.Y, ptRB.X, ptLT.Y); } public static Envelope Wgs84ToXian80(Envelope source, Param4 param = null, int ProjNo = 0) { if (source == null) throw new ArgumentNullException("source"); if (ProjNo == 0) { MapPoint pt = source.GetCenter(); ProjNo = (int)Math.Round(pt.X / 3); } MapPoint ptLT = Wgs84ToXian80(new MapPoint(source.XMin, source.YMax), param, ProjNo); MapPoint ptRB = Wgs84ToXian80(new MapPoint(source.XMax, source.YMin), param, ProjNo); return new Envelope(ptLT.X, ptRB.Y, ptRB.X, ptLT.Y); } /// /// 度分秒转成度 /// /// /// public static MapPoint DFM2D(MapPoint point) { return new MapPoint(DFM2D(point.X), DFM2D(point.Y)); } /// /// 度转成度分秒 /// /// /// public static MapPoint D2DFM(MapPoint point) { return new MapPoint(D2DFM(point.X), D2DFM(point.Y)); } public static double DFM2D(double dfm) { double d = Math.Floor(dfm); double f = Math.Floor((dfm - d) * 100); double m = (dfm - d - f / 100) * 10000; return d + f / 60 + m / 3600; } public static double D2DFM(double degree) { double d = Math.Floor(degree); double f = Math.Floor((degree - d) * 60); double m = (degree - d - f / 60) * 3600; return d + f / 100 + m / 10000; } public static int Get80PrjNo(double x) { if (x > 1000000 && x < 70000000) { return (int)x / 1000000; } return 0; } } public class Param4 : INotifyPropertyChanged { public event PropertyChangedEventHandler PropertyChanged; protected void PropertyChangedNotify(string propName) { if (this.PropertyChanged != null) PropertyChanged(this, new PropertyChangedEventArgs(propName)); } private double _PX; public double PX { get { return _PX; } set { _PX = value; PropertyChangedNotify("PX"); } } private double _PY; public double PY { get { return _PY; } set { _PY = value; PropertyChangedNotify("PY"); } } private double _SX; public double SX { get { return _SX; } set { _SX = value; PropertyChangedNotify("SX"); } } private double _SY; public double SY { get { return _SY; } set { _SY = value; PropertyChangedNotify("SY"); } } public Param4 Reverse() { return new Param4() { PX = -this.PX, PY = -this.PY, SX = -this.SX, SY = -this.SY }; } public bool IsEmpty { get { return PX == 0 && PY == 0 && SX == 0 && SY == 0; } } } public class EarthParam { /// /// 长轴 eg:6378137 /// public double A { get; set; } /// /// 扁心率 eg:1 / 298.257223563 /// public double F { get; set; } public static EarthParam WGS84 { get; private set; } public static EarthParam XIAN80 { get; private set; } public static EarthParam BJ54 { get; private set; } static EarthParam() { WGS84 = new EarthParam() { A = 6378137, F = 1 / 298.257223563 }; XIAN80 = new EarthParam() { A = 6378140, F = 1.0 / 298.3 }; BJ54 = new EarthParam() { A = 6378245.0, F = 1 / 298.257 }; } } public class YPRMethod { public static double GetAngle(double yaw, double pitch, double roll) { yaw = yaw * Math.PI / 180; pitch = pitch * Math.PI / 180; roll = roll * Math.PI / 180; //MapPoint pt1 = new MapPoint(0, 0, 0, 0); //MapPoint pt2 = new MapPoint(0, 1, 0, 0); //MapPoint pt22 = new MapPoint(0, 1, 0, 0); //MapPoint pt3 = new MapPoint(1, 0, 0, 0); //rotate(pt1, yaw, pitch, roll); //rotate(pt2, yaw, pitch, roll); //rotate2(pt22, yaw, pitch, roll); //rotate(pt3, yaw, pitch, roll); //MapPoint vec = get_Normal(pt1, pt2, pt3); //if (vec.X == 0 && vec.Y == 0) //{ // return (360 + yaw / Math.PI * 180) % 360; //} //double angle = Math.Atan2(vec.Y, vec.X); //if (double.IsNaN(angle)) //{ // return angle; //} //else //{ // return (360 + angle * 180 / Math.PI) % 360; //} MapPoint pt = RotateOL(yaw, pitch, roll); double angle = Math.Atan2(pt.X, pt.Z); if (double.IsNaN(angle) || angle == 0) { return yaw * 180 / Math.PI; } else { return -angle * 180 / Math.PI; } //return 0; } private static void rotate(MapPoint pt, double yaw, double pitch, double roll) { double x1 = pt.X * Math.Cos(yaw) - pt.Y * Math.Sin(yaw); double y1 = pt.X * Math.Sin(yaw) + pt.Y * Math.Cos(yaw); double z1 = pt.Z; double y2 = y1 * Math.Cos(roll) - z1 * Math.Sin(roll); double z2 = y1 * Math.Sin(roll) + z1 * Math.Cos(roll); double x2 = x1; pt.Z = z2 * Math.Cos(pitch) - x2 * Math.Sin(pitch); pt.X = z2 * Math.Sin(pitch) + x2 * Math.Cos(pitch); pt.Y = y2; } private static void rotate2(MapPoint pt, double yaw, double pitch, double roll) { double y1 = pt.Y * Math.Cos(roll) - pt.Z * Math.Sin(roll); double z1 = pt.Y * Math.Sin(roll) + pt.Z * Math.Cos(roll); double x1 = pt.X; double x2 = x1 * Math.Cos(yaw) - y1 * Math.Sin(yaw); double y2 = x1 * Math.Sin(yaw) + y1 * Math.Cos(yaw); double z2 = z1; pt.Z = z2 * Math.Cos(pitch) - x2 * Math.Sin(pitch); pt.X = z2 * Math.Sin(pitch) + x2 * Math.Cos(pitch); pt.Y = y2; } private static MapPoint get_Normal(MapPoint p1, MapPoint p2, MapPoint p3) { double a = ((p2.Y - p1.Y) * (p3.Z - p1.Z) - (p2.Z - p1.Z) * (p3.Y - p1.Y)); double b = ((p2.Z - p1.Z) * (p3.X - p1.X) - (p2.X - p1.X) * (p3.Z - p1.Z)); double c = ((p2.X - p1.X) * (p3.Y - p1.Y) - (p2.Y - p1.Y) * (p3.X - p1.X)); return new MapPoint(a, b, c, 0); } private static MapPoint RotateOL(double yaw, double pitch, double roll) { //double[][] mb = GetMaxtrixB(roll); //double[][] mc = GetMaxtrixC(pitch); //double[][] md = GetMaxtrixD(yaw); double[][] ma = GetMaxtrixA(pitch, yaw, roll); double[][] xyz = new double[3][]; xyz[0] = new double[1] { 0 }; xyz[1] = new double[1] { 1 }; xyz[2] = new double[1] { 0 }; //double[][] bc = new double[3][]; //for (int i = 0; i < 3; i++) //{ // bc[i] = new double[3]; //} //CoordinateTrans.multiplication(mb, mc, bc, 3, 3, 3); //double[][] bcd = new double[3][]; //for (int i = 0; i < 3; i++) //{ // bcd[i] = new double[3]; //} //CoordinateTrans.multiplication(bc, md, bcd, 3, 3, 3); double[][] xyz1 = new double[3][]; for (int i = 0; i < 3; i++) { xyz1[i] = new double[1]; } CoordinateTrans.multiplication(ma, xyz, xyz1, 3, 3, 1); return new MapPoint(xyz1[0][0], xyz1[1][0], xyz1[2][0], 0); } private static double[][] GetMaxtrixB(double d) { double[][] mxtrix = new double[3][]; mxtrix[0] = new double[] { Math.Cos(d), Math.Sin(d), 0 }; mxtrix[1] = new double[] { -Math.Sin(d), Math.Cos(d), 0 }; mxtrix[2] = new double[] { 0, 0, 1 }; return mxtrix; } private static double[][] GetMaxtrixC(double d) { double[][] mxtrix = new double[3][]; mxtrix[0] = new double[] { 1, 0, 0 }; mxtrix[1] = new double[] { 0, Math.Cos(d), Math.Sin(d) }; mxtrix[2] = new double[] { 0, -Math.Sin(d), Math.Cos(d) }; return mxtrix; } private static double[][] GetMaxtrixD(double d) { double[][] mxtrix = new double[3][]; mxtrix[0] = new double[] { Math.Cos(d), Math.Sin(d), 0 }; mxtrix[1] = new double[] { -Math.Sin(d), Math.Cos(d), 0 }; mxtrix[2] = new double[] { 0, 0, 1 }; return mxtrix; } private static double[][] GetMaxtrixA(double x, double y, double z) { double[][] mxtrix = new double[3][]; mxtrix[0] = new double[] { Math.Cos(y)*Math.Cos(z)-Math.Sin(x)*Math.Sin(y)*Math.Sin(z), Math.Cos(y)*Math.Sin(z)+Math.Sin(x)*Math.Sin(y)*Math.Cos(z), -Math.Sin(y)*Math.Cos(x) }; mxtrix[1] = new double[] { -Math.Cos(x)*Math.Sin(z), Math.Cos(x)*Math.Cos(z), Math.Sin(x) }; mxtrix[2] = new double[] { Math.Sin(y)*Math.Cos(z)+Math.Sin(x)*Math.Cos(y)*Math.Sin(z), Math.Sin(y)*Math.Sin(z)+Math.Sin(x)*Math.Cos(y)*Math.Cos(z), Math.Cos(x)*Math.Cos(y) }; return mxtrix; } } /// /// 墨卡托与经纬度转换操作 /// public class MercatorToLonLat { /// /// 经纬度转墨卡托坐标 /// /// /// public static ESRI.ArcGIS.Client.Geometry.MapPoint lonLat2Mercator(ESRI.ArcGIS.Client.Geometry.MapPoint lonLat) { double x1 = 0, y1 = 0; transform(lonLat.Y, lonLat.X, out y1, out x1); lonLat = new MapPoint(x1, y1); ESRI.ArcGIS.Client.Geometry.MapPoint mercator = new ESRI.ArcGIS.Client.Geometry.MapPoint(); double x = lonLat.X * 20037508.34 / 180; double y = Math.Log(Math.Tan((90 + lonLat.Y) * Math.PI / 360)) / (Math.PI / 180); y = y * 20037508.34 / 180; mercator.X = x; mercator.Y = y; return mercator; } /// /// 墨卡托坐标转经纬度 /// /// /// public static ESRI.ArcGIS.Client.Geometry.MapPoint Mercator2lonLat(ESRI.ArcGIS.Client.Geometry.MapPoint mercator) { ESRI.ArcGIS.Client.Geometry.MapPoint lonLat = new ESRI.ArcGIS.Client.Geometry.MapPoint(); double x = mercator.X / 20037508.34 * 180; double y = mercator.Y / 20037508.34 * 180; y = 180 / Math.PI * (2 * Math.Atan(Math.Exp(y * Math.PI / 180)) - Math.PI / 2); lonLat.X = x; lonLat.Y = y; return lonLat; } public static ESRI.ArcGIS.Client.Geometry.Envelope Mercator2lonLat(ESRI.ArcGIS.Client.Geometry.Envelope mercator) { double xMax = MercatorX2Lon(mercator.XMax); double xMin = MercatorX2Lon(mercator.XMin); double yMax = MercatorY2Lat(mercator.YMax); double yMin = MercatorY2Lat(mercator.YMin); return new Envelope() { XMax = xMax, XMin = xMin, YMax = yMax, YMin = yMin }; } public static ESRI.ArcGIS.Client.Geometry.Envelope LonLat2Mercator(ESRI.ArcGIS.Client.Geometry.Envelope lonlat) { double xMax = Lon2Mercator(lonlat.XMax); double xMin = Lon2Mercator(lonlat.XMin); double yMax = Lat2Mercator(lonlat.YMax); double yMin = Lat2Mercator(lonlat.YMin); return new Envelope() { XMax = xMax, XMin = xMin, YMax = yMax, YMin = yMin }; } private static double MercatorX2Lon(double x) { return x / 20037508.34 * 180; } private static double MercatorY2Lat(double y) { double lat = y / 20037508.34 * 180; return 180 / Math.PI * (2 * Math.Atan(Math.Exp(lat * Math.PI / 180)) - Math.PI / 2); } private static double Lon2Mercator(double x) { return x * 20037508.34 / 180; } private static double Lat2Mercator(double y) { double mercatorY = Math.Log(Math.Tan((90 + y) * Math.PI / 360)) / (Math.PI / 180); return mercatorY * 20037508.34 / 180; } #region 经纬度转火星坐标 const double pi = 3.14159265358979324; // // Krasovsky 1940 // // a = 6378245.0, 1/f = 298.3 // b = a * (1 - f) // ee = (a^2 - b^2) / a^2; const double a = 6378245.0; const double ee = 0.00669342162296594323; // // World Geodetic System ==> Mars Geodetic System public static void transform(double wgLat, double wgLon, out double mgLat, out double mgLon) { if (outOfChina(wgLat, wgLon)) { mgLat = wgLat; mgLon = wgLon; return; } double dLat = transformLat(wgLon - 105.0, wgLat - 35.0); double dLon = transformLon(wgLon - 105.0, wgLat - 35.0); double radLat = wgLat / 180.0 * pi; double magic = Math.Sin(radLat); magic = 1 - ee * magic * magic; double sqrtMagic = Math.Sqrt(magic); dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi); dLon = (dLon * 180.0) / (a / sqrtMagic * Math.Cos(radLat) * pi); mgLat = wgLat + dLat; mgLon = wgLon + dLon; } static bool outOfChina(double lat, double lon) { if (lon < 72.004 || lon > 137.8347) return true; if (lat < 0.8293 || lat > 55.8271) return true; return false; } static double transformLat(double x, double y) { double 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 * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0; ret += (20.0 * Math.Sin(y * pi) + 40.0 * Math.Sin(y / 3.0 * pi)) * 2.0 / 3.0; ret += (160.0 * Math.Sin(y / 12.0 * pi) + 320 * Math.Sin(y * pi / 30.0)) * 2.0 / 3.0; return ret; } static double transformLon(double x, double y) { double 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 * pi) + 20.0 * Math.Sin(2.0 * x * pi)) * 2.0 / 3.0; ret += (20.0 * Math.Sin(x * pi) + 40.0 * Math.Sin(x / 3.0 * pi)) * 2.0 / 3.0; ret += (150.0 * Math.Sin(x / 12.0 * pi) + 300.0 * Math.Sin(x / 30.0 * pi)) * 2.0 / 3.0; return ret; } #endregion #region 百度坐标转火星坐标 const double x_pi = 3.14159265358979324 * 3000.0 / 180.0; static MapPoint LonLatToBaidu(MapPoint mapPointLon) { MapPoint rst = new MapPoint(); double bd_lat = 0; double bd_lon = 0; LonLatToBaidu(mapPointLon.X, mapPointLon.Y, ref bd_lat, ref bd_lon); rst.X = bd_lon; rst.Y = bd_lat; return rst; } static void LonLatToBaidu(double lon, double lat, ref double bd_lat, ref double bd_lon) { double x = lon, y = lat; double z = Math.Sqrt(x * x + y * y) + 0.00002 * Math.Sin(y * x_pi); double theta = Math.Atan2(y, x) + 0.000003 * Math.Cos(x * x_pi); bd_lon = z * Math.Cos(theta) + 0.0065; bd_lat = z * Math.Sin(theta) + 0.006; } static void BaiduToLonLat(double bd_lat, double bd_lon, ref double lon, ref double lat) { double x = bd_lon - 0.0065, y = bd_lat - 0.006; double z = Math.Sqrt(x * x + y * y) - 0.00002 * Math.Sin(y * x_pi); double theta = Math.Atan2(y, x) - 0.000003 * Math.Cos(x * x_pi); lon = z * Math.Cos(theta); lat = z * Math.Sin(theta); } #endregion } /// /// 四参数实体 /// public class ParamsEntity { /// /// 标识码 /// public int BSM { get; set; } /// /// 行政区代码 /// public string XZQDM { get; set; } /// /// X平移84转80 /// public double X84OFFSET { get; set; } /// /// Y平移84转80 /// public double Y84OFFSET { get; set; } /// /// X旋转84转80 /// public double X84ROTATE { get; set; } /// /// Y旋转84转80 /// public double Y84ROTATE { get; set; } /// /// X平移80转84 /// public double X80OFFSET { get; set; } /// /// Y平移80转84 /// public double Y80OFFSET { get; set; } /// /// X旋转80转84 /// public double X80ROTATE { get; set; } /// /// Y旋转80转84 /// public double Y80ROTATE { get; set; } /// ///代号 /// public string DH { get; set; } } public class ObjectConvertOperate { /// /// Object转string /// /// /// public static string ObjectToString(object obj) { if (obj == null) return ""; return obj.ToString(); } /// /// Object转double /// /// /// public static double ObjectToDouble(object obj) { if (obj == null) return 0.0; double rst = 0.0; double.TryParse(obj.ToString(), out rst); return rst; } /// /// Object转double /// /// /// public static int ObjectToInt(object obj) { if (obj == null) return 0; int rst = 0; int.TryParse(obj.ToString(), out rst); return rst; } /// /// Object转bool /// /// /// public static bool ObjectToBool(object obj) { if (obj == null) return false; bool rst = false; bool.TryParse(obj.ToString(), out rst); return rst; } /// /// Object转datetime /// /// /// public static DateTime ObjectToDate(object obj) { if (obj == null) return new DateTime(); DateTime dt = new DateTime(); DateTime.TryParse(obj.ToString(), out dt); return dt; } } public class JsonRings { public List> rings; public static JsonRings GetTemp() { JsonRings jr = new JsonRings() { rings = new List>() }; for (int i = 0; i < 10; i++) { List ring = new List(); for (int j = 0; j < 10; j++) { ring.Add(new double[] { 1111, 1111 }); } jr.rings.Add(ring); } return jr; } public string GetJsonString() { var serializer = new DataContractJsonSerializer(this.GetType()); var stream = new MemoryStream(); serializer.WriteObject(stream, this); byte[] dataBytes = new byte[stream.Length]; stream.Position = 0; stream.Read(dataBytes, 0, (int)stream.Length); return Encoding.UTF8.GetString(dataBytes); } public static JsonRings FromJson(string json) { try { if (string.IsNullOrWhiteSpace(json)) return null; var serializer = new DataContractJsonSerializer(typeof(JsonRings)); var mStream = new MemoryStream(Encoding.Default.GetBytes(json)); return (JsonRings)serializer.ReadObject(mStream); } catch { return null; } } } }