200分求经纬度BL换算到高斯平面直角坐标XY（高斯投影正反算）的源码

[imho888]

200分求经纬度BL换算到高斯平面直角坐标XY（高斯投影正反算）的源码

[hongqi162]

经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法

在 gis 的帖子"用excel完成gps坐标转换的简易方法 " 基础上,
我给出对应的vb程序段,我在 evb 开发的gps定位功能中，用它实现坐标换算(具体的参数请参考gis 的帖子)。
感觉速度比较快，效果比较好。所以帖上来，希望与名位交流：
=====================================
'经纬度bl换算到高斯平面直角坐标xy（高斯投影正算）

private function bl2xy(byref a2 as double, byref f2 as double, byref e2 as double, _
byref s2 as double, byref t2 as double) as boolean
'a2 输入中央子午线，以度.分形式输入，如115度30分则输入115.30; 起算数据l0
'f2 以度小数形式输入经度值, l
'e2 以度小数形式输入纬度值，b
's2 计算结果，横坐标y
't2 计算结果，纵坐标x
'投影带号计算　n=[l/6]+1 如：测得经度103.xxxx，故n=[103.x/6]+1=17+1=18
'中央经线经度　l0 = n*6-3 = [l/6]*6+3

dim b2 as double
'dim g2 as double
dim h2 as double
dim i2 as double
dim j2 as double
dim k2 as double
dim l2 as double
dim m2 as double
dim n2 as double
dim o2 as double
dim p2 as double
dim q2 as double
dim r2 as double

b2 = int(a2) + (int(a2 * 100) - int(a2) * 100) / 60 + (a2 * 10000 - int(a2 * 100) * 100) / 3600
'把l0化成度(a2)
'g2 = f2 - b2 ' l -l0
'h2 = g2 / 57.2957795130823 '化作弧度
h2 = (f2 - b2) / 57.2957795130823 '将经差的单位化为弧度
i2 = tan(e2 / 57.2957795130823) 'tan (b)
j2 = cos(e2 / 57.2957795130823) ' cos (b)
k2 = 0.006738525415 * j2 * j2
l2 = i2 * i2
m2 = 1 + k2
n2 = 6399698.9018 / sqr(m2)
o2 = h2 * h2 * j2 * j2
p2 = i2 * j2
q2 = p2 * p2
r2 = (32005.78006 + q2 * (133.92133 + q2 * 0.7031))
s2 = ((((l2 - 18) * l2 - (58 * l2 - 14) * k2 + 5) * o2 / 20 + m2 - l2) * o2 / 6 + 1) * n2 * (h2 * j2)
s2 = s2 + 18500000 '在计算的基础上加上了“带号”（18）和“东移”（500km）
'计算结果,横坐标y
t2 = 6367558.49686 * e2 / 57.29577951308 - p2 * j2 * r2 + ((((l2 - 58) * l2 + 61) * _
o2 / 30 + (4 * k2 + 5) * m2 - l2) * o2 / 12 + 1) * n2 * i2 * o2 / 2
'计算结果，纵坐标x
'msgbox "pts2=" & s2 & " pt t2=" & t2

bl2xy = true
end function

[brightyang]

// GaussBL2xy.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "CoorTrans.h"
#include <iostream>

using namespace std;

void main(int argc, char* argv[])
{
double MyL0 = 108; //中央子午线
double MyB = 33.44556666; //33 du 44 fen 55.6666 miao
double MyL = 109.22334444; //3度带，109 d 22 m 33.4444 s

PrjPoint_Krasovsky MyPrj;
MyPrj.SetL0(MyL0);
MyPrj.SetBL(MyB, MyL);
double OutMyX; //结果应该大致是:3736714.783, 627497.303
double OutMyY;
OutMyX = MyPrj.x; //正算结果: 北坐标x
OutMyY = MyPrj.y; //结果：东坐标y

//////////////////反算////////////////////////////////////////
double InputMyX = 3736714.783; //如果是独立计算，应该给出中央子午线L0
double InputMyY = 627497.303;
MyPrj.Setxy(InputMyX, InputMyY);
MyPrj.GetBL(&MyPrj.B, &MyPrj.L); //把计算出的BL的弧度值换算为dms形式
double OutputMyB;
double OutputMyL;
OutputMyB = MyPrj.B; //反算结果：B
OutputMyL = MyPrj.L; //反算结果：L

//分析表明，此程序的结果和Coord4.2的转换结果是一样的，只差到毫米级
//原程序有几个问题，1.Pi的值不对。2.SetBL中多了两行错误代码
}


double Dms2Rad(double Dms)
{
double Degree, Miniute;
double Second;
int Sign;
double Rad;
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;
return Rad;
}
double Rad2Dms(double Rad)
{
double Degree, Miniute;
double Second;
int Sign;
double Dms;
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);
return Dms;
}
///////////////////////////////////////////////////
// Definition of PrjPoint
///////////////////////////////////////////////////
bool PrjPoint::BL2xy()
{
double X, N, t, t2, m, m2, ng2;
double sinB, 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;
return true;
}
bool PrjPoint::xy2BL()
{
double sinB, cosB, t, t2, N ,ng2, V, yN;
double preB0, B0;
double eta;
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;
--------------------------------
接下

[brightyang]

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;
return true;
}
bool PrjPoint::SetL0(double dL0)
{
L0 = Dms2Rad(dL0);
return true;
}
bool PrjPoint::SetBL(double dB, double dL)
{
B = Dms2Rad(dB);
L = Dms2Rad(dL);
//B = dB; //我靠，I wana say fuck
//L = dL; //del it !
BL2xy();
return true;
}
bool PrjPoint::GetBL(double *dB, double *dL)
{
*dB = Rad2Dms(B);
*dL = Rad2Dms(L);
return true;
}
bool PrjPoint::Setxy(double dx, double dy)
{
x = dx;
y = dy;
xy2BL();
return true;
}
bool PrjPoint::Getxy(double *dx, double *dy)
{
*dx = x;
*dy = y;
return true;
}
///////////////////////////////////////////////////
// Definition of PrjPoint_IUGG1975
///////////////////////////////////////////////////
PrjPoint_IUGG1975::PrjPoint_IUGG1975() //在类外定义构造成员函数，要加上类名和域限定符 ::
{
a = 6378140;
f = 298.257;
e2 = 1 - ((f - 1) / f) * ((f - 1) / f);
e12 = (f / (f - 1)) * (f / (f - 1)) - 1;
A1 = 111133.0047; //这几个A是什么意思？
A2 = -16038.5282;
A3 = 16.8326;
A4 = -0.0220;
}
PrjPoint_IUGG1975::~PrjPoint_IUGG1975() //析构函数，释放构造函数占用的内存
{
}
///////////////////////////////////////////////////
// Definition of PrjPoint_Krasovsky
///////////////////////////////////////////////////
PrjPoint_Krasovsky::PrjPoint_Krasovsky()
{
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;
}
PrjPoint_Krasovsky::~PrjPoint_Krasovsky()
{
}



//CoorTrans.h 文件
#ifndef _COORTRANS_H_INCLUDED
#define _COORTRANS_H_INCLUDED

#include "stdafx.h"
#include "math.h"

const double PI = 3.141592653589793;

double Dms2Rad(double Dms);
double Rad2Dms(double Rad);

class PrjPoint
{
public:
double L0; // 中央子午线经度
double B, L; // 大地坐标
double x, y; // 高斯投影平面坐标
public:
bool BL2xy();
bool xy2BL();

protected:
double a, f, e2, e12; // 基本椭球参数
double A1, A2, A3, A4; // 用于计算X的椭球参数

public:
bool SetL0(double dL0);
bool SetBL(double dB, double dL);
bool GetBL(double *dB, double *dL);
bool Setxy(double dx, double dy);
bool Getxy(double *dx, double *dy);
};

class PrjPoint_Krasovsky : virtual public PrjPoint //类的继承，声明基类是 PrjPoint，虚基类
{
public:
PrjPoint_Krasovsky(); //定义构造函数，用来初始化。(函数名和类名相同)
~PrjPoint_Krasovsky();
};

class PrjPoint_IUGG1975 : virtual public PrjPoint
{
public:
PrjPoint_IUGG1975();
~PrjPoint_IUGG1975();
};

#endif /* ndef _COORTRANS_H_INCLUDED */

[hongqi162]

需要我给你转换成delphi的吗？太简单了你自己弄吧

[imho888]

这段代码我也有，但是算出来不太一样。
问题如下：
1.需要输入中央子午线，具体是多少呢，怎么可以算出来？
2.

//================================
s2 = s2 + 18500000 '在计算的基础上加上了“带号”（18）和“东移”（500km）
'计算结果,横坐标y
//======================================
带号我通过
'投影带号计算　n=[l/6]+1 如：测得经度103.xxxx，故n=[103.x/6]+1=17+1=18
这个算出来是20，但是东移是什么，一定是500KM还是要通过计算？



望解答。

[hongqi162]

我只能做到这些，你说的那些概念性的东西，早就忘光了


uses math;


function bl2xy(var a2:double;var f2:double;var e2:double;var s2:double;var
t2:double):Boolean;
{
'a2 输入中央子午线，以度.分形式输入，如115度30分则输入115.30; 起算数据l0
'f2 以度小数形式输入经度值, l
'e2 以度小数形式输入纬度值，b
's2 计算结果，横坐标y
't2 计算结果，纵坐标x
'投影带号计算　n=[l/6]+1 如：测得经度103.xxxx，故n=[103.x/6]+1=17+1=18
'中央经线经度　l0 = n*6-3 = [l/6]*6+3
}
var
//dim g2 as double
b2,h2,i2,j2,k2,l2,m2,n2,o2,p2,q2,r2:double;
begin
b2 := int(a2) + (int(a2 * 100) - int(a2) * 100) / 60 + (a2 * 10000 - int(a2 * 100) * 100) / 3600;
//'把l0化成度(a2)
//'g2 = f2 - b2 ' l -l0
//'h2 = g2 / 57.2957795130823 '化作弧度
h2:= (f2 - b2) / 57.2957795130823;// '将经差的单位化为弧度
i2:= tan(e2 / 57.2957795130823);// 'tan (b)
j2:= cos(e2 / 57.2957795130823);// ' cos (b)
k2:= 0.006738525415 * j2 * j2;
l2:= i2 * i2;
m2:= 1 + k2;
n2:= 6399698.9018 / sqr(m2);
o2:= h2 * h2 * j2 * j2;
p2:= i2 * j2;
q2:= p2 * p2;
r2:= (32005.78006 + q2 * (133.92133 + q2 * 0.7031));
s2:= ((((l2 - 18) * l2 - (58 * l2 - 14) * k2 + 5) * o2 / 20 + m2 - l2) * o2 / 6 + 1) * n2 * (h2 * j2);
s2:= s2 + 18500000;// '在计算的基础上加上了“带号”（18）和“东移”（500km）
//'计算结果,横坐标y
t2:= 6367558.49686 * e2 / 57.29577951308 - p2 * j2 * r2 + ((((l2 - 58) * l2 + 61) * o2 / 30 + (4 * k2 + 5) * m2 - l2) * o2 / 12 + 1) * n2 * i2 * o2 / 2;
//'计算结果，纵坐标x
//'msgbox "pts2=" & s2 & " pt t2=" & t2
Result:= true
end;

[imho888]

哎，继续求解中。。。。。。。。。。。。。。。。。。。。

[wyffyw2000]

学习，关注，研究

[dzq138]

留个印...

[dzq138]

留个印...

[mei2235]

学习学习，这个问题把我折磨死了要

[simle221]

关注，学习

[小雷阵雨]

汗~啊 正遇到这个问题

[jessezappy]

都在搞这个，先试用一下2楼的算法，不过想问一句，没有带入放大级别计算吗？因为Googlemap换算时是有放大级别的，Z值。

[llsshh1985]

[Quote=引用 5 楼 imho888 的回复:]

这段代码我也有，但是算出来不太一样。
问题如下：
1.需要输入中央子午线，具体是多少呢，怎么可以算出来？
2.

//================================
s2 = s2 + 18500000 '在计算的基础上加上了“带号”（18）和“东移”（500km）
'计算结果,横坐标y……
[/Quote]
中央子午线问题：3度带的中央子午线是这样算的，比如东经112度，中央子午线为111，带号为37。3度带中央子午线1度带为3度，范围是东经1.5到4.5。每度带加3，以此类推。那114度就是38度带的中央子午线。38度范围是112.5——115.5。6度带分法：1度带从东经0度到6度，中央子午线是3度。那17度带的范围是96到102，中央子午线是东经99度。以此类推。

[jax_lee]

学习。。