史上最全的高斯投影源代码集合精心收集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
升级会员 