从零开始搭建一个GIS开发小框架(四)——扩展功能:CGCS2000坐标转WGS84坐标

2022-12-01 15:48:15 浏览数 (2)

1

背景知识

Knowledge

2000国家大地坐标系,是我国当前最新的国家大地坐标系,英文名称为China Geodetic Coordinate System 2000,英文缩写为CGCS2000

2000国家大地坐标系的原点为包括海洋和大气的整个地球的质量中心;2000国家大地坐标系的Z轴由原点指向历元2000.0的地球参考极的方向,该历元的指向由国际时间局给定的历元为1984.0的初始指向推算,定向的时间演化保证相对于地壳不产生残余的全球旋转,X轴由原点指向格林尼治参考子午线与地球赤道面(历元2000.0)的交点,Y轴与Z轴、X轴构成右手正交坐标系。采用广义相对论意义下的尺度。

国家大地坐标系是测制国家基本比例尺地图的基础。根据《中华人民共和国测绘法》规定,中国建立全国统一的大地坐标系统。

建国以来,中国于上世纪50年代和80年代分别建立了1954年北京坐标系和1980西安坐标系,测制了各种比例尺地形图,在国民经济、社会发展和科学研究中发挥了重要作用,限于当时的技术条件,中国大地坐标系基本上是依赖于传统技术手段实现的。54坐标系采用的是克拉索夫斯基椭球体。该椭球在计算和定位的过程中,没有采用中国的数据,该系统在中国范围内符合得不好,不能满足高精度定位以及地球科学、空间科学和战略武器发展的需要。上世纪70年代,中国大地测量工作者经过二十多年的艰巨努力,终于完成了全国一、二等天文大地网的布测。经过整体平差,采用1975年IUGG第十六届大会推荐的参考椭球参数,中国建立了1980西安坐标系,1980西安坐标系在中国经济建设、国防建设和科学研究中发挥了巨大作用。

随着社会的进步,国民经济建设、国防建设和社会发展、科学研究等对国家大地坐标系提出了新的要求,迫切需要采用原点位于地球质量中心的坐标系统(以下简称地心坐标系)作为国家大地坐标系。采用地心坐标系,有利于采用现代空间技术对坐标系进行维护和快速更新,测定高精度大地控制点三维坐标,并提高测图工作效率。

2008年3月,由国土资源部正式上报国务院《关于中国采用2000国家大地坐标系的请示》,并于2008年4月获得国务院批准。自2008年7月1日起,中国将全面启用2000国家大地坐标系,国家测绘局授权组织实施。

————摘自百度

2

坐标转换功能实现

Coordinate transformation function

1

功能描述

现在国家层面的地图资源已经全面采用CGCS2000坐标,所以我们拿到的规划条件里的界址点坐标都是CGCS2000坐标。如果想做一些基于GIS的数据挖掘和分析,需要把数据(例如POI数据)加载到底图上,然后根据POI附带的有价值的数据进行后续计算分析工作,但是从第三方获取到的POI数据一般都是基于wgs84、百度坐标或gcj02坐标,所以我们必须要实现坐标的转换功能。专业软件如Arcgis可以实现但是买不起,只能自己开发了。因为我们搭建的框架采用的OpenCycleMap是wgs84坐标,所以我的应用场景是把CGCS2000坐标的资源数据转换成wgs84坐标,以这个转换为例做一个案例,其它的转换都是类似的思路。

2

效果演示

项目规划条件的实际位置(CGCS2000坐标):

把地块的CGCS2000坐标复制到demo程序中:

点击换算按钮得到wgs84坐标,并根据wgs84坐标绘制地块图形,可以看到换算后的位置还比较准确:

3

核心功能代码

Code

这个算法是3度带的,3度带中央经线=3度带带号*3,举个例子,比如:武汉的带号是38,所以武汉的中央子午线=38x3=114,别的城市是多少可以自己百度查一下。

double X = lng - 38000000; 这个38000000也是根据带号来的,注意这个细节。

代码语言:javascript复制
/// <summary>
/// CGCS2000toWGS84
/// lat是纬度,lng是经度
/// </summary>
/// <param name="lng">经度</param>
/// <param name="lat">维度</param>
/// <returns></returns>
public double[] CGCS2000toWGS84(double lng, double lat)
{
    double Y = lat;               //CGCS2000坐标
    double X = lng - 38000000;    //CGCS2000坐标

    double[] output = new double[2];
    double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
    //NN曲率半径,测量学里面用N表示
    //M为子午线弧长,测量学里用大X表示 
    //fai为底点纬度,由子午弧长反算公式得到,测量学里用Bf表示
    //R为底点所对的曲率半径,测量学里用Nf表示
    double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;
    iPI = 0.0174532925199433; //3.1415926535898 / 180.0;
    a = 6378137.0; f = 1 / 298.257222101; //CGCS2000坐标系参数
                                          //a=6378137.0; f=1/298.2572236; //wgs84坐标系参数
    longitude0 = 114.0;//中央子午线
    longitude0 = longitude0 * iPI; //中央经线

    X0 = 500000L;
    Y0 = 0;
    xval = X - X0; yval = Y - Y0; //带内大地坐标
    e2 = 2 * f - f * f;
    e1 = (1.0 - Math.Sqrt(1 - e2)) / (1.0   Math.Sqrt(1 - e2));
    ee = e2 / (1 - e2);
    M = yval;
    u = M / (a * (1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256));
    fai = u   (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * Math.Sin(2 * u)   (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * Math.Sin(
            4 * u)
              (151 * e1 * e1 * e1 / 96) * Math.Sin(6 * u)   (1097 * e1 * e1 * e1 * e1 / 512) * Math.Sin(8 * u);
    C = ee * Math.Cos(fai) * Math.Cos(fai);
    T = Math.Tan(fai) * Math.Tan(fai);
    NN = a / Math.Sqrt(1.0 - e2 * Math.Sin(fai) * Math.Sin(fai));
    R = a * (1 - e2) / Math.Sqrt((1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin(fai) * Math.Sin(fai)) * (1 - e2 * Math.Sin
            (fai) * Math.Sin(fai)));
    D = xval / NN;
    //计算经度(Longitude) 纬度(Latitude)
    longitude1 = longitude0   (D - (1   2 * T   C) * D * D * D / 6   (5 - 2 * C   28 * T - 3 * C * C   8 * ee   24 * T * T) * D
            * D * D * D * D / 120) / Math.Cos(fai);
    latitude1 = fai - (NN * Math.Tan(fai) / R) * (D * D / 2 - (5   3 * T   10 * C - 4 * C * C - 9 * ee) * D * D * D * D / 24
              (61   90 * T   298 * C   45 * T * T - 256 * ee - 3 * C * C) * D * D * D * D * D * D / 720);
    //转换为度 DD
    output[0] = longitude1 / iPI;
    output[1] = latitude1 / iPI;

    return output;
}

4

交流

Communication

如果你有什么GIS的研发项目需要合作或者有什么意见或者建议,请发邮到:heavenleft@foxmail.com

0 人点赞