高斯投影正反算C#代码

1、高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系)，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。/6度带宽54北京坐标系/高斯投影由大地坐标(Unit:Metres

2、)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude)int ProjNo; int ZoneWide; /带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; /3.1415926535898/180.0;a = 6378245.0

3、; f = 1.0/298.3; /54年北京坐标系参数/a=6378140.0; f=1/298.257; /80年西安坐标系参数ZoneWide = 6; /6度带宽ProjNo = (int)(X/1000000L) ; /查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; /中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; /带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sq

4、rt(1-e2)/(1.0+sqrt(1-e2);ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256);fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*

5、sin(fai)*sin(fai);R = a*(1-e2)/sqrt(1-e2*sin(fai)*sin(fai)*(1-e2*sin(fai)*sin(fai)*(1-e2*sin(fai)*sin(fai);D = xval/NN;/计算经度(Longitude)纬度(Latitude)longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D*D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4

6、*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);/转换为度DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;/高斯投影由经纬度(Unit:DD)正算平面坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y)int ProjNo=0; int ZoneWide; /带宽double longitu

7、de1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; /3.1415926535898/180.0;ZoneWide = 6; /6度带宽a=6378245.0; f=1.0/298.3; /54年北京坐标系参数/a=6378140.0; f=1/298.257; /80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWi

8、de + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; /经度转换为弧度latitude1 = latitude * iPI ; /纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1);T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude

9、0)*cos(latitude1);M=a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2*e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35*e2*e2*e2/3072)*sin(6*latitude1);xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 =1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;