VC 系统下,地球坐标转换的源代码,在WGS-84坐标和北京54坐标之间的一套转换参数,可以全国通用的，在每个地方会不一样，因为它们是两个不同的椭球基准。-VC system, the Earth coordinates conversion source code, the WGS-84 coordinates and coordinate between Beijing 54 set conversion parameters can be nationwide, in every place to be different, because they are two different ellipsoid benchmarks.
UTM投影正转算法实现
在Gissky上发现了戴勤奋老师关于坐标转换的程序以及公式，程序的界面和转换精度都值得称赞，只是可惜没有提供实现的源码，最近由于工作需要正好需要实现这部分功能，通过对戴勤奋老师提供的公式以及OGP（欧洲石油勘探组织）最新的《Coordinate Conversions and Transformations including Formulas》文档进行了研究，利用程序将UTM投影正反转实现（反转程序在下一篇中介绍），现发布出来供大家共同学习，也希望大家指出程序中的不足。
/****************************************************************************
** 函数名：LLtoUTM
** 输入： a 椭球体长半轴
** f 椭球体扁率 f=(a-b)/a 其中b代表椭球体的短半轴
** Lat 经过UTM投影之前的纬度
** Long 经过UTM投影之前的经度
** LongOrigin 中央经度线
** FN 纬度起始点，北半球为0，南半球为10000000.0m
** UTMNorthing 经过UTM投影后的纬度方向的坐标
** UTMEasting 经过UTM投影后的经度方向的坐标
** 功能描述：UTM投影正转
** 作者： CiCi
** 单位： 中国地质大学（北京）地球科学与资源学院
** 创建日期：2007年3月28日
** 版本：1.0
***************************************************************************/
LLtoUTM(double a,double f,const double Lat, const double Long, double LongOrigin, double FN,double &UTMNorthing,double &UTMEasting)
{
//e表示WGS84第一偏心率,eSquare表示e的平方,
double eSquare =(2*f-f*f) ;
double k0 = 0.9996;
double e2Square;
double V, T, C, A, M;
//确保longtitude位于-180.00----179.9之间
double LongTemp = (Long+180)-int((Long+180)/360)*360-180;
double LatRad = Lat*PI/180;
double LongRad = LongTemp*PI/180;
double LongOriginRad;
LongOriginRad = LongOrigin * PI/180;
e2Square = (eSquare)/(1-eSquare);
V = a/sqrt(1-eSquare*sin(LatRad)*sin(LatRad));
T = tan(LatRad)*tan(LatRad);
C = e2Square*cos(LatRad)*cos(LatRad);
A = cos(LatRad)*(LongRad-LongOriginRad);
M = a*((1-eSquare/4-3*eSquare*eSquare/64-5*eSquare*eSquare*eSquare/256)*LatRad-(3*eSquare/8+3*eSquare*eSquare/32+45*eSquare*eSquare*eSquare/1024)*sin(2*LatRad)+(15*eSquare*eSquare/256+45*eSquare*eSquare*eSquare/1024)*sin(4*LatRad) -(35*eSquare*eSquare*eSquare/3072)*sin(6*LatRad));
UTMEasting = (double)(k0*V*(A+(1-T+C)*A*A*A/6
+ (5-18*T+T*T+72*C-58*e2Square)*A*A*A*A*A/120)+ 500000.0);
UTMNorthing = (double)(k0*(M+V*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+ (61-58*T+T*T+600*C-330*e2Square)*A*A*A*A*A*A/720)));
//南半球纬度起点为10000000.0m
UTMNorthing=UTMNorthing+FN;
}
-------------------------------------
平常的工作中就需要用到坐标转换的问题，特别是我做的webgis，由于需要在web页面上输入经纬度，而在实际上中采用国家标准的大地坐标来作图，所以就需要考虑坐标转换，否则人家输入到经纬度，想定位，结果坐标无法匹配，而最终无法实现功能，客户肯定会发牢骚的。这个是我开发来转换坐标的部分代码，结果很准确的，项目中应用过了的。
/// <summary>
/// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
/// </summary>
/// <param name="longitude">经度,以度为单位</param>
/// <param name="latitude">纬度,以度为单位</param>
/// <returns></returns>
public Point3D GaussProjCal(double longitude, double latitude)
{
Point3D pnt = new Point3D();
int ProjNo = 0; int ZoneWide; ////带宽
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double e2, ee, NN, T, C, A, M;
ZoneWide = paraDegree;
ProjNo = (int)(longitude / ZoneWide);
longitude0 = originLongitude;
longitude1 = Angle.DmsToRad(longitude);//经度转换为弧度
latitude1 = Angle.DmsToRad(latitude);//纬度转换为弧度
e2 = 2 * base.flat - flat * flat;
ee = e2 * (1.0 - e2);
double sinB = Math.Sin(latitude1);
double tanB = Math.Tan(latitude1);
double cosB = Math.Cos(latitude1);
NN = longRadius / Math.Sqrt(1.0 - e2 * sinB * sinB);
T = tanB * tanB;
C = ee * cosB * cosB;
A = (longitude1 - longitude0) * cosB;
double a0 = 1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256;
double a2 = 3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024;
double a4 = 15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024;
double a6 = 35 * e2 * e2 * e2 / 3072;
M = longRadius * (a0 * latitude1 - a2 * Math.Sin(2 * latitude1) + a4 * Math.Sin(4 * latitude1) - a6 * Math.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 * Math.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 + 500000L;
X0 = 500000L;
Y0 = 0;
xval = xval + X0;
yval = yval + Y0;
pnt.X = xval;
pnt.Y = yval;
return pnt;
}
/// <summary>
/// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
/// </summary>
/// <param name="lola"></param>
/// <returns></returns>
public Point3D GaussProjCal(LoLatude lola)
{
double lo = lola.ToLongitude();
double la = lola.ToLatitude();
return GaussProjCal(lo, la);
}
/// <summary>
/// 高斯投影由经纬度(Unit:DD)反算大地坐标(含带号，Unit:Metres)
/// </summary>
/// <param name="longitude">经度,以度为单位</param>
/// <param name="latitude">纬度,以度为单位</param>
/// <param name="x">X坐标</param>
/// <param name="y">Y坐标</param>
public void GaussProjCal(double longitude, double latitude,out double x,out double y)
{
Point3D pnt = GaussProjCal(longitude, latitude);
x = pnt.X;
y = pnt.Y;
}
/// <summary>
/// 高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)
/// </summary>
/// <param name="X">X坐标</param>
/// <param name="Y">Y坐标</param>
/// <param name="longitude">经度</param>
/// <param name="latitude">纬度</param>
public void GaussProjInvCal(double X, double Y, out double longitude, out double latitude)
{
int ProjNo; int ZoneWide; ////带宽
double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
double e1, e2,ee, NN, T, C, M, D, R, u, fai;
ZoneWide = paraDegree; ////6度带宽
ProjNo = (int)(X / 1000000L); //查找带号
longitude0 = originLongitude; //中央经线
//X0 = ProjNo * 1000000L + 500000L;
X0 = 500000L;
Y0 = 0;
xval = X - X0;
yval = Y - Y0; //带内大地坐标
e2 = 2 * flat - flat * flat;
double sqrtE2 = Math.Sqrt(1 - e2);
e1 = (1.0 - sqrtE2) / (1.0 + sqrtE2);
ee = e2 / (1 - e2);
M = yval;
double q0 = 1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256;
double q2 = 3 * e1 / 2 - 27 * e1 * e1 * e1 / 32;
double q4 = 21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32;
double q6 = 151 * e1 * e1 * e1 / 96;
double q8 = 1097 * e1 * e1 * e1 * e1 / 512;
u = M / (longRadius * q0);
fai = u + q2 * Math.Sin(2 * u) + q4 * Math.Sin(4 * u)
+ q6 * Math.Sin(6 * u) + q8 * Math.Sin(8 * u);
double cosB = Math.Cos(fai);
double tanB = Math.Tan(fai);
double sinB = Math.Sin(fai);
C = ee * cosB * cosB;
T = tanB * tanB;
NN = longRadius / Math.Sqrt(1.0 - e2 * sinB * sinB);
R = longRadius * (1 - e2) / Math.Sqrt((1 - e2 * sinB * sinB) * (1 - e2 * sinB * sinB) * (1 - e2 * sinB * sinB));
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) / cosB;
latitude1 = fai - (NN * tanB / R) * (D * D / 2 - (5 + 3 * T + 10 * C - 4 * 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 / Pi2;
latitude = latitude1 / Pi2;
}
/// <summary>
/// 大地坐标反算经纬度坐标
/// </summary>
/// <param name="x">X坐标</param>
/// <param name="y">Y坐标</param>
/// <returns></returns>
public LoLatude GaussProjInvCal(double x,double y)
{
double lo, la;
GaussProjInvCal(x, y, out lo, out la);
lo = lo * Pi2;
la = la * Pi2;
LoLatude lola = new LoLatude();
lola.Longitude = Angle.RadToDms(lo);
lola.Latitude = Angle.RadToDms(la);
return lola;
}
/// <summary>
/// 大地坐标反算经纬度坐标
/// </summary>
/// <param name="pnt">坐标点</param>
/// <returns></returns>
public LoLatude GaussProjInvCal(Point3D pnt)
{
double lo, la;
GaussProjInvCal(pnt.X, pnt.Y, out lo, out la);
lo = lo * Pi2;
la = la * Pi2;
LoLatude lola = new LoLatude();
lola.Longitude = Angle.RadToDms(lo);
lola.Latitude = Angle.RadToDms(la);
return lola;
}
------------------------------
//高斯投影正、反算
//////6度带宽 54年北京坐标系
// WGS84.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include <math.h>
//高斯投影由大地坐标(Unit:Metres)反算经纬度(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; 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));
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*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 longitude1,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 * 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*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;
*X = xval;
*Y = yval;
}
int main(int argc, char* argv[])
{
double wLong, wLat;
double bLong, bLat;
double x, y;
printf("input WGS Long-Lat:");
scanf("%lf,%lf", &wLong, &wLat);
printf("Your WGS84 Longitude-Latitude: %f, %f\n", wLong, wLat);
GaussProjCal(wLong, wLat, &x, &y);
printf("B-54 X,Y: %f, %f\n", x, y);
GaussProjInvCal(x, y, &bLong, &bLat);
printf("Result: %f, %f\n", bLong, bLat);
return 0;
}
NN卯酉圈曲率半径，测量学里面用N表示
M为子午线弧长，测量学里用大X表示 字串2
fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示 字串4
R为底点所对的曲率半径，测量学里用Nf表示
___________________________________
'************************************************************************
'高斯克吕格与经纬度坐标值转换代码
'Writen by Rodger Yuan 9 5 2006
'参考文献
'v0 0.1
'用于在经纬度坐标和高斯克吕格坐标之间的转换。
'高斯克吕格为一种投影，根据椭球体和基准面不同又有所区分，常用的北京54和西安80即
'采用这种投影方式，投影后的坐标为平面坐标系，单位为米
'现在参数的坐标系采用测绘坐标系，x为纵坐标，y为横坐标
'返回参数为自定义类型，双精度点
'调用转换函数前需要调用初始化过程进行初始化
'-------------------------------------------------------------------------------
'Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)
'说明: 用于初始化转换参数
'TuoqiuCanshu 枚举类型，提供北京54、西安80和WGS84三个椭球参数
'Daihao 整型 为高斯克吕格投影六度分带带号，取值为1～60
'-------------------------------------------------------------------------------
'Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal k0 As Double)
'说明: 用于进一步初始化转换参数（暂不提供）
'dE 东偏移
'dE 北偏移
'k0 比例因子
'-------------------------------------------------------------------------------
'Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD
'************************************************************************
'基本变量定义
Dim a As Double '椭球体长半轴
Dim b As Double '椭球体短半周
Dim f As Double '扁率
Dim e As Double '第一偏心率
Dim eq As Double '第二偏心率
Dim dh As Integer '带号
Dim FE As Double '东偏移
Dim FN As Double '北偏移
Dim L0 As Double '中央经度
Dim k0 As Double '比例因子
Const PI As Double = 3.14159265358979
Public Enum Canshu
Beijing54 = 0
Xian80 = 1
WGS84 = 2
End Enum
Type PointD
X As Double
Y As Double
End Type
Public Sub init(ByVal TuoqiuCanshu As Canshu, ByVal Daihao As Integer)
Select Case TuoqiuCanshu
'Krassovsky （北京54采用） 6378245 6356863.0188
'IAG 75（西安80采用） 6378140 6356755.2882
'WGS 84 6378137 6356752.3142
Case 0: '北京五四
a = 6378245
b = 6356863.0188
Case 1: '西安八零
a = 6378140
b = 6356755.2882
Case 2: 'WGS84
a = 6378137
b = 6356752.3142
End Select
f = (a - b) / a
e = Sqr(1 - (b / a) ^ 2)
eq = Sqr((a / b) ^ 2 - 1)
If Daihao < 1 Or Daihao > 60 Then Exit Sub
dh = Daihao
L0 = (6 * dh - 3) * PI / 180
k0 = 1
FE = 500000 + dh * 1000000
FN = 0
End Sub
Public Sub initEx(ByVal dE As Double, ByVal dN As Double, ByVal dk0 As Double)
End Sub
Public Function JWgetGK(ByVal W As Double, ByVal J As Double) As PointD
'给出经纬度坐标，转换为高克投影坐标
Dim BY As Double
Dim LX As Double
Dim TC As Double
Dim CC As Double
Dim AC As Double
Dim MC As Double
Dim NC As Double
Dim rx As Double
Dim ry As Double
Dim resultP As PointD
BY = W * PI / 180
LX = J * PI / 180
TC = Math.Tan(BY) ^ 2
CC = eq ^ 2 * Cos(BY) ^ 2
AC = (LX - L0) * Cos(BY)
MC = a * ((1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256) * BY - (3 * e ^ 2 / 8 + 3 * e ^ 4 / 32 + 45 * e ^ 6 / 1024) * Sin(2 * BY) + (15 * e ^ 4 / 256 + 45 * e ^ 6 / 1024) * Sin(4 * BY) - (35 * e ^ 6 / 3072) * Sin(6 * BY))
NC = a / Sqr(1 - e ^ 2 * (Sin(BY)) ^ 2)
rx = k0 * (MC + NC * Tan(BY) * (AC ^ 2 / 2 + (5 - TC + 9 * CC + 4 * CC ^ 2) * AC ^ 4 / 24) + (61 - 58 * TC + T ^ 2 + 270 * CC - 330 * TC * CC) * AC ^ 6 / 720)
ry = FE + k0 * NC * (AC + (1 - TC + CC) * AC ^ 3 / 6 + (5 - 18 * TC + TC ^ 2 + 14 * CC - 58 * TC * CC) * AC ^ 5 / 120)
resultP.X = rx
resultP.Y = ry
JWgetGK = resultP
End Function
Public Function GKgetJW(ByVal X As Double, ByVal Y As Double) As PointD
'给出高克投影坐标，转换为经纬度坐标
Dim BY As Double
Dim LX As Double
Dim e1 As Double
Dim FI As Double
Dim Mf As Double
Dim Bf As Double
Dim Tf As Double
Dim Cf As Double
Dim Nf As Double
Dim Rf As Double
Dim D As Double
Dim RW As Double '纬度
Dim RJ As Double '经度
Dim resultP As PointD
YE = Y
XN = X
e1 = (1 - b / a) / (1 + b / a)
Mf = (XN - FN) / k0
FI = Mf / (a * (1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256))
Bf = FI + (3 * e1 / 2 - 27 * e1 ^ 3 / 32) * Sin(2 * FI) + (21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32) * Sin(4 * FI) + (151 * e1 ^ 3 / 96) * Sin(6 * FI)
Tf = Tan(Bf) ^ 2
Cf = eq ^ 2 * Cos(Bf) ^ 2
Nf = a / Sqr(1 - e ^ 2 * Sin(Bf) ^ 2)
Rf = a * (1 - e ^ 2) / Sqr((1 - e ^ 2 * Sin(Bf) ^ 2) ^ 3)
D = (YE - FE) / (k0 * Nf)
RW = Bf - (Nf * Tan(Bf) / Rf) * (D ^ 2 / 2 - (5 + 3 * Tf + Cf - 9 * Tf * Cf) * D ^ 4 / 24 + (61 + 90 * Tf + 45 * Tf ^ 2) * D ^ 6 / 720)
RJ = L0 + 1 / Cos(Bf) * (D - (1 + 2 * Tf + Cf) * D ^ 3 / 6 + (5 + 28 * Tf + 6 * Cf + 8 * Tf * Cf + 24 * Tf ^ 2) * D ^ 5 / 120)
resultP.X = RW * 180 / PI
resultP.Y = RJ * 180 / PI
GKgetJW = resultP
End Function
Sub demo()
init 0, 20
Dim p As PointD
p = JWgetGK(23.5, 120.5)
Debug.Print p.X, ",", p.Y
p = GKgetJW(p.X, p.Y)
Debug.Print p.X, ",", p.Y
End Sub