座標轉換-大地轉高斯平面&平面座標轉換
阿新 • • 發佈：2018-12-07
楔子
以前呢，總感覺大學老師教的都沒用，基本上都用不上。得了，這兩天碰上了。
我們公司給甲方做了個小程式，因為業務原因必須使用深圳獨立座標系。
一個業務需求-導航。因為我們是獲取手機的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
*/ 注：這裡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);
// 轉換為度 DD
output[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.000
const DY = 390000.000
const T = -0.01900
const K = 0.90000000
let coord = [0.0, 0.0]
let X
let Y
let a = DX
let b = DY
let t = T
let 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
console.log('轉換前：' + x + ',' + y + '\t轉換後：' + coord)
return coord
}
/**
* 84高斯平面座標轉深圳54高斯平面座標
*
* @param x 縱座標
* @param y 橫座標
* @return
*/
function convert84GaussTo54Gauss (x, y) {
// Coord gm 控制點計算結果
const DX = -2470000.000
const DY = -410000.000
const T = 0.01900
const K = 0.90000000
let coord = [0.0, 0.0]
let X
let Y
let a = DX
let b = DY
let t = T
let 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
console.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.2572236
let X0, Y0, xval, yval
let e2, ee, NN, T, C, A, M
// 經度轉換為弧度
longitude = longitude * (Math.PI / 180)
// 緯度轉換為弧度
latitude = latitude * (Math.PI / 180)
e2 = 2 * f - f * f
ee = 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)