坐标转换-大地转高斯平面平面坐标转换
楔子
以前呢，总感觉大学老师教的都没用，基本上都用不上。得了，这两天碰上了。
我们公司给甲方做了个小程序，因为业务原因必须使用深圳独立坐标系。
一个业务需求-导航。因为我们是获取手机的gps坐标，我们起先是使用的甲方提供的坐标转换服务，发现84大地转到深圳独立坐标（深圳高斯平面坐标）时精度误差极大。于是甲方给我们提供了部分控制点信息，让我们自己去完成坐标的转换，不再使用甲方提供的服务。
一把眼泪，大一学的地图学和测量学早把这些东西忘光了。
一步步来吧，随让自己当初不好好学呢。
数据
49990.201 63839.358 2521448.83 496112.77 113.86475295492812 22.79054504954949
49839.360 97644.067 2511136.78 476638.23 113.86988381063978 22.779607662363563
注：这是这里修改了数据 反正大概长这样 大概位数是这么多 至于具体的值。我已经改得面目全非了
用于四参数计算已经足够了
解决思路
- 这些都是什么鬼数据？
- 这些数据都怎么转换？
- 转换公式是什么?
- 转换代码如何编写？
ps:折腾了快一天后，事实证明，一开始设计的解决问题思路是对的
转换逻辑 84大地坐标系 <–>84平面直角坐标系 <–>深圳平面直角坐标系
数据分析
// 深圳高斯平面数据(至于其他的参数就不给了)
49990.201 63839.358
// 84高斯平面数据2521448.83 496112.77
// 84大地坐标
113.86475295492812 22.79054504954949
软件转换+手工验算
首先呢，理论知识都忘的差不多了，还是先计算清楚把。
这里使用的coord软件
第一步 平面坐标到平面坐标的转换
设置->计算四参数
输入上面两组数据源 单击 -> 计算
单击 ->导出 得到如下四参数
注:这里是深圳平面转84平面的四参数
如何要获取 84平面转深圳平面的四参数 重新计算即可
公式如下
代码如下
/*** 深圳平面坐标转84平面** @param x 纵坐标* @param y 横坐标* @return*/public double[] convert54to84(double x, double y) {double[] coord = {0.0, 0.0};double X;double Y;// Coord gm 控制点计算结果// 核心参数final double DX = 2470000.759939;final double DY = 390000.400211;final double T = -0.0174591152;final double K = 0.999887111016;double a = DX;double b = DY;double t = T;double k = K;// 核心代码X = a + x * k * Math.cos(t) - y * k * Math.sin(t);Y = b + y * k * Math.cos(t) + x * k * Math.sin(t);coord[0] = X;coord[1] = Y;return coord;}
注：这里84平面转深圳平面是一样的 只是核心参数不同罢了
核心参数已经改的面目全非 反正大概数据长这样罢了
第二步 大地坐标于高斯平面坐标的转换
其实这个用专业的术语是 高斯正反算
高斯正算
大地坐标转平面坐标
高斯反算
平面坐标转大地坐标
注：别看了，我正儿八经也没全看懂，主要是在网上找到了现成的代码。具体要看公式的话这里看肯定是不行的，建议下本大地测量学的书，那里会手把手教你如何计算出来。记得大一还是大二的时候老师让我们算过，我可能算过也可能没算过，反正快4年了我忘了。
代码
/*** 由经纬度反算成高斯投影坐标(高斯正算)** @param longitude* @param latitude*/
public double[] BLToGauss(double longitude, double latitude) {int ProjNo = 0;// 带宽int ZoneWide = 3;double longitude1, latitude1, longitude0, X0, Y0, xval, yval;double a, f, e2, ee, NN, T, C, A, M, iPI;// 3.1415926535898/180.0;iPI = 0.0174532925199433;// 84年北京坐标系参数a = 6378137;f = 1.0 / 298.2572236;ProjNo = 0;longitude0 = ProjNo * ZoneWide ;longitude0 = 114 * iPI;// 经度转换为弧度longitude1 = longitude * iPI;// 纬度转换为弧度latitude1 = latitude * iPI;e2 = 2 * f - f * f;ee = e2 * (1.0 - e2);NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude1) * Math.sin(latitude1));T = Math.tan(latitude1) * Math.tan(latitude1);C = ee * Math.cos(latitude1) * Math.cos(latitude1);A = (longitude1 - longitude0) * Math.cos(latitude1);M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude1 - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.sin(2 * latitude1) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude1));xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120);yval = M + NN * Math.tan(latitude1) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720);X0 = 1000000 * (ProjNo ) + 500000;Y0 = 0;xval = xval + X0;yval = yval + Y0;return new double[] { xval, yval };
}/*** 由高斯投影坐标反算成经纬度 84坐标系** @param X* @param Y* @return*/
public double[] GaussToBL(double X, double Y) {int ProjNo;// 带宽int ZoneWide;double[] output = new double[2];// latitude0,double longitude1, latitude1, longitude0, X0, Y0, xval, yval;double e1, e2, f, a, ee, NN, T, C, M, D, R, u, fai, iPI;//3.1415926535898/180.0;iPI = 0.0174532925199433;a = 6378137;f = 1 / 298.2572236;// 3度带宽ZoneWide = 3;// 查找带号ProjNo = 0;// 中央经线longitude0 = 114 * Math.PI / 180;X0 = ProjNo * 1000000L + 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);// 转换为度 DDoutput[0] = longitude1 / iPI;output[1] = latitude1 / iPI;return output;
}注：我在这里埋了个坑，网上找的代码是6度带的。我这里需要的是三度带的，而且我的业务区域很小，只需要114度的中央经线即可，我就偷懒写死了。有会计算三度带的小伙伴可以@我。我从书上看的三度带的计算公式很简单，问题是套在代码里不好使。
回答：今天研究了一下，发现我的y坐标少了两位。所以我套用网上的代码不好使。
我这里的y坐标（横坐标）496112.77
真实的y坐标(横坐标)应该是 38496112.77
Es6版本高斯正反算坐标转换 和 四参数平面坐标转换
注意：四参数平面坐标转换 已经将核心参数改的面目全非了 根据自己的控制点 填写 const参数
/*** 深圳54平面坐标转84平面** @param x 纵坐标* @param y 横坐标* @return*/
function convert54GaussTo84Gauss (x, y) {// Coord gm 控制点计算结果const DX = 2470000.000const DY = 390000.000const T = -0.01900const K = 0.90000000let coord = [0.0, 0.0]let Xlet Ylet a = DXlet b = DYlet t = Tlet k = KX = a + x * k * Math.cos(t) - y * k * Math.sin(t)Y = b + y * k * Math.cos(t) + x * k * Math.sin(t)coord[0] = Xcoord[1] = Yconsole.log('转换前：' + x + ',' + y + '\t转换后：' + coord)return coord
}/*** 84高斯平面坐标转深圳54高斯平面坐标** @param x 纵坐标* @param y 横坐标* @return*/
function convert84GaussTo54Gauss (x, y) {// Coord gm 控制点计算结果const DX = -2470000.000const DY = -410000.000const T = 0.01900const K = 0.90000000let coord = [0.0, 0.0]let Xlet Ylet a = DXlet b = DYlet t = Tlet k = KX = a + x * k * Math.cos(t) - y * k * Math.sin(t)Y = b + y * k * Math.cos(t) + x * k * Math.sin(t)coord[0] = Xcoord[1] = Yconsole.log('转换前：' + x + ',' + y + '\t转换后：' + coord)return coord
}/*** 由经纬度反算成高斯投影坐标** @param longitude 经度* @param latitude 纬度*/
function convert84BLToGauss (longitude, latitude) {// 带宽const ZONE_WITH = 3// 带号// 三度带计算公式 带号 = （经度 + 1.5度） / 带宽// 注意 3度带和6度带的 带号和中央经线的计算方式不同const PROJ_NO = Math.floor(((longitude + 1.5) / ZONE_WITH))// 中央经线（弧度制）// 三度带计算公式： 中央经线(弧度制) = 带号 * 带宽 * (π/180)const longitude0 = (PROJ_NO * ZONE_WITH) * (Math.PI / 180)// 长半径（84）const a = 6378137// 偏率（84）const f = 1 / 298.2572236let X0, Y0, xval, yvallet e2, ee, NN, T, C, A, M// 经度转换为弧度longitude = longitude * (Math.PI / 180)// 纬度转换为弧度latitude = latitude * (Math.PI / 180)e2 = 2 * f - f * fee = e2 * (1.0 - e2)NN = a / Math.sqrt(1.0 - e2 * Math.sin(latitude) * Math.sin(latitude))T = Math.tan(latitude) * Math.tan(latitude)C = ee * Math.cos(latitude) * Math.cos(latitude)A = (longitude - longitude0) * Math.cos(latitude)M = a * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256) * latitude - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2 / 1024) * Math.sin(2 * latitude) + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024) * Math.sin(4 * latitude) - (35 * e2 * e2 * e2 / 3072) * Math.sin(6 * latitude))xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120)yval = M + NN * Math.tan(latitude) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720)X0 = 1000000 * PROJ_NO + 500000Y0 = 0xval = xval + X0yval = yval + Y0let coord = [xval, yval]console.log('转换前：' + longitude * 180 / Math.PI + ',' + latitude * 180 / Math.PI + '\t转换后：' + coord)return coord
}/*** 由高斯投影坐标反算成经纬度 84坐标系** @param x 纵坐标* @param y 横坐标* @return*/
function convert84GaussToBL (x, y) {// 带宽const ZONE_WITH = 3// 带号// 三度带计算公式 带号 = 横坐标前两位// 注意 3度带和6度带的 带号和中央经线的计算方式不同const PROJ_NO = Math.floor(y / 1000000)// 中央经线（弧度制）// 三度带计算公式： 中央经线(弧度制) = 带号 * 带宽 * (π/180)const longitude0 = (PROJ_NO * ZONE_WITH) * (Math.PI / 180)// 长半径（84）const a = 6378137// 偏率（84）const f = 1 / 298.2572236let longitude, latitude, X0, Y0, xval, yvallet e1, e2, ee, NN, T, C, M, D, R, u, faiX0 = 0Y0 = PROJ_NO * 1000000 + 500000xval = x - X0yval = y - Y0e2 = 2 * f - f * fe1 = (1.0 - Math.sqrt(1 - e2)) / (1.0 + Math.sqrt(1 - e2))ee = e2 / (1 - e2)M = xvalu = 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 = yval / NN// 计算经度(Longitude) 纬度(Latitude)longitude = 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)latitude = 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)// 转换为度let coord = [0.0, 0.0]coord[0] = longitude * 180 / Math.PIcoord[1] = latitude * 180 / Math.PIconsole.log('转换前：' + x + ',' + y + '\t转换后：' + coord)return coord
}
console
这里参考了一些其他的博文：
https://blog.csdn.net/u010534192/article/details/77196064 坐标转换流程与公式 七参数 四参数
https://blog.csdn.net/m0_37970224/article/details/79677868 坐标高斯正反算
坐标转换-大地转高斯平面平面坐标转换相关推荐
- 经纬度与高斯-克吕格平面坐标转换
原作者链接:https://blog.csdn.net/jianyi7659/article/details/7583339 前言 支持将地理坐标(经纬度坐标)转换到高斯-克吕格投影下的平面坐标,如北 ...
- 大地坐标和高斯平面坐标转换
图幅理论面积与图斑椭球面积计算公式及要求 一. 图幅理论面积计算公式 (1) 式中: a-椭球长半轴(单位:米),α-椭球扁率,b-椭球短半轴(单位:米). е²﹦(a²﹣b²)/a². A﹦1﹢(3 ...
- Arcgis经纬度到平面坐标转换
从网上下载的DEM文件,导入Arcgis中发现,显示的是经纬度值(117.121 40.164),如下图.与平面坐标有很大偏差. 尝试了各种办法,最后终于找到了解决办法: 1.设置DEM数据文件的坐标 ...
- 空间直角坐标系、大地坐标系、平面坐标系、高斯平面直角坐标系（转）
本篇学习了空间直角坐标系.大地坐标系.平面坐标系.高斯平面直角坐标系.这个个坐标系有时很容易弄混淆! ( 一)空间直角坐标系 空间直角坐标系的坐标原点位于参考椭球的中心,Z轴指向参考椭球的北极 ...
- 通过IDL编辑器实现空间直角到高斯平面坐标系的转化
目录 空间直角坐标系: 大地坐标系: 高斯平面直角坐标系: 转化方法过程细节讲述: 坐标系简介: 空间直角坐标系: 空间直角坐标系,如图所示,空间中的任意点坐标用(X,Y,Z)表示,坐标原点位于地球质 ...
- android 卫星地图,推荐一款亲测好用，可显示卫星地图，高斯平面直角坐标和计算图幅编号等功能的安卓定位导航软件～...
推荐一款亲测好用,可显示卫星地图,高斯平面直角坐标和计算图幅编号等功能的安卓定位导航软件-步行者坐标导航. 一.软件下载 小米应用商店搜:步行者坐标导航或 http://appcdn.wapx.cn/ ...
- 坐标转换-两坐标系间平面坐标转换(附软件下载）
用测量助理进行同基准条件下两个坐标系之间平面直角坐标转换,原理是通过计算四参数,再用四参数分别计算每个点转换后坐标,具体操作步骤如下: 1.打开测量助理软件,左侧栏点击选择二维坐标转换,右侧栏将出现待 ...
- 经纬度BL换算到高斯平面直角坐标XY的简单算法
找到这个经纬度BL换算到高斯平面直角坐标XY的简单算法(她用EXCEL算的,正好是一个个步骤,翻过来就是了): 从经纬度BL换算到高斯平面直角坐标XY(高斯投影正算),或从XY换算成BL(高斯投影反算 ...
- c语言大地坐标转换空间坐标,大地坐标与空间直角坐标转换_C程序
大地坐标与空间直角坐标转换的C程序 #include #include double HD(double a,double b,double c) {b=b+c/60; a=a+b/60; a=a/1 ...
- 坐标转换 计算机图形学_计算机图形学的转换类型
坐标转换 计算机图形学 什么是转型? (What is Transformation?) Transformation refers to the mathematical operations or ...
最新文章
- android之lint警告This Handler class should be static or leaks might occur
- java 连接两个arraylist,java – 在两个线程之间共享一个ArrayList？
- 数据可视化组队学习：《Task02 - 艺术画笔见乾坤》笔记
- 大鱼吃光小鱼，绝不可能!盘点2016存储行业发生的大事件
- 爬虫简单入门：第一个简单爬虫
- python中如何删除字典中的元素_python中字典删除元素
- socket编程 —— 非阻塞socket (转)---例子已上传至文件中
- 整数快速幂（原理+模板）
- if 语句 写了return 报错
- Linux下通过ssh上传、下载文件或者文件夹
- Python快捷键大全（PyCharm常用）
- ajax 的data,ajax请求的data数据格式
- ArcGIS导出shape地图边界点数据
- 为什么我偏爱用GitHub来写书？
- matlab生成word文档
- 青春使命网页制作html,青春使命句子
- 古诗词鉴赏，断句的重要性
- 目前比较流行的网站设计风格有哪些
- Preftest测试
- 【HTML5】H5新标签大实例