GPS坐标转大地坐标

  1 double B = 39.3926;
  2 double L = 117.4514;
  3 //double B = 0;
  4 //double L = 117;
  5 double L0 = 117;
  6 
  7 /*
  8 ************************************************************************************************************************
  9 * Description: This function convert longitude/latitude(xx度xx分xxxxxx秒) to longitude/latitude(xx.xxxxxxxx度).
 10 *
 11 * Arguments  :     gps_x_y
 12                 longitude/latitude(xx度mm分ssssss秒) 输入形式为： xx.mmssssss
 13 *
 14 * Note(s)    : none
 15 ************************************************************************************************************************
 16 */
 17 double CALC_GPS_angle(float gps_x_y)
 18 {
 19     /*
 20     double min,sec;
 21     min = (floor(gps_x_y*100)-floor(gps_x_y)*100)/60;
 22     sec = (gps_x_y*10000-floor(gps_x_y*100)*100)/3600;
 23     printf("min = %.6lf\r\n",min);
 24     printf("sec = %.6lf\r\n",sec);
 25     */
 26     return (floor(gps_x_y)+(floor(gps_x_y*100)-floor(gps_x_y)*100)/60+(gps_x_y*10000-floor(gps_x_y*100)*100)/3600);
 27 }
 28 
 29 /*
 30 ************************************************************************************************************************
 31 * Description: 角度转换为弧度.
 32 *
 33 * Arguments  : angle 
 34                 输入形式(xx.xxxxxxxx度)
 35 *
 36 * Note(s)    : none
 37 ************************************************************************************************************************
 38 */
 39 double Angle2RAD(double angle)
 40 {
 41     return (angle*3.1415926/180);    
 42 }
 43 
 44 /*
 45 ************************************************************************************************************************
 46 * Description:     该函数调用以上2个函数完成计算功能.
 47 *
 48 * Arguments  :    input L :    经度longitude -> (xx度mm分ssssss秒) 输入形式为： xx.mmssssss
 49                 input B :    纬度latitude -> (xx度mm分ssssss秒) 输入形式为： xx.mmssssss
 50                 input L0:    中央子午线Central meridian -> (xx度mm分ssssss秒) 输入形式为： xx.mmssssss
 51                 output X:    坐标系的y轴
 52                 output Y    坐标系的x轴
 53 *
 54 * Note(s)    :     none
 55 ************************************************************************************************************************
 56 */
 57 void GpsGrid2GaussGrid(double L0, double B, double L, double* X, double* Y)
 58 {
 59     double E2,F2,G2,H2,I2,J2,K2,L2,M2,N2,O2,P2,Q2,R2;
 60     
 61     E2 = CALC_GPS_angle(B);
 62     #ifdef PRINT_DEG
 63     printf("E2 = %.9lf\r\n",E2);
 64     #endif
 65     F2 = CALC_GPS_angle(L);
 66     #ifdef PRINT_DEG
 67     printf("F2 = %.9lf\r\n",F2);
 68     #endif
 69     G2 = F2 - CALC_GPS_angle(L0);
 70     #ifdef PRINT_DEG
 71     printf("2G = %.9lf\r\n",G2);
 72     #endif
 73     H2 = 0.01744444444*G2;
 74     #ifdef PRINT_DEG
 75     printf("H2 = %.9lf\r\n",H2);
 76     #endif
 77     I2 = tan(Angle2RAD(fabs(E2)));
 78     #ifdef PRINT_DEG
 79     printf("I2 = %.9lf\r\n",I2);
 80     #endif
 81     J2 = cos(Angle2RAD(fabs(E2)));
 82     #ifdef PRINT_DEG
 83     printf("J2 = %.9lf\r\n",J2);
 84     #endif
 85     K2 = 0.006738525415*J2*J2;
 86     #ifdef PRINT_DEG
 87     printf("K2 = %.9lf\r\n",K2);
 88     #endif
 89     L2 = I2*I2;
 90     #ifdef PRINT_DEG
 91     printf("L2 = %.9lf\r\n",L2);
 92     #endif
 93     M2 = 1 + K2;
 94     #ifdef PRINT_DEG
 95     printf("M2 = %.9lf\r\n",M2);
 96     #endif
 97     N2 = 6399698.9018/sqrt(M2);
 98     #ifdef PRINT_DEG
 99     printf("N2 = %.9lf\r\n",N2);
100     #endif
101     O2 = H2*H2*J2*J2;
102     #ifdef PRINT_DEG
103     printf("O2 = %.9lf\r\n",O2);
104     #endif
105     P2 = I2*J2;
106     #ifdef PRINT_DEG
107     printf("P2 = %.9lf\r\n",P2);
108     #endif
109     Q2 = P2*P2;
110     #ifdef PRINT_DEG
111     printf("Q2 = %.9lf\r\n",Q2);
112     #endif
113     R2 = (32005.78006+Q2*(133.92133+Q2*0.7031));
114     #ifdef PRINT_DEG
115     printf("R2 = %.9lf\r\n",R2);
116     #endif
117     
118     *X = 6367558.49686*E2/57.29577951308-P2*J2*R2+((((L2-58)*L2+61)* O2/30+(4*K2+5)*M2-L2)*O2/12+1)*N2*I2*O2/2;
119     *Y = ((((L2-18)*L2-(58*L2-14)*K2+5)*O2/20+M2-L2)*O2/6+1)*N2*(H2*J2)+500000;
120 
121 }