高斯-克吕格投影

简称 高斯投影

百度 可以得到 计算公式

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 }