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;
}
}
}
}