在前面的文章中系統的闡述了工程坐標的轉換類別和轉換的方法。關于轉換代碼實現,有很多的類庫:
這里針對GPS接收的WGS84橢球的經緯度轉換為地方坐標系的問題,利用C#,對工程坐標轉換方法和步驟做出詳細的解答。不基于任何類庫和函數庫,也未使用矩陣庫,可以便利的將代碼移植到任何語言。
根據上一篇文章中對七參數、四參數、高程擬合在坐標轉換的作用和使用條件的闡述,我們可以將上一篇文章第7節的總結圖,按照計算的流程重新繪制。
根據上圖可知,預將WGS84橢球的GPS坐標需要經過5次轉換。其中,
因此,根據計算原理,直接可以利用C#代碼實現。
5個轉換是對點的操作,不妨構建自定義點類MyPoint
,在這個類中定義轉換方法。在實現轉換方法之前,需要定義數據屬性,以承載轉換參數和轉換數據。代碼框架如下:
internal class MyPoint
{// 定義橢球類型。這里僅列舉了4中國內常見的橢球類型
// 國際橢球可以增加自行定義
public enum EllipsoidType
{ WGS84,
CGCS2000,
西安80,
北京54
}
//大地坐標經度、維度、高度
public double L {get; set; }
public double B {get; set; }
public double H {get; set; }
//空間坐標系
public double X {get; set; }
public double Y {get; set; }
public double Z {get; set; }
//七參數轉換后的空間坐標
public double X2 {get; set; }
public double Y2 {get; set; }
public double Z2 {get; set; }
private double a = 0, f = 0, b = 0, e = 0, e2 = 0; //橢球參數
private readonly double rho = 180 / Math.PI;
private readonly double d2r = Math.PI / 180;
public double Xs {get; set; }
public double Ys {get; set; }
public double Hs {get; set; }
//七參數 三個線性平移量-單位米 三個旋轉平移量-十進制秒為單位(運算時注意轉換為度) 比例因子-單位百萬分率 (ppm)
//測量隊給出的七參數單位與計算的單位不同,要進行單位轉化 1 秒=0.0000048481373323 弧度
//尺度因子有兩種單位的表示形式,一種結果約為1,如1.0000045,用k表示;
//另一種就是ppm的表示形式,稍微比1大一點,如4.5,用m表示。k=m/1000000
private double dx = 0, dy = 0, dz = 0, rx = 0, ry = 0, rz = 0, m = 0, k = 0;
}
3.2 橢球參數賦值常見的橢球參數值在我的文章經緯度坐標轉換為工程坐標可以找到,這里選取與上述代碼對應的4類橢球,并在上述MyPoint
類中增加函數EllipsoidParam(EllipsoidType type)
。
////// 橢球參數設置
/// ///橢球類型private void EllipsoidParam(EllipsoidType type)
{// CGCS2000 橢球參數
if (type == EllipsoidType.CGCS2000)
{this.a = 6378137;
this.f = 1 / 298.257222101;
}
// 西安 80
else if (type == EllipsoidType.西安80)
{this.a = 6378140;
this.f = 1 / 298.257;
}
// 北京 54
else if (type == EllipsoidType.北京54)
{this.a = 6378245;
this.f = 1 / 298.3;
}
// WGS-84
else
{this.a = 6378137;
this.f = 1 / 298.257223563;
}
this.b = this.a * (1 - this.f);
this.e = Math.Sqrt(this.a * this.a - this.b * this.b) / this.a; //第一偏心率
this.e2 = Math.Sqrt(this.a * this.a - this.b * this.b) / this.b; //第二偏心率
}
3.3 轉換1、3(大地經緯度坐標與地心地固坐標的轉換)charlee44的博客有C++代碼的實現,現在利用C#重構即可。上述MyPoint
類中增加BLH2XYZ(EllipsoidType type)
和XYZ2BLH(EllipsoidType type)
兩個函數。
////// 經緯度坐標轉空間直角坐標
/// ///橢球類型public void BLH2XYZ(EllipsoidType type = EllipsoidType.WGS84)
{EllipsoidParam(type);
double sB = Math.Sin(this.B * d2r);
double cB = Math.Cos(this.B * d2r);
double sL = Math.Sin(this.L * d2r);
double cL = Math.Cos(this.L * d2r);
double N = this.a / (Math.Sqrt(1 - this.e * this.e * sB * sB));
this.X = (N + this.H) * cB * cL;
this.Y = (N + this.H) * cB * sL;
this.Z = (N * (1 - this.e * this.e) + this.H) * sB;
this.X2 = this.X;
this.Y2 = this.Y;
this.Z2 = this.Z;
}
////// 空間直角坐標轉經緯度坐標
/// ///橢球類型public void XYZ2BLH(EllipsoidType type)
{EllipsoidParam(type);
// 這里轉出來的B L是弧度
this.L = Math.Atan(this.Y2 / this.X2) + Math.PI;
this.L = this.L * 180 / Math.PI;
// B需要迭代計算
double B2 = Math.Atan(Z2 / Math.Sqrt(X2 * X2 + Y2 * Y2));
double B1;
double N;
while (true)
{N = a / Math.Sqrt(1 - f * (2 - f) * Math.Sin(B2) * Math.Sin(B2));
B1 = Math.Atan((Z2 + N * f * (2 - f) * Math.Sin(B2)) / Math.Sqrt(X2 * X2 + Y2 * Y2));
if (Math.Abs(B1 - B2)< 1e-12)
break;
B2 = B1;
}
this.B = B2 * 180 / Math.PI;
double sB = Math.Sin(this.B * d2r);
double cB = Math.Cos(this.B * d2r);
this.H = this.Z2 / sB - N * (1 - this.e * this.e);
}
3.4 投影轉換此處僅實現了常見的高斯-克里格投影。上述MyPoint
類中增加GaussProjection(EllipsoidType type, ProjectionSetting prjSetting)
函數。
////// 利用高斯投影將指定橢球類型的經緯度坐標轉為投影坐標
/// ///橢球類型///投影設置實例public void GaussProjection(EllipsoidType type, ProjectionSetting prjSetting)
{this.EllipsoidParam(type);
double l = (this.L - prjSetting.CenterL) / this.rho;
double cB = Math.Cos(this.B * this.d2r);
double sB = Math.Sin(this.B * this.d2r);
double s2b = Math.Sin(this.B * this.d2r * 2);
double s4b = Math.Sin(this.B * this.d2r * 4);
double s6b = Math.Sin(this.B * this.d2r * 6);
double s8b = Math.Sin(this.B * this.d2r * 8);
double N = this.a / Math.Sqrt(1 - this.e * this.e * sB * sB); // 卯酉圈曲率半徑
double t = Math.Tan(this.B * this.d2r);
double eta = this.e2 * cB;
double m0 = this.a * (1 - this.e * this.e);
double m2 = 3.0 / 2.0 * this.e * this.e * m0;
double m4 = 5.0 / 4.0 * this.e * this.e * m2;
double m6 = 7.0 / 6.0 * this.e * this.e * m4;
double m8 = 9.0 / 8.0 * this.e * this.e * m6;
double a0 = m0 + 1.0 / 2.0 * m2 + 3.0 / 8.0 * m4 + 5.0 / 16.0 * m6 + 35.0 / 128.0 * m8;
double a2 = 1.0 / 2.0 * m2 + 1.0 / 2.0 * m4 + 15.0 / 32.0 * m6 + 7.0 / 16.0 * m8;
double a4 = 1.0 / 8.0 * m4 + 3.0 / 16.0 * m6 + 7.0 / 32.0 * m8;
double a6 = 1.0 / 32.0 * m6 + 1.0 / 16.0 * m8;
double a8 = 1.0 / 128.0 * m8;
// X1為自赤道量起的子午線弧長
double X1 = a0 * (this.B * this.d2r) - 1.0 / 2.0 * a2 * s2b + 1.0 / 4.0 * a4 * s4b - 1.0 / 6.0 * a6 * s6b + 1.0 / 8.0 * a8 * s8b;
this.Xs = X1 + N / 2 * t * cB * cB * l * l + N / 24 * t * (5 - t * t + 9 * Math.Pow(eta, 2) + 4 * Math.Pow(eta, 4)) * Math.Pow(cB, 4) * Math.Pow(l, 4)
+ N / 720 * t * (61 - 58 * t * t + Math.Pow(t, 4)) * Math.Pow(cB, 6) * Math.Pow(l, 6);
this.Ys = N * cB * l + N / 6 * (1 - t * t + eta * eta) * Math.Pow(cB, 3) * Math.Pow(l, 3)
+ N / 120 * (5 - 18 * t * t + Math.Pow(t, 4) + 14 * Math.Pow(eta, 2) - 58 * eta * eta * t * t) * Math.Pow(cB, 5) * Math.Pow(l, 5);
this.Hs = this.H;
// 假東 假北偏移
this.Xs += prjSetting.PseudoNorth;
this.Ys += prjSetting.PseudoEast;
}
其中,ProjectionSetting
是一個投影參數設置類,獨立于MyPoint
類,用于設定中央經線、東偏等投影參數。
internal class ProjectionSetting
{private double _centerL;
public double CenterL
{ get {return _centerL; }
set {_centerL = value; }
}
private double _centerB;
public double CenterB
{ get {return _centerB; }
set {_centerB = value; }
}
private double _pseudoEast;
public double PseudoEast
{ get {return _pseudoEast; }
set {_pseudoEast = value; }
}
private double _pseudoNorth;
public double PseudoNorth
{ get {return _pseudoNorth; }
set {_pseudoNorth = value; }
}
private double _prjScale;
public double PrjScale
{ get {return _prjScale; }
set {_prjScale = value; }
}
////// 設置全部的投影參數
/// ///////////////public ProjectionSetting(double centerL, double centerB,
double pseudoEast, double pseudoNorth,
double prjScale)
{ CenterL = centerL;
CenterB = centerB;
PseudoEast = pseudoEast;
PseudoNorth = pseudoNorth;
PrjScale = prjScale;
}
////// 僅設置中央經線和東偏
/// //////public ProjectionSetting(double centerL, double pseudoEast)
{CenterL = centerL;
CenterB = 0.0;
PseudoEast = pseudoEast;
PseudoNorth = 0.0;
PrjScale = 1.0;
}
////// 默認常用投影參數,中央經線120°,東偏500000
/// public ProjectionSetting()
{CenterL = 120.0;
CenterB = 0.0;
PseudoEast = 500000;
PseudoNorth = 0.0;
PrjScale = 1.0;
}
}
3.5 轉換2的實現(三參數、七參數)上述MyPoint
類中增加SevenParamTrans(Datum7Paras datum7Paras)
和TreeParamTrans(Datum3Paras datum3Paras)
函數。
////// 利用7參數進行坐標系之間轉換
/// ///7參數實例public void SevenParamTrans(Datum7Paras datum7Paras)
{this.dx = datum7Paras.Dx;
this.dy = datum7Paras.Dy;
this.dz = datum7Paras.Dz;
this.rx = datum7Paras.Rx * 0.0000048481373323; //1 秒=0.0000048481373323 弧度
this.ry = datum7Paras.Ry * 0.0000048481373323;
this.rz = datum7Paras.Rz * 0.0000048481373323;
this.m = datum7Paras.PPM;
this.k = this.m / 1000000;
this.X2 = (1 + k) * (this.X + this.rz * this.Y - this.ry * this.Z) + this.dx;
this.Y2 = (1 + k) * (-this.rz * this.X + this.Y + this.rx * this.Z) + this.dy;
this.Z2 = (1 + k) * (this.ry * this.X - this.rx * this.Y + this.Z) + this.dz;
}
////// 利用3參數進行坐標系之間轉換
/// ///3參數實例public void TreeParamTrans(Datum3Paras datum3Paras)
{this.dx = datum3Paras.Dx;
this.dy = datum3Paras.Dy;
this.dz = datum3Paras.Dz;
this.X2 = this.X + this.dx;
this.Y2 = this.Y + this.dy;
this.Z2 = this.Z + this.dz;
}
Datum3Paras
和Datum7Paras
是獨立于MyPoint
類,用于設定坐標轉換參數。
////// 7參數
/// internal class Datum7Paras
{private double _dx;
public double Dx
{ get {return _dx; }
set {_dx = value; }
}
private double _dy;
public double Dy
{get {return _dy; }
set {_dy = value; }
}
private double _dz;
public double Dz
{get {return _dz; }
set {_dz = value; }
}
private double _rx;
public double Rx
{get {return _rx; }
set {_rx = value; }
}
private double _ry;
public double Ry
{get {return _ry; }
set {_ry = value; }
}
private double _rz;
public double Rz
{get {return _rz; }
set {_rz = value; }
}
private double _ppm;
public double PPM
{get {return _ppm; }
set {_ppm = value; }
}
public Datum7Paras(double dx, double dy, double dz,
double rx, double ry, double rz,
double ppm)
{_dx = dx;
_dy = dy;
_dz = dz;
_rx = rx;
_ry = ry;
_rz = rz;
_ppm = ppm;
}
}
internal class Datum3Paras
{private double _dx;
public double Dx
{get {return _dx; }
set {_dx = value; }
}
private double _dy;
public double Dy
{get {return _dy; }
set {_dy = value; }
}
private double _dz;
public double Dz
{get {return _dz; }
set {_dz = value; }
}
public Datum3Paras(double dx, double dy, double dz)
{Dx = dx;
Dy = dy;
Dz = dz;
}
}
3.6 轉換5的實現(四參數+高程擬合)上述MyPoint
類中增加Transform4Para(Trans4Paras transPara)
函數。此處,高程擬合僅實現了已知一個測點的固定改正差。
////// 投影坐標獲取后,進一步利用4參數轉換坐標
/// ///public void Transform4Para(Trans4Paras transPara)
{var X1 = transPara.Dx;
var Y1 = transPara.Dy;
var cosAngle = Math.Cos(transPara.A);
var sinAngle = Math.Sin(transPara.A);
X1 += transPara.K * (cosAngle * this.Xs - sinAngle * this.Ys);
Y1 += transPara.K * (sinAngle * this.Xs + cosAngle * this.Ys);
this.Xs = X1;
this.Ys = Y1;
// 固定改正差
this.Hs += transPara.Dh;
}
Trans4Paras
是獨立于MyPoint
類,用于設定坐標轉換參數。
internal class Trans4Paras
{private double _dx;
public double Dx
{get {return _dx; }
set {_dx = value; }
}
private double _dy;
public double Dy
{get {return _dy; }
set {_dy = value; }
}
private double _a;
public double A
{get {return _a; }
set {_a = value; }
}
private double _k;
public double K
{get {return _k; }
set {_k = value; }
}
private double _dh;
public double Dh
{get {return _dh; }
set {_dh = value; }
}
public Trans4Paras(double dx, double dy, double a, double k, double dh)
{Dx = dx;
Dy = dy;
A = a;
K = k;
Dh = dh;
}
public Trans4Paras()
{}
}
3.7 調用過程里面的參數,因為保密原因,做出了隨機更改,實際使用時可根據自己情況賦值。
3.7.1 一步法// 實例化計算參數
MyPoint p = new MyPoint();.
p.L=113.256;
p.B=31.565;
p.H=5.216;
// 經緯度轉空間坐標
p.BLH2XYZ();
// 實例化七參數
Datum7Paras datum7Paras = new Datum7Paras(
489.2994563566, 141.1525159753, 15.74421120568,
-0.164423, 4.141573, -4.808299,
-6.56482989958);
p.SevenParamTrans(datum7Paras);
// 空間坐標轉回經緯度
p.XYZ2BLH(EllipsoidType.WGS84);
// 高斯投影 經緯度轉平面坐標
// 實例化投影參數類
ProjectionSetting projectionSetting = new ProjectionSetting(120,500000);
p.GaussProjection(EllipsoidType.WGS84, projectionSetting);
3.7.2 兩步法// 實例化計算參數
MyPoint p = new MyPoint();.
p.SetLBH(113.256,31.565,5.216);
// 經緯度轉空間坐標
p.BLH2XYZ();
// 實例化七參數
Datum7Paras datum7Paras = new Datum7Paras(
489.2994563566, 141.1525159753, 15.74421120568,
-0.164423, 4.141573, -4.808299,
-6.56482989958);
p.SevenParamTrans(datum7Paras);
// 空間坐標轉回經緯度
p.XYZ2BLH(EllipsoidType.WGS84);
// 高斯投影 經緯度轉平面坐標
// 實例化投影參數類
ProjectionSetting projectionSetting = new ProjectionSetting(120,500000);
p.GaussProjection(EllipsoidType.WGS84, projectionSetting);
Trans4Paras transformPara = new(6456.15957352521, -134618.390707439, 0.011104964500129, 1.00002537583871, 5.788);
p.Transform4Para(transformPara);
至此,關于工程坐標系轉化,即GPS接收的WGS84橢球的經緯度轉換為地方坐標系的問題,基本全部實現。代碼正確性和準確性的驗證是與 南方GPS工具箱做對比。例如,采用上述的一步法,在設定好坐標、7參數、投影參數后,計算發現,與南方GPS工具箱在y方向偏差1mm。結果如下圖:
你是否還在尋找穩定的海外服務器提供商?創新互聯www.cdcxhl.cn海外機房具備T級流量清洗系統配攻擊溯源,準確流量調度確保服務器高可用性,企業級服務器適合批量采購,新人活動首月15元起,快前往官網查看詳情吧
分享名稱:工程坐標轉換方法C#代碼實現-創新互聯
分享URL:http://vcdvsql.cn/article18/cescgp.html
成都網站建設公司_創新互聯,為您提供域名注冊、網站導航、做網站、響應式網站、自適應網站、網站維護
聲明:本網站發布的內容(圖片、視頻和文字)以用戶投稿、用戶轉載內容為主,如果涉及侵權請盡快告知,我們將會在第一時間刪除。文章觀點不代表本網站立場,如需處理請聯系客服。電話:028-86922220;郵箱:631063699@qq.com。內容未經允許不得轉載,或轉載時需注明來源: 創新互聯