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