求助:经纬度坐标和高斯坐标转换

Hawker101 2004-05-07 04:26:43
在将经纬度坐标转换成高斯坐标时(用的是北京1954参考球和6度带投影法), 发现在跨带时Y坐标的增量不是平滑的,而是跳跃式的,
比如:
(32.30, 118) -> (3486844, 20595004)
(32.30, 119) -> (3488144, 20690022) DelatY = 95018
(32.30, 120) -> (3490312, 20785066) DelatY = 95044
(32.30, 121) -> (3488144, 21309977) DelatY = 524911 (???)

为何会出现这种情况,应该是连续的值?
希有大虾告知,小弟不胜感激!



...全文
226 2 点赞 打赏 收藏 举报
写回复
2 条回复
切换为时间正序
当前发帖距今超过3年,不再开放新的回复
发表回复
Hawker101 2004-05-08
这是小弟用的高斯投影(6度带)的代码,各位看看还有问题.

double CFlyMapView::N(double B)
{
double BB, AX, BX, E, ES;
AX = 6378245.0;
BX = 6356863.0188;
E = sqrt(AX * AX - BX * BX) / AX;
ES = E * sin(B);
BB = AX / sqrt(1 - ES * ES);
return BB;
}

double CFlyMapView::S(double B)
{
double SS, S1, S2, S3;
S1 = 6367558.497 * B - 16036.4803 * sin(2 * B);
S2 = 16.8281 * sin(4 * B)- 0.02198 * sin(6 * B);
S3 = 0.000031 * sin(8 * B);
SS = S1 + S2 + S3;
return SS;

}
/***********************************************************************

JwdToGauss - convert the GPS coordinate to the gauss coordinate

Entry:
fa - the latitude of GPS (in)
leda - the longitude of the GPS (in)
xx - the x vector of the gauss coordinate (out)
yy - the y vector of the gauss coordinate (out)

Exit: None

/**********************************************************************/
void CFlyMapView::JwdToGauss(double fa, double leda, double &xx, double &yy)
{
double pi = 3.1415926;
double AX, BX, E1, B, L,N1, S1, AT, T, X, Y;
int numOfProjectStrap = (int)(leda/6) + 1;
leda = leda + 3 - numOfProjectStrap * 6.0;

AX = 6378245.0;
BX = 6356863.0188;
E1 = (AX * AX- BX * BX) / pow(BX, 2.0);
B = fa * pi / 180.0;
L = leda * pi / 180.0;
N1 = N(B);
S1 = S(B);
AT = E1 * pow(cos(B), 2.0);
T = tan(B);
X = S1 + L* L* N1 * sin(B) * cos(B)/2 + pow(L, 4.0)
* N1 * sin(B) * pow(cos(B), 3.0) *
(5- T * T + 9 * AT + 4 * AT * AT) / 24;
Y = L * N1 * cos(B) + pow(L, 3.0) * N1 * pow(cos(B), 3.0)
* (1- T * T + AT) / 6 + pow(L, 5.0) * N1 *
pow(cos(B), 5.0) * (5 - 18 * T * T + pow(T, 4.0)) / 120;
Y = Y + 500 * 1000.0 + (numOfProjectStrap) * 1000.0 * 1000.0;

xx = X;
yy = Y;
}

void CFlyMapView::GaussToJwd(double x, double y, double &B, double &L)
{
double e2, e12, tmp, bf, nf, Nf, tf,fa;
double yNf,yy;
int l0,n;
double fa0 =(0.1570460641219e-6); double b2 =(0.25184647783e-2);
double b4 =(0.36998873e-5);
double b6 =(0.74449e-8);
double b8 =(0.1828e-10);
double a =6378245;
double b =6356863.02;
double pi =3.1415926;

e2=(a-b)*(a+b)/a/a; //¼ÆËãÆ«ÂÊ
e12=(a-b)*(a+b)/b/b;
n=(int)(y/1000.0/1000.0); //¼ÆËã´øºÅ
//n=(int)(y/1000.0);
yy=y-(n*1000.0+500.0)*1000.0; //¼ÆËãyÖµ
l0=n*6-3; //¼ÆËãÖÐÑë¾­¶È

fa=fa0*x; //¼ÆËãbf
bf=fa+b2*sin(2.0*fa)+b4*sin(4.0*fa)+b6*sin(6.0*fa)+b8*sin(8.0*fa);

tf=tan(bf); //¼ÆËãtf,nf,Nf
nf=sqrt(e12)*cos(bf);
Nf=sqrt(1.0-e2*pow(sin(bf),2.0));
Nf=a/Nf;

yNf=yy/Nf; //¼ÆËãB,L(ÒÔ¶ÈΪµ¥Î»)
tmp=1.0-yNf*yNf*(5.0+3.0*tf*tf+nf*nf-9.0*tf*tf*nf*nf)/12.0;
tmp=tmp+pow(yNf,4.0)*(61.0+90.0*tf*tf+45.0*pow(tf,4.0))/360.0;
tmp=tmp*yNf*yNf*(1.0+nf*nf)*tf/2.0;
B=(bf-tmp)*180.0/pi;
tmp=1.0-yNf*yNf*(1.0+2.0*tf*tf+nf*nf)/6.0;
tmp=tmp+pow(yNf,4.0)*(5.0+28.0*tf*tf+24.0*pow(tf,4.0)+6.0*nf*nf+8.0*nf*nf*tf*tf)/120.0;
L=(yy/(Nf*cos(bf))*tmp)*180.0/pi;
L=L+l0;

}
  • 打赏
  • 举报
回复
klbt 2004-05-07
计算误差?关注.
  • 打赏
  • 举报
回复
相关推荐
发帖
地理信息系统
加入

1814

社区成员

它是一种特定的十分重要的空间信息系统。它是在计算机硬、软件系统支持下,对整个或部分地球表层(包括大气层)空间中的有关地理分布数据进行采集、储存、管理、运算、分析、显示和描述的技术系统。
申请成为版主
帖子事件
创建了帖子
2004-05-07 04:26
社区公告
暂无公告