标签：克吕格 高斯 54 半轴 double 投影 tanB pow e1
高斯-克吕格投影
简称 高斯投影
百度 可以得到 计算公式
https://wenku.baidu.com/view/a02f613183c4bb4cf7ecd181.html
大致如下：
代码
https://blog.csdn.net/mysonismysun/article/details/8802437
1 #include <cmath>
2
3 //高斯平面坐标系
4 struct CRDCARTESIAN
5 {
6 double x;
7 double y;
8 double z;
9 } ;
10
11 //大地坐标系（可以是 北京54坐标系，西安80坐标系，WGS84坐标系（GPS 坐标））
12 struct CRDGEODETIC
13 {
14 double longitude; //经度
15 double latitude; //纬度
16 double height; //大地高,可设为0
17 };
18
19 #define WGS84 84 //WGS84坐标系（GPS 坐标）
20 #define BJ54 54 //北京54坐标系
21 #define XIAN80 80 //西安80坐标系
22 #define ZONEWIDE3 3 //投影带宽度 3
23 #define ZONEWIDE6 6 //投影带宽度 6
24
25 //---------------------------------------------------------------------------
26 void BLTOXY(CRDCARTESIAN * pcc, CRDGEODETIC * pcg, int Datum, int zonewide)
27 {
28 double B = pcg->latitude; //纬度
29 double L = pcg->longitude; //经度//纬度度数
30
31 double L0; //中央经线度数
32 double N; //卯酉圈曲率半径
33 double q2;
34 double x; //高斯平面纵坐标
35 double y; //高斯平面横坐标
36 double s; //赤道至纬度B的经线弧长
37 double f; //参考椭球体扁率
38 double e1; //椭球第一偏心率
39 double a; //参考椭球体长半轴
40 //double b; //参考椭球体短半轴
41 double a1, a2, a3, a4;
42 double b1, b2, b3, b4;
43 double c0, c1, c2, c3;
44
45 const double IPI = 0.0174532925199433333333; //3.1415926535898/180.0
46
47 int prjno = 0; //投影带号
48 // zonewide 投影带宽带， 3 或者是 6
49 if (zonewide == 6)
50 {
51 prjno = (int) (L / zonewide) + 1;
52 L0 = prjno * zonewide - 3;
53 }
54 else
55 {
56 prjno = (int) ((L - 1.5) / 3) + 1;
57 L0 = prjno * 3;
58 }
59
60 /*
61 * 北京 54
62 * 长半轴a=6378245m
63 * 短半轴b=6356863.0188m
64 * 扁率α=1/298.3
65 * 第一偏心率平方 =0.006693421622966
66 * 第二偏心率平方 =0.006738525414683
67 *
68 * 西安80
69 * 长半轴a=6378140±5（m）
70 * 短半轴b=6356755.2882m
71 * 扁率α=1/298.257
72 * 第一偏心率平方 =0.00669438499959
73 * 第二偏心率平方=0.00673950181947
74 *
75 * WGS84
76 * 长半轴a=6378137± 2（m）
77 * 短半轴b=6356752.3142m
78 * 扁率α=1/298.257223563
79 * 第一偏心率平方 =0.00669437999013
80 * 第二偏心率平方 =0.00673949674223
81 *
82 */
83
84 //Datum 投影基准面类型：北京54基准面为54，西安80基准面为80，WGS84基准面为84
85 if (Datum == 84)
86 {
87 a = 6378137;
88 f = 1 / 298.257223563;
89 }
90 else if (Datum == 54)
91 {
92 a = 6378245;
93 f = 1 / 298.3;
94 }
95 else if (Datum == 80)
96 {
97 a = 6378140;
98 f = 1 / 298.257;
99 }
100
101 e1 = 2 * f - f*f; //(a*a-b*b)/(a*a) 椭球第一偏心率
102
103 L0 = L0*IPI; // 转为弧度
104 L = L*IPI; // 转为弧度
105 B = B*IPI; // 转为弧度
106
107 double sinB = sin(B); //sinB
108 double cosB = cos(B); //cosB
109 double tanB = tan(B); //tanB
110
111 double l = L - L0; //L-L0l
112 double m = l * cosB; //ltanB
113
114 N = a / sqrt(1 - e1 * pow(sinB, 2));
115 q2 = e1 / (1 - e1) * pow(cosB, 2);
116
117 a1 = 1 + 3.0 / 4.0 * e1 + 45.0 / 64.0 * pow(e1, 2) + 175.0 / 256.0 * pow(e1, 3)
118 + 11025.0 / 16384.0 * pow(e1, 4) + 43659.0 / 65536.0 * pow(e1, 5);
119
120 a2 = 3.0 / 4.0 * e1 + 15.0 / 16.0 * pow(e1, 2) + 525.0 / 512.0 * pow(e1, 3)
121 + 2205.0 / 2048.0 * pow(e1, 4) + 72765.0 / 65536.0 * pow(e1, 5);
122
123 a3 = 15.0 / 64.0 * pow(e1, 2) + 105.0 / 256.0 * pow(e1, 3) + 2205.0 / 4096.0
124 * pow(e1, 4) + 10359.0 / 16384.0 * pow(e1, 5);
125
126 a4 = 35.0 / 512.0 * pow(e1, 3) + 315.0 / 2048.0 * pow(e1, 4) + 31185.0 / 13072.0
127 * pow(e1, 5);
128 b1 = a1 * a * (1 - e1);
129 b2 = -1.0 / 2.0 * a2 * a * (1 - e1);
130 b3 = 1.0 / 4.0 * a3 * a * (1 - e1);
131 b4 = -1.0 / 6.0 * a4 * a * (1 - e1);
132 c0 = b1;
133 c1 = 2 * b2 + 4 * b3 + 6 * b4;
134 c2 = -(8 * b3 + 32 * b4);
135 c3 = 32 * b4;
136 s = c0 * B + cosB * (c1 * sinB + c2 * pow(sinB, 3) + c3 * pow(sinB, 5));
137
138 x = s + 0.5 * N * tanB * pow(m, 2) + 1.0 / 24.0 * (5 - pow(tanB, 2) + 9 * q2 + 4
139 * pow(q2, 2)) * N * tanB * pow(m, 4) + 1.0 / 720.0 * (61 - 58 * pow(tanB, 6))
140 * N * tanB * pow(m, 6);
141
142 y = N * m + 1.0 / 6.0 * (1 - pow(tanB, 2) + q2) * N * pow(m, 3) + 1.0 / 120.0
143 * (5 - 18 * tanB * tanB + pow(tanB, 4) - 14 * q2 - 58 * q2 * pow(tanB, 2)) * N * pow(m, 5);
144
145 y = y + 1000000 * prjno + 500000;
146 pcc->x = x;
147 pcc->y = y - 38000000;
148 pcc->z = 0;
149 }
150
151 double BLDistance(CRDGEODETIC *pcg1, CRDGEODETIC *pcg2, int Datum, int zonewide)
152 {
153 CRDCARTESIAN pcc1, pcc2;
154 BLTOXY(&pcc1, pcg1, Datum, zonewide);
155 BLTOXY(&pcc2, pcg2, Datum, zonewide);
156
157 double xdes = fabs(pcc1.x - pcc2.x);
158 double ydes = fabs(pcc1.y - pcc2.y);
159 double des = sqrt(xdes * xdes + ydes * ydes);
160
161 return des;
162 }
标签：克吕格,高斯,54,半轴,double,投影,tanB,pow,e1 来源： https://www.cnblogs.com/zhangzhiwei122/p/15912481.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。