大地座標BLH轉平面座標xyh(高斯投影座標正算) Java版
阿新 • • 發佈：2018-11-19
技術背景
做過位置資料處理的小夥伴基本上都會遇到座標轉換，而基於高斯投影原理的大地座標轉平面座標就是其中一種座標轉換，座標轉換的目的就是方便後面資料的處理工作，大地座標轉高斯平面座標常用的有兩種，即3°帶和6°帶，具體採用哪種根據實際情況而定。
計算原理
6°帶帶號n與相應的中央子午線L0經度的關係為：
3°帶帶號n’與相應的中央子午線L0’經度的關係為：
設參考橢球的長半軸為 a，第一偏心率為 e，並令：
設中央子午線的經度為 L0，再記：
則高斯投影正算公式為：
其中：
沒學過測繪的同學對以上原理不是很理解，也很正常，大家可以查閱相關《大地測量學》書籍，具體武大版還是礦大版的，沒什麼區別。
具體實現
具體實現平臺依然是IBM的Eclipse軟體，程式語言為Java，下面是以3°帶為例，進行高斯投影座標正算，具體內容請看程式碼
1 package package1; 2 3 public class BLH_xyh { 4 5 public static double a = 6378137; 6 public static double e = Math.sqrt(0.0066943799013); 7 public static double scale_wide = 3; 8 publicstatic double π = 3.14159265358979323846; 9 10 public static void main(String[] args) { 11 // TODO 自動生成的方法存根 12 Point3d xyh = the_coordinates_are_counting(36.0307523111,120.184664478,9.7065); 13 System.out.println(xyh.getX()+","+xyh.getY()+","+xyh.getZ()); 14 } 15 public static Point3d the_coordinates_are_counting(doubleB,double L,double H){ 16 double A_ = 1 17 +3*e*e/4 18 +45*e*e*e*e/64 19 +175*e*e*e*e*e*e/256 20 +11025*e*e*e*e*e*e*e*e/16384 21 +43659*e*e*e*e*e*e*e*e*e*e/65536; 22 double B_ = 3*e*e/4+15*e*e*e*e/16 23 +525*e*e*e*e*e*e/512 24 +2205*e*e*e*e*e*e*e*e/2048 25 +72765*e*e*e*e*e*e*e*e*e*e/65536; 26 double C_ = 15*e*e*e*e/64 27 +105*e*e*e*e*e*e/256 28 +2205*e*e*e*e*e*e*e*e/4096 29 +10395*e*e*e*e*e*e*e*e*e*e/16384; 30 double D_ = 35*e*e*e*e*e*e/512 31 +315*e*e*e*e*e*e*e*e/2048 32 +31185*e*e*e*e*e*e*e*e*e*e/13072; 33 /* 34 double E_ = 315*e*e*e*e*e*e*e*e/16384 35 +3465*e*e*e*e*e*e*e*e*e*e/65536; 36 double F_ = 693*e*e*e*e*e*e*e*e*e*e/13072; 37 */ 38 39 double α = A_*a*(1-e*e); 40 double β = -B_*a*(1-e*e)/2; 41 double γ = C_*a*(1-e*e)/4; 42 double δ = -D_*a*(1-e*e)/6; 43 /* 44 double ε = E_*a*(1-e*e)/8; 45 double ζ = -F_*a*(1-e*e)/10; 46 */ 47 48 double C0 = α; 49 double C1 = 2*β+4*γ+6*δ; 50 double C2 = -8*γ-32*δ; 51 double C3 = 32*δ; 52 53 double x,y,sign; 54 double scale_number = Math.floor(L/scale_wide); 55 if(L > (scale_number * scale_wide + scale_wide/2)){ 56 scale_number =scale_number + 1; 57 sign = -1; 58 }else{ 59 sign = 1; 60 } 61 62 double L0 = scale_wide*scale_number; 63 double l = Math.abs(L-L0); 64 B = B*π/180; 65 l = l*π/180; 66 double t = Math.tan(B); 67 double m0 = Math.cos(B)*l; 68 double η = Math.sqrt(e*e*Math.pow(Math.cos(B),2)/(1-e*e)); 69 double N = a/Math.sqrt(1-e*e*Math.pow(Math.sin(B), 2)); 70 71 double X0 = C0*B+Math.cos(B)*(C1*Math.sin(B)+C2*Math.pow(Math.sin(B),3)+C3*Math.pow(Math.sin(B), 5)); 72 73 x = X0 74 +N*t*m0*m0/2 75 +N*t*m0*m0*m0*m0*(5-t*t+9*η*η+4*η*η*η*η)/24 76 +N*t*m0*m0*m0*m0*m0*m0*(61-58*t*t+t*t*t*t)/720; 77 y = N*m0+ 78 N*m0*m0*m0*(1-t*t+η*η)/6+ 79 N*m0*m0*m0*m0*m0*m0*(5-18*t*t+t*t*t*t+14*η*η-58*η*η*t*t)/120; 80 81 y = y*sign+500000; 82 83 double h = H; 84 85 Point3d xyh = new Point3d(x,y,h); 86 return xyh; 87 } 88 }
其中Point3d的定義如下：
1 package package1; 2 3 public class Point3d { 4 private double x; 5 private double y; 6 private double z; 7 public Point3d(double x,double y,double z){ 8 this.x=x; 9 this.y=y; 10 this.z=z; 11 } 12 public double getX() { 13 return x; 14 } 15 public void setX(double x) { 16 this.x = x; 17 } 18 public double getY() { 19 return y; 20 } 21 public void setY(double y) { 22 this.y = y; 23 } 24 public double getZ() { 25 return z; 26 } 27 public void setZ(double z) { 28 this.z = z; 29 } 30 }
測試資料為36.0307523111,120.184664478,9.7065，執行結果如下：
至此結束
致謝
感謝山東科技大學北斗星光創客興趣學習小組的王老師對於原理文件的整理以及鄭** C++程式碼的技術分享！
參考文件
1、山東科技大學”北斗星光創客”興趣學習小組GNSS技術文件