高斯投影正反算最好给出VC++程序源代码!-查字典问答网
分类选择

来自方流的问题

  高斯投影正反算最好给出VC++程序源代码!

  高斯投影正反算

  最好给出VC++程序源代码!

1回答
2020-10-13 05:51
我要回答
请先登录
盛志伟

  //高斯投影正、反算

  //////6度带宽54年北京坐标系

  //高斯投影由经纬度(Unit:DD)反算大地坐标(含带号,Unit:Metres)

  voidGaussProjCal(doublelongitude,doublelatitude,double*X,double*Y)

  {

  intProjNo=0;intZoneWide;////带宽

  doublelongitude1,latitude1,longitude0,latitude0,X0,Y0,xval,yval;

  doublea,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*ZoneWide+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-longitude0)*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*l

  atitude1));

  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;

  *X=xval;

  *Y=yval;

  }

  //高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)

  voidGaussProjInvCal(doubleX,doubleY,double*longitude,double*latitude)字串9

  {

  intProjNo;intZoneWide;////带宽

  doublelongitude1,latitude1,longitude0,latitude0,X0,Y0,xval,yval;

  doublee1,e2,f,a,ee,NN,T,C,M,D,R,u,fai,iPI;

  iPI=0.0174532925199433;////3.1415926535898/180.0;

  a=6378245.0;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-sqrt(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*sin(fai)*sin(fai));字串1

  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

2020-10-13 05:54:21

最新问答

推荐文章

猜你喜欢

附近的人在看

推荐阅读

拓展阅读

  • 大家都在看
  • 小编推荐
  • 猜你喜欢
  •