经纬度 to 平面直角坐标转换
2015-04-15 15:20阅读：
在大地测量学中，坐标系分为两大类：地心坐标系和参心坐标系。
地心坐标系是坐标系原点与地球质心重合的坐标系，参心坐标系是坐标系原点位于参考椭球体中心，但不与地球质心重合的坐标系。
我国使用的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_AX
IS, 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
