经纬度转换为xy，xy转换为经纬度

#include "iostream" #include "math.h" #define PI 3.1415926; void GuaussProj_6(double latd_d,double lontd_d,double &Guss6_x,double &Guss6_y) //高斯坐标6度带正算 { double a = 6378137; //椭球长轴，1975年参数 没有用现在元数据中给的a=6378245 const double e1_2 = 0.0066943799; //第一偏心率e1的平方 double e2_2 = 0.0067394967; //第二偏心率e2的平方 double z = 0.0174532925; // pi/180 double ro = 57.2957795; // 180/pi //const double r = 206265; double p = 1e-10; double A = 1.0050517739; //A,B,C,D为沿经线弧长公式的系数 double B = 0.00506237764; double C = 0.0000106245; double D = 0.00000002081; int zone; //带号 int center_lontd; //中央经线 double X, Y; //高斯坐标 double N; //卯酉圈半径 double eta_2; //eta的平方，计算时需用到的参数 double s; //沿经线方向的弧长 //double lontd_d; //角度表示的经度 //double latd_d; //角度表示的纬度 double latd_r; //弧度表示的纬度 double lontd_cha; //所给经度与中央经线的经差 double lontd_2, lontd_3, lontd_4, lontd_5; //经差的三、四、五次方 double sin_latd, cos_latd, tan_latd, cos_latd_3, cos_latd_5, tan_latd_4; //纬度的三角函数值 latd_r = latd_d/ro;//纬度弧度制 // latd_r = latd_d;//纬度弧度制 // latd_d =latd_d*ro ;//纬度角度制 // lontd_d =lontd_d*ro ;//经度角度值 if (lontd_d/6 != 0) //求带号 zone = int(lontd_d/6) + 1; else zone = int(lontd_d/6); center_lontd = 6*zone - 3; //求中央经线 lontd_cha = lontd_d - center_lontd; //求经差 lontd_cha = lontd_cha * z; //经差转换为弧度 lontd_2 = lontd_cha*lontd_cha; lontd_3 = lontd_cha*lontd_cha*lontd_cha; lontd_4 = lontd_cha*lontd_cha*lontd_cha*lontd_cha; lontd_5 = lontd_cha*lontd_cha*lontd_cha*lontd_cha*lontd_cha; sin_latd = sin(latd_r); cos_latd = cos(latd_r); tan_latd = tan(latd_r); cos_latd_3 = cos_latd*cos_latd*cos_latd; //三次方 cos_latd_5 = cos_latd*cos_latd*cos_latd*cos_latd*cos_latd; //五次方 tan_latd_4 = tan_latd*tan_latd*tan_latd*tan_latd; N = a / sqrt((1 - e1_2 * sin_latd *sin_latd)); //计算卯酉圈半径 eta_2 = e2_2 * cos_latd * cos_latd; s = a * (1 - e1_2) * ((A/ro)*latd_d - (B/2)*sin(2*latd_r) + (C/4)*sin(4*latd_r) - (D/6)*sin(6*latd_r)); X = s + (lontd_2*N/2)*sin_latd*cos_latd + (lontd_4*N/24)*sin_latd*cos_latd_3 * (5 - tan_latd*tan_latd + 9*eta_2 + 4*eta_2*eta_2); Y = lontd_cha*N*cos_latd + lontd_3*N*cos_latd_3*(1 - tan_latd*tan_latd + eta_2)/6 + lontd_5*N*cos_latd_5*(5 - 18*tan_latd*tan_latd + tan_latd_4); Y = Y + 500000; double temp_Y = zone * 1000000; //m_X = X; //m_Y = temp_Y + Y; //在Y坐标前加上带号 Guss6_x=X; Guss6_y=Y; } /*---------------- '函数名称：Trf_XYtoBL '功能:高斯克吕格坐标系统反算 '接口：[IN]x,y，高斯坐标 (按值传递) '[IN]L0,中央子午线,CorD 高斯坐标带类型(3度/6度) '[OUT]b,l 站点大地坐标值 (按地址传递) */ void xytobl(double x,double y,double &B ,double &L , double L0 , double CorD ) { double R= 6367558.496863 ;//'地球椭球半径 double e2 = 0.00673852415 ;//'第二偏心率值 double G2 , H2 , I2 , J2 , k2 , L2 , N2 , M2 ,O2 , P2 , Q2 , R2 , S2 , T2 , U2 , V2 , W2 , X2 , ABSY ,tmp1 , tmp2 , tmp3 , nDebtNo ; G2 = x/ R; H2 = cos(G2); I2 = sin(G2); J2 = I2 *I2; tmp1 = J2 * 2.38209 * pow(10 ,-7); tmp2 = J2 * (2.983718 * pow(10 ,-5)); tmp3 = I2 * H2 * (5.051773759 * pow(10 ,-3)); k2 = G2 + tmp3 - tmp2 - tmp1; L2 = cos(k2); M2 = e2 * L2 *L2; N2 = 1 + M2; O2 = 6399698.9018 / sqrt(N2); P2 = tan(k2); Q2 = P2 *P2; //nDebtNo = IIf(CorD = 6, ((L0 + 3) / 6), IIf(CorD = 3, (L0 / 3), 0)) nDebtNo =(CorD = 6, (L0 + 3) / 6)?((L0 + 3) / 6):(CorD = 3)?(L0 / 3): 0; R2 = nDebtNo * pow(10 ,6) + 5 * pow(10,5); if ( y > 500000 && y < nDebtNo * pow(10 ,6)) y = y - 500000 * 2; ABSY = abs(y); S2 = (ABSY < 500000)?(ABSY - 500000) / O2:(ABSY - R2) / O2; T2 = S2*S2; B = (k2 - ((((45 * Q2 + 90) * Q2 + 61) * T2 / 30 - (3 - 9 * M2) * Q2 - 5 - M2) * T2 / 12 + 1) * T2 * P2 * N2 / 2) * 180 / (4 * atan(1.000001)); L = ((((24 * Q2 + 28) * Q2 + (8 * Q2 + 6) * M2 + 5) * T2 / 20 - 2 * Q2 - N2) * T2 / 6 + 1) * S2 / L2 * 180 / (4 * atan(1.00000001)); if (y<0) L = L0 - L; else L = L0 + L; } void main() { int m=30,n=70;//---所求行列号 const int Ix = 100;//---图像行 const int Iy = 100;//---图像列 double x_resolution=3,y_resolution=3; double x_center; double y_center; int m_L,n_L; double G_X[Ix][Iy],G_Y[Ix][Iy]; double B,L,H,X,Y,Z; double L0=117, CorD=6 ; double B_center=40,L_center=116; GuaussProj_6(B_center,L_center,x_center,y_center); m_L=m-Ix/2; n_L=n-Iy/2; G_X[m][n]=x_center-(m_L)*x_resolution; G_Y[m][n]=y_center+(n_L)*y_resolution; xytobl(G_X[m][n],G_Y[m][n],B ,L , L0,CorD ); }