程序定義一個(gè)投影的Transform的類，橢球ellipsoid為傳入的參數(shù)，橢球相關(guān)的內(nèi)容可見這篇博客

高斯投影正算是傳入大地坐標(biāo)與中央經(jīng)線經(jīng)度，計(jì)算得到該投影帶獨(dú)立坐標(biāo)系下的坐標(biāo)。若有要求可對(duì)y坐標(biāo)東移，+500000，加以帶號(hào)表示得到Y(jié)的通用坐標(biāo)。

本程序高斯投影反算是在獨(dú)立坐標(biāo)系的前提下，輸入點(diǎn)的xy坐標(biāo)和中央經(jīng)線經(jīng)度，迭代計(jì)算得到大地經(jīng)緯度。
計(jì)算時(shí)要分清是通用坐標(biāo)還是獨(dú)立坐標(biāo)系下的坐標(biāo)，不然會(huì)計(jì)算錯(cuò)誤。
計(jì)算公式是由高斯投影正反算公式經(jīng)過推導(dǎo)得到的適于電算的公式。

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;//輔助計(jì)算公式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);//子午線弧長(zhǎng)計(jì)算公式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;//子午弧長(zhǎng)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);//經(jīng)差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"><"高斯平面直角坐標(biāo)系坐標(biāo)">/// <param name="L0"><"中央子午線經(jīng)度">/// <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;}}

測(cè)繪學(xué)科中的經(jīng)緯度坐標(biāo)一般為dd.mmssssss格式，如114.123345表示114°，12′，33.45″，用以下AngleRadian類進(jìn)行弧度角度轉(zhuǎn)換。

class AngleRadian{/// <summary>/// dd.mmssssss格式的角度轉(zhuǎn)換為弧度/// </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>/// 弧度轉(zhuǎn)換為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);//返回四位小數(shù)double degrees = d + m / 100 + s / 10000;return degrees;}/// <summary>/// dd.mmssssss格式的角度轉(zhuǎn)換為°′″形式的字符串/// </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>/// 弧度轉(zhuǎn)換為°′″形式的字符串/// </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);//返回四位小數(shù)string dms = symbol + d.ToString() + "°" + m.ToString() + "′" + s.ToString() + "″";return dms;}}

淺談高斯投影

1、分帶

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

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

為了指明該點(diǎn)屬于何帶，還規(guī)定在橫坐標(biāo)y值之前，要寫上帶號(hào)。
因此坐標(biāo)值表現(xiàn)形式有三種：自然值、+500KM值、通用值。

所以拿到一個(gè)坐標(biāo)應(yīng)當(dāng)進(jìn)行判斷它到底是哪一種類型的坐標(biāo)值，本程序的轉(zhuǎn)換是基于自然值的轉(zhuǎn)換

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。