史上最全的高斯投影源代码集合精心收集9种.docx

1、史上最全的高斯投影源代码集合精心收集9种高斯投影C+源码 1CCoordinate CEarthPlaneDialog:EarthToPlane(const CEarth & earth) double B = earth.m_dB.m_realRadians; double L = earth.m_dL.m_realRadians; L = SetRad_PIPI(L); double L0 = SetL0(L); double X, N, t, t2, m, m2, ng2; double sinB, cosB; X = A1 * B * 180.0 / PI + A2 * sin(2 *

2、 B) + A3 * sin(4 * B) + A4 * sin(6 * B); sinB = sin(B); cosB = cos(B); t = tan(B); t2 = t * t; N = a / sqrt(1 - e2 * sinB * sinB); m = cosB * (L - L0); m2 = m * m; ng2 = cosB * cosB * e2 / (1 - e2); double x = X + N * t * (0.5 + (5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) *

3、m2 / 720.0) * m2) * m2); double y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0); y += 500000; CCoordinate coor(x,y); return coor;/将大地坐标系转换为平面直角坐标系CEarth CEarthPlaneDialog:PlaneToEarth(const CCoordinate & coor, UINT Sign) double x = c

4、oor.m_dX; double y = coor.m_dY; double L0 = (6.0 * Sign - 3.0) / 180 * PI; double sinB, cosB, t, t2, N ,ng2, V, yN; double preB0, B0; double eta; y -= 500000; B0 = x / A1; do preB0 = B0; B0 = B0 * PI / 180.0; B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0) / A1; eta = fabs(B0 - pre

5、B0); while(eta 0.0000000000000001); B0 = B0 * PI / 180.0; sinB = sin(B0); cosB = cos(B0); t = tan(B0); t2 = t * t; N = a / sqrt(1 - e2 * sinB * sinB); ng2 = cosB * cosB * e2 / (1 - e2); V = sqrt(1 + ng2); yN = y / N; double B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN /

6、12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2; double L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB; CEarth earth(B,L,0,r); return earth;/将平面直角坐标系转换

7、为大地坐标系void CEarthPlaneDialog:OnRADIOPlane1954() / TODO: Add your control notification handler code here a = 6378245; f = 298.3; e2 = 1 - (f - 1) / f) * (f - 1) / f); e12 = (f / (f - 1) * (f / (f - 1) - 1; A1 = 111134.8611; A2 = -16036.4803; A3 = 16.8281; A4 = -0.0220; m_iReferenceFrame = 0; void CEa

8、rthPlaneDialog:OnRADIOPlane1980() / TODO: Add your control notification handler code here a = 6378140; f = 298.257; e2 = 1 - (f - 1) / f) * (f - 1) / f); e12 = (f / (f - 1) * (f / (f - 1) - 1; A1 = 111133.0047; A2 = -16038.5282; A3 = 16.8326; A4 = -0.0220; m_iReferenceFrame = 1; void CEarthPlaneDial

9、og:OnRADIOPlaneWGS84() / TODO: Add your control notification handler code here a = 6378137; f = 298.257223563; e2 = 1 - (f - 1) / f) * (f - 1) / f); e12 = (f / (f - 1) * (f / (f - 1) - 1; A1 = 111133.0047; A2 = -16038.5282; A3 = 16.8326; A4 = -0.0220; m_iReferenceFrame = 2; double CEarthPlaneDialog:

10、SetL0(double L) L = L / PI * 180; if( L 0 ) for(m_iCincture = 1 ; !( L = (6 * m_iCincture - 6 ) & ( L ( 6 * m_iCincture ); m_iCincture+) L0 = 6 * m_iCincture - 3; else if( L (6 * m_iCincture - 6 ) & ( L 0) double i = floor(Rad / PI / 2); Rad -= 2 * PI * i; if(Rad PI) Rad -= 2 * PI; return Rad;高斯投影C+

11、源码 2/ GaussBL2xy.cpp : Defines the entry point for the console application. / #include stdafx.h #include math.h #include CoorTrans.h #include using namespace std; void main(int argc, char* argv) double MyL0 = 108; /中央子午线 double MyB = 33.44556666; /33 du 44 fen 55.6666 miao double MyL = 109.22334444;

12、 /3度带，109 d 22 m 33.4444 s PrjPoint_Krasovsky MyPrj; MyPrj.SetL0(MyL0); MyPrj.SetBL(MyB, MyL); double OutMyX; /结果应该大致是:3736714.783, 627497.303 double OutMyY; OutMyX = MyPrj.x; /正算结果: 北坐标x OutMyY = MyPrj.y; /结果：东坐标y /反算/ double InputMyX = 3736714.783; /如果是独立计算，应该给出中央子午线L0 double InputMyY = 627497.303

13、; MyPrj.Setxy(InputMyX, InputMyY); MyPrj.GetBL(&MyPrj.B, &MyPrj.L); /把计算出的BL的弧度值换算为dms形式 double OutputMyB; double OutputMyL; OutputMyB = MyPrj.B; /反算结果：B OutputMyL = MyPrj.L; /反算结果：L /分析表明，此程序的结果和Coord4.2的转换结果是一样的，只差到毫米级 /原程序有几个问题，1.Pi的值不对。2.SetBL中多了两行错误代码 double Dms2Rad(double Dms) double Degree, M

14、iniute; double Second; int Sign; double Rad; if(Dms = 0) Sign = 1; else Sign = -1; Dms = fabs(Dms); Degree = floor(Dms); Miniute = floor(fmod(Dms * 100.0, 100.0); Second = fmod(Dms * 10000.0, 100.0); Rad = Sign * (Degree + Miniute / 60.0 + Second / 3600.0) * PI / 180.0; return Rad; double Rad2Dms(do

15、uble Rad) double Degree, Miniute; double Second; int Sign; double Dms; if(Rad = 0) Sign = 1; else Sign = -1; Rad = fabs(Rad * 180.0 / PI); Degree = floor(Rad); Miniute = floor(fmod(Rad * 60.0, 60.0); Second = fmod(Rad * 3600.0, 60.0); Dms = Sign * (Degree + Miniute / 100.0 + Second / 10000.0); retur

16、n Dms; / / Definition of PrjPoint / bool PrjPoint:BL2xy() double X, N, t, t2, m, m2, ng2; double sinB, cosB; X = A1 * B * 180.0 / PI + A2 * sin(2 * B) + A3 * sin(4 * B) + A4 * sin(6 * B); sinB = sin(B); cosB = cos(B); t = tan(B); t2 = t * t; N = a / sqrt(1 - e2 * sinB * sinB); m = cosB * (L - L0); m

17、2 = m * m; ng2 = cosB * cosB * e2 / (1 - e2); /x,y的计算公式见孔祥元等主编武汉大学出版社2002年出版的控制测量学的第72页 /书的的括号有问题， ( 和 应该交换 x = X + N * t * (0.5 + (5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24.0 + (61 - 58 * t2 + t2 * t2) * m2 / 720.0) * m2) * m2); y = N * m * ( 1 + m2 * ( (1 - t2 + ng2) / 6.0 + m2 * ( 5 - 18 * t2 + t2 *

18、 t2 + 14 * ng2 - 58 * ng2 * t2 ) / 120.0); y += 500000; return true; bool PrjPoint:xy2BL() double sinB, cosB, t, t2, N ,ng2, V, yN; double preB0, B0; double eta; y -= 500000; B0 = x / A1; do preB0 = B0; B0 = B0 * PI / 180.0; B0 = (x - (A2 * sin(2 * B0) + A3 * sin(4 * B0) + A4 * sin(6 * B0) / A1; eta

19、 = fabs(B0 - preB0); while(eta 0.000000001); B0 = B0 * PI / 180.0; B = Rad2Dms(B0); sinB = sin(B0); cosB = cos(B0); t = tan(B0); t2 = t * t; N = a / sqrt(1 - e2 * sinB * sinB); ng2 = cosB * cosB * e2 / (1 - e2); V = sqrt(1 + ng2); yN = y / N; B = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * y

20、N * yN * yN * yN / 12.0 + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360.0) * V * V * t / 2; L = L0 + (yN - (1 + 2 * t2 + ng2) * yN * yN * yN / 6.0 + (5 + 28 * t2 + 24 * t2 * t2 + 6 * ng2 + 8 * ng2 * t2) * yN * yN * yN * yN * yN / 120.0) / cosB; return true; bool PrjPoint:SetL0(do

21、uble dL0) L0 = Dms2Rad(dL0); return true; bool PrjPoint:SetBL(double dB, double dL) B = Dms2Rad(dB); L = Dms2Rad(dL); /B = dB; /我靠，I wana say fuck /L = dL; /del it ! BL2xy(); return true; bool PrjPoint:GetBL(double *dB, double *dL) *dB = Rad2Dms(B); *dL = Rad2Dms(L); return true; bool PrjPoint:Setxy

22、(double dx, double dy) x = dx; y = dy; xy2BL(); return true; 高斯投影VB源码 1坐标转换代码double a, f, e2, e12; / 基本椭球参数 double A1, A2, A3, A4; / 用于计算X的椭球参数 Const PI = 3.14159265358979Const Y0 = 500000#Dim A1#, A2#, A3#, A4#, a#, f#高斯正算求XPublic Function X(B#, L#, L0#, zbx%) As Double / B纬度L经度 L0中央子午线zbcs zbxe2#

23、= 1 - (f - 1) / f) * (f - 1) / f)e12# = (f / (f - 1) * (f / (f - 1) - 1B = D_R(B): L = D_R(L): L0 = D_R(L0)Dim n, t, t2, m, m2, ng2Dim sinB, cosBX = A1 * B * 180# / PI + A2 * Sin(2 * B) + A3 * Sin(4 * B) + A4 * Sin(6 * B)sinB = Sin(B)cosB = Cos(B)t = Tan(B)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)

24、m = cosB * (L - L0)m2 = m * mng2 = cosB * cosB * e2 / (1 - e2)X = X + n * t * (0.5 + (5 - t2 + 9 * ng2 + 4 * ng2 * ng2) / 24# + (61 - 58 * t2 + t2 * t2) * m2 / 720#) * m2) * m2)End Function高斯正算求YPublic Function Y(B#, L#, L0#, zbx%) As Doublezbcs zbxe2# = 1 - (f - 1) / f) * (f - 1) / f)e12# = (f / (f

25、 - 1) * (f / (f - 1) - 1B = D_R(B): L = D_R(L): L0 = D_R(L0)Dim n, t, t2, m, m2, ng2Dim sinB, cosBsinB = Sin(B)cosB = Cos(B)t = Tan(B)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)m = cosB * (L - L0)m2 = m * mng2 = cosB * cosB * e2 / (1 - e2)Y = n * m * (1 + m2 * (1 - t2 + ng2) / 6# + m2 * (5 - 18 * t2

26、 + t2 * t2 + 14 * ng2 - 58 * ng2 * t2) / 120#)Y = Y + Y0End Function高斯反算求B（纬度）Public Function B(X#, Y#, L0#, zbx%) As Doublezbcs zbxe2# = 1 - (f - 1) / f) * (f - 1) / f)e12# = (f / (f - 1) * (f / (f - 1) - 1L0 = D_R(L0)Dim sinB, cosB, t, t2, n, ng2, V, yNDim preB0, B0Dim etaY = Y - Y0B0 = X / A1Do W

27、hile True preB0 = B0 B0 = B0 * PI / 180# B0 = (X - (A2 * Sin(2 * B0) + A3 * Sin(4 * B0) + A4 * Sin(6 * B0) / A1 eta = Abs(B0 - preB0) If eta 0.0000001 Then Exit DoLoopB0 = B0 * PI / 180#sinB = Sin(B0)cosB = Cos(B0)t = Tan(B0)t2 = t * tn = a / Sqr(1 - e2 * sinB * sinB)ng2 = cosB * cosB * e2 / (1 - e2)V = Sqr(1 + ng2)yN = Y / nB = B0 - (yN * yN - (5 + 3 * t2 + ng2 - 9 * ng2 * t2) * yN * yN * yN * yN / 12# + (61 + 90 * t2 + 45 * t2 * t2) * yN * yN * yN * yN * yN * yN / 360#) * V * V * t / 2B = R_D(B)End Function高斯反算求L（经度）Public Function L(X#, Y#, L0#, zbx%) As Doublezbcs zbx