在大地测量学中，坐标系分为两大类：地心坐标系和参心坐标系。
地心坐标系是坐标系原点与地球质心重合的坐标系，参心坐标系是坐标系原点位于参考椭球体中心，但不与地球质心重合的坐标系。
我国使用的1954北京坐标系，1980西安坐标系都属于参心坐标系。GPS中使用的世界大地坐标系WGS－84属于地心坐标系，我国最近开始启用的中国大地坐标系2000（即CGCS2000），也属于地心坐标系。
以上两大类坐标系都有下列几种表达形式：
1.空间大地坐标系，即大地经纬度（B，L，H）形式
2.空间直角坐标系，即三维空间坐标（X，Y，Z）形式
3.投影平面直角坐标系。即二维平面坐标（x，y，h）形式
在工程测量和施工中，我国普遍使用的是1954北京或1980西安的高斯投影平面直角坐标系。
但为满足工程施工精度要求，通常会在测区建立独立的地方坐标系，且独立地方坐标系都能够通过转换公式换算为国家统一的坐标系上，如1954北京坐标系或1980西安坐标系。楼主说的施工图纸上面标的那个是测量坐标可能是
国家平面直角坐标系和独立的地方平面坐标系之一。
如何将GPS坐标转换为XY平面坐标（简易转换）
本文根据《GPS经纬度坐标转平面坐标的简化计算方法及精度分析》这篇文章中的的方法将GPS经纬度坐标转换为以地平面上平面直角坐标系中的X、Y坐标。这个转换方法的前提是在一定的范围以内。具体的转化公式请参考该文，下面是坐标转换的代码：
public class PlaneCoordinate {
private static final double MACRO_AXIS = 6378137; //
赤道圆的平均半径
private static final double MINOR_AXIS = 6356752; //
半短轴的长度，地球两极距离的一半
// 返回Y坐标
private static double turnY(GePoint basePoint, GePoint point)
{
double a = Math.pow(MACRO_AXIS, 2.0);
double b = Math.pow(MINOR_AXIS, 2.0);
double c = Math.pow(Math.tan(basePoint.getLatitude()),
2.0);
double d =
Math.pow(1/Math.tan(basePoint.getLatitude()),2.0);
double x = a/Math.sqrt(a + b*c);
double y = b/Math.sqrt(b + a*d);
c = Math.pow(Math.tan(point.getLatitude()), 2.0);
d = Math.pow(1/Math.tan(point.getLatitude()), 2.0);
double m = a/Math.sqrt(a + b*c);
double n = b/Math.sqrt(b + a*d);
return new PePoint(x, y).distanceBetween(new PePoint(m,
n));
}
// 返回X坐标
private static double turnX(GePoint basePoint, GePoint point)
{
double a = Math.pow(MACRO_AXIS, 2.0);
double b = Math.pow(MINOR_AXIS, 2.0);
double c = Math.pow(Math.tan(basePoint.getLatitude()),
2.0);
double x = a/Math.sqrt(a + b*c);
return x * (point.getLongtitude() -
basePoint.getLongtitude());
}
}
public class GePoint {
private double latitude; // 纬度坐标
private double longtitude; // 经度坐标
public GePoint() {
}
public GePoint(double latitude, double longtitude) {
this.latitude = latitude;
this.longtitude = longtitude;
}
public double getLatitude() {
return 2 * latitude * Math.PI / 360 ;
}
public void setLatitude(double latitude) {
this.latitude = latitude;
}
public double getLongtitude() {
return 2 * longtitude * Math.PI / 360;
}
public void setLongtitude(double longtitude) {
this.longtitude = longtitude;
}
}
另一种方法：
经纬度坐标与高斯坐标的转换代码
//没钱又丑发表于2007-12-6
13:08:00
// double y;
输入参数: 高斯坐标的横坐标，以米为单位
// double x; 输入参数:
高斯坐标的纵坐标，以米为单位
// short DH;
输入参数: 带号，表示上述高斯坐标是哪个带的
// double *L;
输出参数: 指向经度坐标的指针，其中经度坐标以秒为单位
// double *B;
输出参数: 指向纬度坐标的指针，其中纬度坐标以秒为单位
void GaussToGeo(double y, double x, short DH,
double* L, double* B, double LP)
{
double l0;
// 经差
double tf;
// tf =
tg(Bf0),注意要将Bf转换成以弧度为单位
double nf;
// n = y *
sqrt( 1 + etf ** 2) / c, 其中etf = e'**2 * cos(Bf0) ** 2
double
t_l0; // l0，经差，以度为单位
double
t_B0; // B0，纬度，以度为单位
double
Bf0; //
Bf0
double
etf; //
etf,其中etf = e'**2 * cos(Bf0) ** 2
double
X_3;
double PI
= 3.14159265358979;
double
b_e2 = 0.0067385254147;
double b_c
= 6399698.90178271;
X_3 = x /
1000000.00 - 3;
// 以兆米（）为单位
//
对于克拉索夫斯基椭球，计算Bf0
Bf0 =
27.11115372595 + 9.02468257083 * X_3 - 0.00579740442 * pow(X_3,
2)
- 0.00043532572 * pow(X_3, 3)
+ 0.00004857285 * pow(X_3, 4)
+ 0.00000215727 * pow(X_3, 5)
- 0.00000019399 * pow(X_3, 6);
tf =
tan(Bf0 * PI / 180);
// tf =
tg(Bf),注意这里将Bf转换成以弧度为单位
etf = b_e2
* pow(cos(Bf0 * PI / 180), 2); //
etf = e'**2 * cos(Bf) ** 2
nf = y *
sqrt(1 + etf) / b_c; //
n = y * sqrt( 1 + etf ** 2) / c
//
计算纬度，注意这里计算出来的结果是以度为单位的
t_B0 = Bf0
- (1.0 + etf) * tf / PI * (90.0 * pow(nf, 2)
- 7.5 * (5.0 + 3 * pow(tf, 2) + etf - 9 * etf *
pow(tf, 2)) * pow(nf, 4)
+ 0.25 * (61 + 90 * pow(tf, 2) + 45 * pow(tf, 4))
* pow(nf, 6));
//
计算经差，注意这里计算出来的结果是以度为单位的
t_l0 =
(180 * nf - 30 * (1 + 2 * pow(tf, 2) + etf) * pow(nf, 3)
+ 1.5 * (5 + 28 * pow(tf, 2) +
24 * pow(tf, 4)) * pow(nf, 5))
/ (PI * cos(Bf0 * PI /
180));
l0 = (t_l0
* 3600.0);
// 将经差转成秒
if (LP ==
-1000)
{
*L = (double)((DH * 6 - 3) *
3600.0 + l0); // 根据带号计算出以秒为单位的绝对经度，返回指针
}
else
{
*L = LP * 3600.0 + l0;
// 根据带号计算出以秒为单位的绝对经度，返回指针
}
//----------------------------------
*B =
(double)(t_B0 * 3600.0); //
将纬差转成秒，并返回指针
}
// double jd;
输入参数:
地理坐标的经度，以秒为单位
// double wd;
输入参数:
地理坐标的纬度，以秒为单位
// short DH;
输入参数:
三度带或六度带的带号
void GeoToGauss(double jd, double wd, short DH,
short DH_width, double* y, double* x, double LP)
{
double t;
//
t=tgB
double L;
//
中央经线的经度
double l0;
// 经差
double
jd_hd, wd_hd; //
将jd、wd转换成以弧度为单位
double
et2; // et2
= (e' ** 2) * (cosB ** 2)
double N;
// N = C /
sqrt(1 + et2)
double X;
//
克拉索夫斯基椭球中子午弧长
double m;
// m = cosB
* PI/180 * l0
double
tsin, tcos; // sinB,cosB
double PI
= 3.14159265358979;
double
b_e2 = 0.0067385254147;
double b_c
= 6399698.90178271;
jd_hd = jd
/ 3600.0 * PI / 180.0; //
将以秒为单位的经度转换成弧度
wd_hd = wd
/ 3600.0 * PI / 180.0; //
将以秒为单位的纬度转换成弧度
//
如果不设中央经线（缺省参数: -1000），则计算中央经线，
//
否则，使用传入的中央经线，不再使用带号和带宽参数
//L = (DH
- 0.5) * DH_width ;
// 计算中央经线的经度
if (LP ==
-1000)
{
L = (DH - 0.5) * DH_width;
//
计算中央经线的经度
}
else
{
L = LP;
}
l0 = jd /
3600.0 - L;
// 计算经差
tsin =
sin(wd_hd);
// 计算sinB
tcos =
cos(wd_hd);
// 计算cosB
//
计算克拉索夫斯基椭球中子午弧长X
X =
111134.8611 / 3600.0 * wd - (32005.7799 * tsin + 133.9238 *
pow(tsin, 3)
+ 0.6976 *
pow(tsin, 5) + 0.0039 * pow(tsin, 7)) * tcos;
et2 = b_e2
* pow(tcos, 2);
// et2 = (e' ** 2) * (cosB **
2)
N = b_c /
sqrt(1 + et2);
// N = C / sqrt(1 + et2)
t =
tan(wd_hd);
//
t=tgB
m = PI /
180 * l0 * tcos;
// m = cosB * PI/180 *
l0
*x = X + N
* t * (0.5 * pow(m, 2)
+ (5.0 - pow(t, 2) + 9.0 *
et2 + 4 * pow(et2, 2)) * pow(m, 4) / 24.0
+ (61.0 - 58.0 * pow(t, 2) +
pow(t, 4)) * pow(m, 6) / 720.0);
*y = N *
(m + (1.0 - pow(t, 2) + et2) * pow(m, 3) / 6.0
+ (5.0 - 18.0 * pow(t, 2) +
pow(t, 4) + 14.0 * et2
- 58.0 * et2 * pow(t, 2)) * pow(m, 5) /
120.0);
}
来自:http://gisempire.com/blog/user1/28/200712613822.html
加载中…
喜欢