gps.txt

WGS84下大地坐标和UTM坐标系统转换

UserInput_Lon = InputBox("请输入你所需地点的经度值"&Chr(13) & Chr(10)&Chr(13) & Chr(10)&"格式：dd.dddddd")
UserInput_Lat = InputBox("请输入你所需地点的纬度值"&Chr(13) & Chr(10)&Chr(13) & Chr(10)&"格式：dd.dddddd")
'计算用户UTM区域号
UTM_Zone=((Int(UserInput_Lon)-(Int(UserInput_Lon) Mod 6)+180)/6+1)
UTM_Zone_Center_Lon=UTM_Zone*6-180-3
'以下定义了WGS84椭球的参数
PI=3.1415926
WGS84_a=6378137
WGS84_1e2=0.00669437999013
WGS84_2e2=0.006739496742227
WGS84_N=WGS84_a/((1-WGS84_1e2*(sin(UserInput_Lat*PI/180))^2)^(1/2))
WGS84_m=cos(UserInput_Lat*PI/180)*(UserInput_Lon-UTM_Zone_Center_Lon)*PI/180
WGS84_nita2=(WGS84_2e2)*(cos(UserInput_Lat*PI/180)^2)
WGS84_t=tan(UserInput_Lat*PI/180)
'采用IUGG1975椭球计算子午线长度，如有WGS84参数则将下式替换为WGS84数值将会更精确
IUGG1975_PRO1=111133.0047*UserInput_Lat
IUGG1975_PRO2=16038.5282*sin(2*UserInput_Lat*PI/180)
IUGG1975_PRO3=16.8326*sin(4*UserInput_Lat*PI/180)
IUGG1975_PRO4=0.0220*sin(6*UserInput_Lat*PI/180)
IUGG1975_X=IUGG1975_PRO1-IUGG1975_PRO2+IUGG1975_PRO3-IUGG1975_PRO4
'计算WGS84与IUGG1975混合椭球下的UTM坐标
POR1=(WGS84_m^2)/2
PRO2=(5-WGS84_t^2+9*WGS84_nita2+4*WGS84_nita2*WGS84_nita2)*(WGS84_m^4)/24
PRO3=(6*(UserInput_Lon-UTM_Zone_Center_Lon)-58*WGS84_t^2+WGS84_t^4)*WGS84_m^6/720
PRO4=(1-WGS84_t^2+WGS84_nita2)*WGS84_m^3/6
PRO5=(5-18*WGS84_t^2+WGS84_t^4+14*WGS84_nita2-58*WGS84_nita2*WGS84_t^2)*WGS84_m^5/120
'计算出的用户输入点的UTM坐标
UserInput_Y=IUGG1975_X+WGS84_N*WGS84_t*(PRO1+PRO2+PRO3)
UserInput_X=WGS84_N*(WGS84_m+PRO4+PRO5)+500000
MyVar=MsgBox("UTM_Zone="&UTM_Zone&" UTM_X="&UserInput_X&" UTM_Y="&UserInput_Y,0,"计算结果：")