最近忙，今天上来一看，GIS版很冷清嘛，我发点热帖暖暖手（经纬度坐标到上海坐标系的转换）(0分)

吕雪松 · 2001-12-02

//Lat：纬度，注：是度.度格式，不是度.分格式
//Lon：经度，同上。
procedure ConvertCoordinate(Lat,Lon:do
uble;var xx,yy :do
uble);
const R2D=57.2957795131;
TOLAT=(31.0+(14.0+7.55996/60.0)/60.0)/R2D;
TOLON=(121.0+(28.0+1.80651/60.0)/60)/R2D;
REARTH=6371006.84;
var
Horizerr,Frlat,Frlon,gcdist,gcbrg :do
uble;
clatf,clatt,slatf,slatt,dlon,cdlon,sdlon,sdist,cdist,sbrg,cbrg,temp :do
uble;
begin

// Frlat:=lat/R2D;
// Frlon:=lon/R2D;

// FrLat := (StrToFloat(Copy(Lat,1,2))+StrToFloat(Copy(Lat,3,7))/60)/R2D;
// FrLon := (StrToFloat(Copy(Lon,1,3))+StrToFloat(Copy(Lon,4,7))/60)/R2D;

FrLat := (Trunc(Lat)+(Lat - Trunc(Lat))/0.6)/R2D;
FrLon := (Trunc(Lon)+(Lon - Trunc(Lon))/0.6)/R2D;

clatt:=cos(Frlat);
clatf:=cos(TOLAT);
slatt:=sin(Frlat);
slatf:=sin(TOLAT);
dlon:=-TOLON+Frlon;
cdlon:=cos(dlon);
sdlon:=sin(dlon);
cdist:=slatf*slatt+clatf*clatt*cdlon;
temp:=sqr((clatt*sdlon))+sqr((clatf*slatt-slatf*clatt*cdlon));
sdist:=sqrt(abs(temp));
if((abs(sdist)>0.0000001) or (abs(cdist)>0.0000001)) then

gcdist :=arctan2(sdist,cdist)
else
gcdist :=0.0;
sbrg :=sdlon*clatt;
cbrg :=(clatf*slatt-slatf*clatt*cdlon);
if((abs(sbrg)>0.0000001) or (abs(cbrg)>0.0000001)) then
begin

temp := arctan2(sbrg,cbrg);
while (temp < 0.0)do
temp := temp+2.0*PI;
gcbrg :=temp;
end
else

gcbrg :=0.0;

Horizerr :=gcdist*REARTH;

xx := Horizerr*sin(gcbrg);
yy := Horizerr*cos(gcbrg);
end;

dingbaosheng · 2001-12-02

kem · 2001-12-02

吕雪松，您好。
我现在有上海的道路地图（由于项目需要，由某官方单位提供），是以人民广场为原点的XY坐标（米单位），要想把GPS接受的坐标数据与该地图匹配起来的话，请问，可以您这个函数吗？？
关键：
1.该地图以人民广场为原点
2.该地图以XY坐标系统，米单位（即在AutoCAD中，画一米的线条，在实际的地图上就是距离1米）
由于本人不是地理信息出身，所以，请您解释一下我的疑问，好吗？
多谢！！！

吕雪松 · 2001-12-02

正好了，就是用这个函数。
上海坐标系就是以米为单位。

kem · 2001-12-03

多谢您的函数，真是雪中送炭。
请问，这个函数的误差大致是多少呢？？

吕雪松 · 2001-12-03

我用的时候，没有感觉到有误差。到上海很边上的如宝山松江等地方，GPS定位的车辆轨迹
和上海坐标系的道路中心线吻合得都很好，说明没有误差扩散的趋势。

kem · 2001-12-03

我的图是从AutoCAD-->DXF--->MapInfo，这样过来的，但愿不会有大的误差。
吕雪松，您是不是也采用AutoCAD-->DXF-->***呢？？
在基础地图制作的精度上，您能给些建议吗？？
多谢！！
java1977@21cn.com

吕雪松 · 2001-12-03

数据误差主要来自于你的基础图的测绘误差和手工数字化误差。
数据格式转换不会增加地图误差，但可能会丢失部分属性和制图数据；但坐标系转换
会导致误差的！
我的数据是ESRI系列的，ArcInfo/ArcView

kem · 2001-12-03

多谢！
以后没有问题的话，还可能要请教您。
能给我六个e-mail吗？

吕雪松 · 2001-12-03

六个Email？我没有这么多！

leisure.lv@263.net

YB_unique · 2001-12-05

//经纬度坐标到上海坐标系的转换
呵呵！有没有经纬度坐标到长沙坐标系的转换？？？

我也灌水！
通过地区和经度计算距离。

USES MATH;
var
StartLat:do
uble;

StartLong:do
uble;

EndLat:do
uble;

EndLong:do
uble;

fPhimean:do
uble;

fdLambda:do
uble;

fdPhi:do
uble;

fAlpha:do
uble;

fRho:do
uble;

fNu:do
uble;

fR:do
uble;

fz:do
uble;

fTemp:do
uble;

Distance:do
uble;

Bearing:do
uble;

const
D2R:do
uble = 0.017453;

R2D:do
uble = 57.295781;

a:do
uble = 6378137.0;

b:do
uble = 6356752.314245;

e2:do
uble = 0.006739496742337;

f:do
uble = 0.003352810664747;

begin

fdLambda := (StartLong - EndLong) * D2R;
fdPhi := (StartLat - EndLat) * D2R;
fPhimean := ((StartLat + EndLat) / 2.0) * D2R;
fTemp := 1 - e2 * (Power(Sin(fPhimean),2));
fRho := (a * (1 - e2)) / Power(fTemp, 1.5);
fNu := a / (Sqrt(1 - e2 * (Sin(fPhimean) * Sin(fPhimean))));
fz := Sqrt(Power(Sin(fdPhi/2.0),2)+Cos(EndLat*D2R)*Cos(StartLat*D2R)*Power(Sin(fdLambda/2.0),2))
fz := 2 * ArcSin(fz);
fAlpha := Cos(EndLat * D2R) * Sin(fdLambda) * 1 / Sin(fz);
fAlpha := ArcSin(fAlpha);
fR := (fRho * fNu) / ((fRho * Power(Sin(fAlpha),2)) + (fNu * Power(Cos(fAlpha),2)));
Distance := (fz * fR);
if((StartLat < EndLat) and (StartLong < EndLong)) then

Bearing := Abs(fAlpha * R2D)
else
if ((StartLat < EndLat) and (StartLong > EndLong)) then

Bearing := 360 - Abs(fAlpha * R2D)
else
if ((StartLat > EndLat) and (StartLong > EndLong)) then

Bearing := 180 + Abs(fAlpha * R2D)
else
if ((StartLat > EndLat) and (StartLong < EndLong)) then

Bearing := 180 - Abs(fAlpha * R2D);
end;

GBlueMan · 2001-12-09

关注

ljx · 2002-01-24

这是我们公司的位置E110.289722 N25.278333，怎么改你的函数呢？

beachboy · 2002-01-25

你们公司不在上海，这要看你需要什么坐标系统，如果仅仅是需要x，y坐标，那可以直接作Gauss投影就可以了
这里是C的源码：
我是搞GPS的，GPS和GIS不分家吗，各位看官可交流，yong.john@163.com
#include
#include
#include
#include "CoorTrans.h"
////////////////////////////////////////////
// Common functions
////////////////////////////////////////////
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 = 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 * t
2 + 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;

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,do
uble dL)
{
B = Dms2Rad(dB);

L = Dms2Rad(dL);

B = dB;

L = dL;

BL2xy();

return TRUE;

}
BOOL PrjPoint::GetBL(double *dB,do
uble *dL)
{
*dB = Rad2Dms(B);

*dL = Rad2Dms(L);

return TRUE;

}
BOOL PrjPoint::Setxy(double dx,do
uble dy)
{
x = dx;

y = dy;

xy2BL();

return TRUE;

}
BOOL PrjPoint::Getxy(double *dx,do
uble *dy)
{
*dx = x;

*dy = y;

return TRUE;

}
///////////////////////////////////////////////////
// Definition of PrjPoint_IUGG1975
///////////////////////////////////////////////////
PrjPoint_IUGG1975:

rjPoint_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;

A2 = -16038.5282;

A3 = 16.8326;

A4 = -0.0220;

}
PrjPoint_IUGG1975::~PrjPoint_IUGG1975()
{
}
///////////////////////////////////////////////////
// Definition of PrjPoint_Krasovsky
///////////////////////////////////////////////////
PrjPoint_Krasovsky:

rjPoint_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
constdo
uble PI = 3.14159265353846;

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,do
uble dL);

BOOL GetBL(double *dB,do
uble *dL);

BOOL Setxy(double dx,do
uble dy);

BOOL Getxy(double *dx,do
uble *dy);

};

class PrjPoint_Krasovsky : virtual public PrjPoint
{
public:
PrjPoint_Krasovsky();

~PrjPoint_Krasovsky();

};

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

~PrjPoint_IUGG1975();

};

double Dms2Rad(double Dms);

double Rad2Dms(double Rad);

#endif /* ndef _COORTRANS_H_INCLUDED */

pqx · 2002-02-07

经纬度坐标到北京54坐标系如何转换？
北京54坐标系如何转换成WGS84经纬度坐标系?

coronor · 2002-02-24

你的图是什么格式呢?
如果是arcinfo的格式,可以用arctools进行转换,但是如果想利用别人的工具
首先得先知道几个参数,比如经差等

kisse · 2002-03-22

我想问一下吕雪松大虾，关于gps中各坐标转换的关系是怎样的？

吕雪松 · 2002-03-23

你得知道你要转换的两种坐标系的投影公式，才能建立二者之前的转换关系。

minifly · 2002-04-19

我用了 beachboy 提供的原码，但结果不对呀！
我需要可以把经纬度（beijing1954）和xy(高斯)做相互转换。
参考《地图投影学》中的公式已经解决了经纬度--〉xy，但xy-->经纬度就没有公式参考了
难道只有反推经纬度--〉xy的公式了吗（好麻烦）？

beachboy · 2002-04-22

minifly:
经纬度到Gauss和Gauss到经纬度是一样的，你没找到公式吗？
超星书库里有大地测量学这本书，里面很详细的内容都有～～～～
在看看吧～～～～

最近忙，今天上来一看，GIS版很冷清嘛，我发点热帖暖暖手（经纬度坐标到上海坐标系的转换）(0分)

吕雪松

Unregistered / Unconfirmed

dingbaosheng

Unregistered / Unconfirmed

kem

Unregistered / Unconfirmed

吕雪松

Unregistered / Unconfirmed

kem

Unregistered / Unconfirmed

吕雪松

Unregistered / Unconfirmed

kem

Unregistered / Unconfirmed

吕雪松

Unregistered / Unconfirmed

kem

Unregistered / Unconfirmed

吕雪松

Unregistered / Unconfirmed

YB_unique

Unregistered / Unconfirmed

GBlueMan

Unregistered / Unconfirmed

ljx

Unregistered / Unconfirmed

beachboy

Unregistered / Unconfirmed

pqx

Unregistered / Unconfirmed

coronor

Unregistered / Unconfirmed

kisse

Unregistered / Unconfirmed

吕雪松

Unregistered / Unconfirmed

minifly

Unregistered / Unconfirmed

beachboy

Unregistered / Unconfirmed