从零开始搭建一个GIS开发小框架（四）——扩展功能：CGCS2000坐标转WGS84坐标
从零开始搭建一个GIS开发小框架（四）——扩展功能：CGCS2000坐标转WGS84坐标
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也是根据带号来的，注意这个细节。
/// <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
腾讯云开发者
