程序定義一個投影的Transform的類，橢球ellipsoid為傳入的參數，橢球相關的內容可見這篇博客

高斯投影正算是傳入大地坐標與中央經線經度，計算得到該投影帶獨立坐標系下的坐標。若有要求可對y坐標東移，+500000，加以帶號表示得到Y的通用坐標。

本程序高斯投影反算是在獨立坐標系的前提下，輸入點的xy坐標和中央經線經度，迭代計算得到大地經緯度。
計算時要分清是通用坐標還是獨立坐標系下的坐標，不然會計算錯誤。
計算公式是由高斯投影正反算公式經過推導得到的適于電算的公式。

class Transform{public double a;public double ec;public double ecc;public Transform(Ellipsoid ellipsoid){a = ellipsoid.a;ec = ellipsoid.ec;ecc = ellipsoid.ecc;}//高斯投影正算public Pointxy BLToxy(PointBL bL,double L0){double B = bL.B;double L = bL.L;//輔助計算公式double W = Math.Sqrt(1 - ec * Math.Sin(B) * Math.Sin(B));double n2 = ecc * Math.Cos(B) * Math.Cos(B);double t = Math.Tan(B);//曲率半徑double N = a / W;double M = a * (1 - ec) / Math.Pow(W, 3);double M0 = a * (1 - ec);//子午線弧長計算公式double Ac = 1 + 3 / 4d * ec + 45 / 64d * Math.Pow(ec, 2) + 175 / 256d * Math.Pow(ec, 3) + 11025 / 16384d * Math.Pow(ec, 4) + 43659 / 65536d * Math.Pow(ec, 5);double Bc = 3 / 4d * ec + 15 / 16d * Math.Pow(ec, 2) + 525 / 512d * Math.Pow(ec, 3) + 2205 / 2048d * Math.Pow(ec, 4) + 72765 / 65536d * Math.Pow(ec, 5);double Cc = 15 / 64d * Math.Pow(ec, 2) + 105 / 256d * Math.Pow(ec, 3) + 2205 / 4096d * Math.Pow(ec, 4) + 10395 / 16384d * Math.Pow(ec, 5);double Dc = 35 / 512d * Math.Pow(ec, 3) + 315 / 2048d * Math.Pow(ec, 4) + 31185 / 131072d * Math.Pow(ec, 5);double Ec = 315 / 16384d * Math.Pow(ec, 4) + 3465 / 65536d * Math.Pow(ec, 5);double Fc = 693 / 131072d * Math.Pow(ec, 5);double Alpha = Ac * M0;double Beta = -1 / 2d * Bc * M0;double Gamma = 1 / 4d * Cc * M0;double Delte = -1 / 6d * Dc * M0;double Epsilon = 1 / 8d * Ec * M0;double Zeta = -1 / 10d * Fc * M0;//子午弧長double X = Alpha * B + Beta * Math.Sin(2 * B) + Gamma * Math.Sin(4 * B) + Delte * Math.Sin(6 * B) + Epsilon * Math.Sin(8 * B) + Zeta * Math.Sin(10 * B);//經差double l = L - L0 / 180 * Math.PI;//輔助量double a0 = X;double a1 = N * Math.Cos(B);double a2 = 1 / 2d * N * Math.Pow(Math.Cos(B), 2) * t;double a3 = 1 / 6d * N * Math.Pow(Math.Cos(B), 3) * (1 - Math.Pow(t, 2) + n2);double a4 = 1 / 24d * N * Math.Pow(Math.Cos(B), 4) * (5 - Math.Pow(t, 2) + 9 * n2 + 4 * Math.Pow(n2, 2));double a5 = 1 / 120d * N * Math.Pow(Math.Cos(B), 5) * (5 - 18 * Math.Pow(t, 2) + Math.Pow(t, 4) + 14 * n2 - 58 * n2 * Math.Pow(t, 2));double a6 = 1 / 720d * N * Math.Pow(Math.Cos(B), 6) * (61 - 58 * Math.Pow(t, 2) + Math.Pow(t, 4) + 270 * n2 - 330 * n2 * Math.Pow(t, 2)) * t;Pointxy xy = new Pointxy();xy.x = a0 + a2 * Math.Pow(l, 2) + a4 * Math.Pow(l, 4) + a6 * Math.Pow(l, 6);xy.y = a1 * l + a3 * Math.Pow(l, 3) + a5 * Math.Pow(l, 5);return xy;}/// <summary>/// 高斯投影反算/// </summary>/// <param name="xy"><"高斯平面直角坐標系坐標">/// <param name="L0"><"中央子午線經度">/// <returns></returns>public PointBL xyToBL(Pointxy xy,double L0){double x = xy.x;double y = xy.y;double Ac = 1 + 3 / 4d * ec + 45 / 64d * Math.Pow(ec, 2) + 175 / 256d * Math.Pow(ec, 3) + 11025 / 16384d * Math.Pow(ec, 4) + 43659 / 65536d * Math.Pow(ec, 5);double Bc = 3 / 4d * ec + 15 / 16d * Math.Pow(ec, 2) + 525 / 512d * Math.Pow(ec, 3) + 2205 / 2048d * Math.Pow(ec, 4) + 72765 / 65536d * Math.Pow(ec, 5);double Cc = 15 / 64d * Math.Pow(ec, 2) + 105 / 256d * Math.Pow(ec, 3) + 2205 / 4096d * Math.Pow(ec, 4) + 10395 / 16384d * Math.Pow(ec, 5);double Dc = 35 / 512d * Math.Pow(ec, 3) + 315 / 2048d * Math.Pow(ec, 4) + 31185 / 131072d * Math.Pow(ec, 5);double Ec = 315 / 16384d * Math.Pow(ec, 4) + 3465 / 65536d * Math.Pow(ec, 5);double Fc = 693 / 131072d * Math.Pow(ec, 5);double M0 = a * (1 - ec);double Alpha = Ac * M0;double Beta = -1 / 2d * Bc * M0;double Gamma = 1 / 4d * Cc * M0;double Delte = -1 / 6d * Dc * M0;double Epsilon = 1 / 8d * Ec * M0;double Zeta = -1 / 10d * Fc * M0;double X = x;double B0 = X / Alpha;double Bf = 0;while (true){double dert = Beta * Math.Sin(2 * B0) + Gamma * Math.Sin(4 * B0) + Delte * Math.Sin(6 * B0) + Epsilon * Math.Sin(8 * B0) + Zeta * Math.Sin(10 * B0);Bf = (X - dert) / Alpha;if (Math.Abs(Bf - B0) < 0.00000001)break;elseB0 = Bf;}//輔助公式double Wf = Math.Sqrt(1 - ec * Math.Pow(Math.Sin(Bf), 2));double n2 = ecc * Math.Pow(Math.Cos(Bf), 2);double tf = Math.Tan(Bf);double Nf = a / Wf;double Mf = a * (1 - ec) / Math.Pow(Wf, 3);double b0 = Bf;double b1 = 1 / (Nf * Math.Cos(Bf));double b2 = -tf / (2 * Nf * Mf);double b3 = -(1 + 2 * tf * tf + n2) * b1 / (6 * Nf * Nf);double b4 = -(5 + 3 * tf * tf + n2 - 9 * n2 * tf * tf) * b2 / (12 * Nf * Nf);double b5 = -(5 + 28 * tf * tf + 24 * Math.Pow(tf, 4) + 6 * n2 + 8 * n2 * tf * tf) * b1 / (120 * Math.Pow(Nf, 4));double b6 = (61 + 90 * tf * tf + 45 * Math.Pow(tf, 4)) * b2 / (360 * Math.Pow(Nf, 4));PointBL BL = new PointBL();BL.B = b0 + b2 * Math.Pow(y, 2) + b4 * Math.Pow(y, 4) + b6 * Math.Pow(y, 6);BL.L = b1 * Math.Pow(y, 1) + b3 * Math.Pow(y, 3) + b5 * Math.Pow(y, 5) + L0 * Math.PI / 180;return BL;}}

測繪學科中的經緯度坐標一般為dd.mmssssss格式，如114.123345表示114°，12′，33.45″，用以下AngleRadian類進行弧度角度轉換。

class AngleRadian{/// <summary>/// dd.mmssssss格式的角度轉換為弧度/// </summary>/// <param name="degrees"></param>/// <returns></returns>public double ConvertDegreesToRadians(double degrees){double d = Math.Truncate(degrees);double m = Math.Truncate((degrees - d) * 100);double s = ((degrees - d) * 100 - m) * 100;double radians = (d + m / 60 + s / 3600) / 180 * Math.PI;return radians;}/// <summary>/// 弧度轉換為dd.mmssssss格式的角度/// </summary>/// <param name="radians"></param>/// <returns></returns>public double ConvertRadiansToDegrees(double radians){double dd = radians / Math.PI * 180;double d = Math.Truncate(dd);double m = Math.Truncate((dd - d) * 60);double s = Math.Round(((dd - d) * 60 - m) * 60, 4);//返回四位小數double degrees = d + m / 100 + s / 10000;return degrees;}/// <summary>/// dd.mmssssss格式的角度轉換為°′″形式的字符串/// </summary>/// <param name="degrees"></param>/// <returns></returns>public string ConvertDegreesToString(double degrees){string symbol = "";if (degrees < 0){degrees = Math.Abs(degrees);symbol = "-";}double d = Math.Truncate(degrees);double m = Math.Truncate((degrees - d) * 100);double s = ((degrees - d) * 100 - m) * 100;string dms = symbol + d.ToString() + "°" + m.ToString() + "′" + s.ToString() + "″";return dms;}/// <summary>/// 弧度轉換為°′″形式的字符串/// </summary>/// <param name="radians"></param>/// <returns></returns>public string ConvertRadiansToString(double radians){string symbol = "";if (radians < 0){radians = Math.Abs(radians);symbol = "-";}double dd = radians / Math.PI * 180;double d = Math.Truncate(dd);double m = Math.Truncate((dd - d) * 60);double s = Math.Round(((dd - d) * 60 - m) * 60, 4);//返回四位小數string dms = symbol + d.ToString() + "°" + m.ToString() + "′" + s.ToString() + "″";return dms;}}

淺談高斯投影

1、分帶

我國采用6、3度帶。中央子午線與帶號的關系
6度帶
自格林威治零度經線起，每6度分為一個投影帶，自西向東分帶，全球共分為60個投影帶，帶號依次編為第 1、2…60帶。我國6°帶中央子午線的經度，由73°起每隔6°而至135°，共計11帶，帶號用n表示，中央子午線的經度用Ｌ０表示。
L=6n-3（n為帶號，L為中央經線經度）
3度帶
是在6度帶的基礎上分成的，它的中央子午線與六度帶的中央子午線和分帶子午線重合，即自 1.5度子午線起每隔經差3度自西向東分帶，帶號依次編為三度帶第 1、2…120帶。中央子午線經度依次為3°, 6°, 9°, … , 360°。
L=3n（n為帶號，L為中央經線經度）

2、在我國x坐標都是正的，y坐標的最大值(在赤道上)約為330km。為了避免出現負的y坐標，則無論3°或6°帶，每帶的縱坐標軸要西移500 km，即在每帶的橫（y）坐標上加500 km。

為了指明該點屬于何帶，還規定在橫坐標y值之前，要寫上帶號。
因此坐標值表現形式有三種：自然值、+500KM值、通用值。

所以拿到一個坐標應當進行判斷它到底是哪一種類型的坐標值，本程序的轉換是基于自然值的轉換

總結

以上是生活随笔為你收集整理的地图投影（一）高斯克吕格投影的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。