史上最全的高斯投影源代码集合精心收集9种.docx
《史上最全的高斯投影源代码集合精心收集9种.docx》由会员分享,可在线阅读,更多相关《史上最全的高斯投影源代码集合精心收集9种.docx(43页珍藏版)》请在冰豆网上搜索。
史上最全的高斯投影源代码集合精心收集9种
高斯投影C++源码1
CCoordinateCEarthPlaneDialog:
:
EarthToPlane(constCEarth&earth)
{
doubleB=earth.m_dB.m_realRadians;
doubleL=earth.m_dL.m_realRadians;
L=SetRad_PIPI(L);
doubleL0=SetL0(L);
doubleX,N,t,t2,m,m2,ng2;
doublesinB,cosB;
X=A1*B*180.0/PI+A2*sin(2*B)+A3*sin(4*B)+A4*sin(6*B);
sinB=sin(B);
cosB=cos(B);
t=tan(B);
t2=t*t;
N=a/sqrt(1-e2*sinB*sinB);
m=cosB*(L-L0);
m2=m*m;
ng2=cosB*cosB*e2/(1-e2);
doublex=X+N*t*((0.5+((5-t2+9*ng2+4*ng2*ng2)/24.0+(61-58*t2+t2*t2)*m2/720.0)*m2)*m2);
doubley=N*m*(1+m2*((1-t2+ng2)/6.0+m2*(5-18*t2+t2*t2+14*ng2-58*ng2*t2)/120.0));
y+=500000;
CCoordinatecoor(x,y);
returncoor;
}
//将大地坐标系转换为平面直角坐标系
CEarthCEarthPlaneDialog:
:
PlaneToEarth(constCCoordinate&coor,UINTSign)
{
doublex=coor.m_dX;
doubley=coor.m_dY;
doubleL0=(6.0*Sign-3.0)/180*PI;
doublesinB,cosB,t,t2,N,ng2,V,yN;
doublepreB0,B0;
doubleeta;
y-=500000;
B0=x/A1;
do
{
preB0=B0;
B0=B0*PI/180.0;
B0=(x-(A2*sin(2*B0)+A3*sin(4*B0)+A4*sin(6*B0)))/A1;
eta=fabs(B0-preB0);
}while(eta>0.0000000000000001);
B0=B0*PI/180.0;
sinB=sin(B0);
cosB=cos(B0);
t=tan(B0);
t2=t*t;
N=a/sqrt(1-e2*sinB*sinB);
ng2=cosB*cosB*e2/(1-e2);
V=sqrt(1+ng2);
yN=y/N;
doubleB=B0-(yN*yN-(5+3*t2+ng2-9*ng2*t2)*yN*yN*yN*yN/12.0+(61+90*t2+45*t2*t2)*yN*yN*yN*yN*yN*yN/360.0)*V*V*t/2;
doubleL=L0+(yN-(1+2*t2+ng2)*yN*yN*yN/6.0+(5+28*t2+24*t2*t2+6*ng2+8*ng2*t2)*yN*yN*yN*yN*yN/120.0)/cosB;
CEarthearth(B,L,0,'r');
returnearth;
}
//将平面直角坐标系转换为大地坐标系
voidCEarthPlaneDialog:
:
OnRADIOPlane1954()
{
//TODO:
Addyourcontrolnotificationhandlercodehere
a=6378245;
f=298.3;
e2=1-((f-1)/f)*((f-1)/f);
e12=(f/(f-1))*(f/(f-1))-1;
A1=111134.8611;
A2=-16036.4803;
A3=16.8281;
A4=-0.0220;
m_iReferenceFrame=0;
}
voidCEarthPlaneDialog:
:
OnRADIOPlane1980()
{
//TODO:
Addyourcontrolnotificationhandlercodehere
a=6378140;
f=298.257;
e2=1-((f-1)/f)*((f-1)/f);
e12=(f/(f-1))*(f/(f-1))-1;
A1=111133.0047;
A2=-16038.5282;
A3=16.8326;
A4=-0.0220;
m_iReferenceFrame=1;
}
voidCEarthPlaneDialog:
:
OnRADIOPlaneWGS84()
{
//TODO:
Addyourcontrolnotificationhandlercodehere
a=6378137;
f=298.257223563;
e2=1-((f-1)/f)*((f-1)/f);
e12=(f/(f-1))*(f/(f-1))-1;
A1=111133.0047;
A2=-16038.5282;
A3=16.8326;
A4=-0.0220;
m_iReferenceFrame=2;
}
doubleCEarthPlaneDialog:
:
SetL0(doubleL)
{
L=L/PI*180;
if(L>0)
{
for(m_iCincture=1;!
((L>=(6*m_iCincture-6))&(L<(6*m_iCincture)));m_iCincture++)
{}
L0=6*m_iCincture-3;
}
elseif(L<0)
{
L=-L;
for(m_iCincture=1;!
((L>(6*m_iCincture-6))&(L<=(6*m_iCincture)));m_iCincture++)
{}
m_iCincture=61-m_iCincture;
L0=m_iCincture*6-363;
}
else
{
L0=3;
m_iCincture=1;
}
doubletemp=((double)L0)/180*PI;
returntemp;
}
doubleCEarthPlaneDialog:
:
SetRad_PIPI(doubleRad)
{
if(Rad>0)
{
doublei=floor(Rad/PI/2);
Rad-=2*PI*i;
}
if(Rad<0)
{
Rad=-Rad;
doublei=floor(Rad/PI/2);
Rad-=2*PI*i;
Rad=-Rad+PI*2;
}
if(Rad>PI)
Rad-=2*PI;
returnRad;
}
高斯投影C++源码2
//GaussBL2xy.cpp:
Definestheentrypointfortheconsoleapplication.
//
#include"stdafx.h"
#include"math.h"
#include"CoorTrans.h"
#include
usingnamespacestd;
voidmain(intargc,char*argv[])
{
doubleMyL0=108;//中央子午线
doubleMyB=33.44556666;//33du44fen55.6666miao
doubleMyL=109.22334444;//3度带,109d22m33.4444s
PrjPoint_KrasovskyMyPrj;
MyPrj.SetL0(MyL0);
MyPrj.SetBL(MyB,MyL);
doubleOutMyX;//结果应该大致是:
3736714.783,627497.303
doubleOutMyY;
OutMyX=MyPrj.x;//正算结果:
北坐标x
OutMyY=MyPrj.y;//结果:
东坐标y
//////////////////反算////////////////////////////////////////
doubleInputMyX=3736714.783;//如果是独立计算,应该给出中央子午线L0
doubleInputMyY=627497.303;
MyPrj.Setxy(InputMyX,InputMyY);
MyPrj.GetBL(&MyPrj.B,&MyPrj.L);//把计算出的BL的弧度值换算为dms形式
doubleOutputMyB;
doubleOutputMyL;
OutputMyB=MyPrj.B;//反算结果:
B
OutputMyL=MyPrj.L;//反算结果:
L
//分析表明,此程序的结果和Coord4.2的转换结果是一样的,只差到毫米级
//原程序有几个问题,1.Pi的值不对。
2.SetBL中多了两行错误代码
}
doubleDms2Rad(doubleDms)
{
doubleDegree,Miniute;
doubleSecond;
intSign;
doubleRad;
if(Dms>=0)
Sign=1;
else
Sign=-1;
Dms=fabs(Dms);
Degree=floor(Dms);
Miniute=floor(fmod(Dms*100.0,100.0));
Second=fmod(Dms*10000.0,100.0);
Rad=Sign*(Degree+Miniute/60.0+Second/3600.0)*PI/180.0;
returnRad;
}
doubleRad2Dms(doubleRad)
{
doubleDegree,Miniute;
doubleSecond;
intSign;
doubleDms;
if(Rad>=0)
Sign=1;
else
Sign=-1;
Rad=fabs(Rad*180.0/PI);
Degree=floor(Rad);
Miniute=floor(fmod(Rad*60.0,60.0));
Second=fmod(Rad*3600.0,60.0);
Dms=Sign*(Degree+Miniute/100.0+Second/10000.0);
returnDms;
}
///////////////////////////////////////////////////
//DefinitionofPrjPoint
///////////////////////////////////////////////////
boolPrjPoint:
:
BL2xy()
{
doubleX,N,t,t2,m,m2,ng2;
doublesinB,cosB;
X=A1*B*180.0/PI+A2*sin(2*B)+A3*sin(4*B)+A4*sin(6*B);
sinB=sin(B);
cosB=cos(B);
t=tan(B);
t2=t*t;
N=a/sqrt(1-e2*sinB*sinB);
m=cosB*(L-L0);
m2=m*m;
ng2=cosB*cosB*e2/(1-e2);
//x,y的计算公式见孔祥元等主编武汉大学出版社2002年出版的《控制测量学》的第72页
//书的的括号有问题,(和[应该交换
x=X+N*t*((0.5+((5-t2+9*ng2+4*ng2*ng2)/24.0+(61-
58*t2+t2*t2)*m2/720.0)*m2)*m2);
y=N*m*(1+m2*((1-t2+ng2)/6.0+m2*(5-18*t2+t2*t2
+14*ng2-58*ng2*t2)/120.0));
y+=500000;
returntrue;
}
boolPrjPoint:
:
xy2BL()
{
doublesinB,cosB,t,t2,N,ng2,V,yN;
doublepreB0,B0;
doubleeta;
y-=500000;
B0=x/A1;
do
{
preB0=B0;
B0=B0*PI/180.0;
B0=(x-(A2*sin(2*B0)+A3*sin(4*B0)+A4*sin(6*B0)))/A1;
eta=fabs(B0-preB0);
}while(eta>0.000000001);
B0=B0*PI/180.0;
B=Rad2Dms(B0);
sinB=sin(B0);
cosB=cos(B0);
t=tan(B0);
t2=t*t;
N=a/sqrt(1-e2*sinB*sinB);
ng2=cosB*cosB*e2/(1-e2);
V=sqrt(1+ng2);
yN=y/N;
B=B0-(yN*yN-(5+3*t2+ng2-9*ng2*t2)*yN*yN*yN*yN/
12.0+(61+90*t2+45*t2*t2)*yN*yN*yN*yN*yN*yN/360.0)
*V*V*t/2;
L=L0+(yN-(1+2*t2+ng2)*yN*yN*yN/6.0+(5+28*t2+24
*t2*t2+6*ng2+8*ng2*t2)*yN*yN*yN*yN*yN/120.0)/cosB;
returntrue;
}
boolPrjPoint:
:
SetL0(doubledL0)
{
L0=Dms2Rad(dL0);
returntrue;
}
boolPrjPoint:
:
SetBL(doubledB,doubledL)
{
B=Dms2Rad(dB);
L=Dms2Rad(dL);
//B=dB;//我靠,Iwanasayfuck
//L=dL;//delit!
BL2xy();
returntrue;
}
boolPrjPoint:
:
GetBL(double*dB,double*dL)
{
*dB=Rad2Dms(B);
*dL=Rad2Dms(L);
returntrue;
}
boolPrjPoint:
:
Setxy(doubledx,doubledy)
{
x=dx;
y=dy;
xy2BL();
returntrue;
}
高斯投影VB源码1
坐标转换代码
doublea,f,e2,e12;//基本椭球参数
doubleA1,A2,A3,A4;//用于计算X的椭球参数
ConstPI=3.14159265358979
ConstY0=500000#
DimA1#,A2#,A3#,A4#,a#,f#
'高斯正算求X
PublicFunctionX(B#,L#,L0#,zbx%)AsDouble//B纬度L经度L0中央子午线
zbcszbx
e2#=1-((f-1)/f)*((f-1)/f)
e12#=(f/(f-1))*(f/(f-1))-1
B=D_R(B):
L=D_R(L):
L0=D_R(L0)
Dimn,t,t2,m,m2,ng2
DimsinB,cosB
X=A1*B*180#/PI+A2*Sin(2*B)+A3*Sin(4*B)+A4*Sin(6*B)
sinB=Sin(B)
cosB=Cos(B)
t=Tan(B)
t2=t*t
n=a/Sqr(1-e2*sinB*sinB)
m=cosB*(L-L0)
m2=m*m
ng2=cosB*cosB*e2/(1-e2)
X=X+n*t*((0.5+((5-t2+9*ng2+4*ng2*ng2)/24#+(61-58*t2+t2*t2)*m2/720#)*m2)*m2)
EndFunction
'高斯正算求Y
PublicFunctionY(B#,L#,L0#,zbx%)AsDouble
zbcszbx
e2#=1-((f-1)/f)*((f-1)/f)
e12#=(f/(f-1))*(f/(f-1))-1
B=D_R(B):
L=D_R(L):
L0=D_R(L0)
Dimn,t,t2,m,m2,ng2
DimsinB,cosB
sinB=Sin(B)
cosB=Cos(B)
t=Tan(B)
t2=t*t
n=a/Sqr(1-e2*sinB*sinB)
m=cosB*(L-L0)
m2=m*m
ng2=cosB*cosB*e2/(1-e2)
Y=n*m*(1+m2*((1-t2+ng2)/6#+m2*(5-18*t2+t2*t2+14*ng2-58*ng2*t2)/120#))
Y=Y+Y0
EndFunction
'高斯反算求B(纬度)
PublicFunctionB(X#,Y#,L0#,zbx%)AsDouble
zbcszbx
e2#=1-((f-1)/f)*((f-1)/f)
e12#=(f/(f-1))*(f/(f-1))-1
L0=D_R(L0)
DimsinB,cosB,t,t2,n,ng2,V,yN
DimpreB0,B0
Dimeta
Y=Y-Y0
B0=X/A1
DoWhileTrue
preB0=B0
B0=B0*PI/180#
B0=(X-(A2*Sin(2*B0)+A3*Sin(4*B0)+A4*Sin(6*B0)))/A1
eta=Abs(B0-preB0)
Ifeta<0.0000001ThenExitDo
Loop
B0=B0*PI/180#
sinB=Sin(B0)
cosB=Cos(B0)
t=Tan(B0)
t2=t*t
n=a/Sqr(1-e2*sinB*sinB)
ng2=cosB*cosB*e2/(1-e2)
V=Sqr(1+ng2)
yN=Y/n
B=B0-(yN*yN-(5+3*t2+ng2-9*ng2*t2)*yN*yN*yN*yN/12#+(61+90*t2+45*t2*t2)*yN*yN*yN*yN*yN*yN/360#)*V*V*t/2
B=R_D(B)
EndFunction
'高斯反算求L(经度)
PublicFunctionL(X#,Y#,L0#,zbx%)AsDouble
zbcszbx