加载中…
个人资料
Qrcode
Qrcode
  • 博客等级:
  • 博客积分:0
  • 博客访问:14,985
  • 关注人气:2
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

经纬度 to 平面直角坐标转换

(2015-04-15 15:20:58)
标签:

it

在大地测量学中,坐标系分为两大类:地心坐标系和参心坐标系。

地心坐标系是坐标系原点与地球质心重合的坐标系,参心坐标系是坐标系原点位于参考椭球体中心,但不与地球质心重合的坐标系。

我国使用的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

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有