|
@@ -0,0 +1,75 @@
|
|
|
+
|
|
|
+
|
|
|
+
|
|
|
+pub fn gauss_proj_cal(longitude : f64, latitude:f64) -> (f64,f64){
|
|
|
+ let i_pi = 0.0174532925199433; ////3.1415926535898/180.0;
|
|
|
+ let zone_wide = 6; ////6度带宽
|
|
|
+ let ap = 6378245.0;
|
|
|
+ let f = 1.0 / 298.3; //54年北京坐标系参数
|
|
|
+ let lnint = longitude as i32;
|
|
|
+
|
|
|
+ let proj_no : i32 = lnint / zone_wide;
|
|
|
+ let mut longitude0 : f64 = proj_no as f64 * zone_wide as f64 + zone_wide as f64/ 2.0;
|
|
|
+ longitude0 = longitude0 * i_pi;
|
|
|
+ // let mut latitude0: f64 = 0.0;
|
|
|
+ let longitude1: f64 = longitude * i_pi; //经度转换为弧度
|
|
|
+ let latitude1: f64 = latitude * i_pi; //纬度转换为弧度
|
|
|
+ let e2: f64 = 2.0 * f - f * f;
|
|
|
+ let ee: f64 = e2 * (1.0 - e2);
|
|
|
+ let nn: f64 = ap / (1.0 - e2 * latitude1.sin()*latitude1.sin()).sqrt();
|
|
|
+ let t: f64= latitude1.tan()*latitude1.tan();
|
|
|
+ let c: f64 = ee * latitude1.cos()*latitude1.cos();
|
|
|
+ let a: f64 = (longitude1 - longitude0)*latitude1.cos();
|
|
|
+ let m: f64 = ap * ((1.0 - e2 / 4.0 - 3.0 * e2*e2 / 64.0 - 5.0 * e2*e2*e2 / 256.0)*latitude1 - (3.0 * e2 / 8.0 + 3.0 * e2*e2 / 32.0 + 45.0 * e2*e2
|
|
|
+ *e2 / 1024.0)*((2.0 * latitude1).sin())
|
|
|
+ + (15.0 * e2*e2 / 256.0 + 45.0 * e2*e2*e2 / 1024.0)*((4.0 * latitude1).sin()) - (35.0 * e2*e2*e2 / 3072.0)*((6.0 * latitude1).sin()));
|
|
|
+ // println!(" lat: {} m: {}",latitude1,m);
|
|
|
+ let mut xval: f64 = nn * (a + (1.0 - t + c)*a*a*a / 6.0 + (5.0 - 18.0 * t + t * t + 72.0 * c - 58.0 * ee)*a*a*a*a*a / 120.0);
|
|
|
+ let mut yval: f64 = m + nn * latitude1.tan()*(a*a / 2.0 + (5.0 - t + 9.0 * c + 4.0 * c*c)*a*a*a*a / 24.0
|
|
|
+ + (61.0 - 58.0 * t + t * t + 600.0 * c - 330.0 * ee)*a*a*a*a*a*a / 720.0);
|
|
|
+ let x0: i32 = 1000000 * (proj_no + 1) + 500000;
|
|
|
+ let y0: i32 = 0;
|
|
|
+ xval = xval + x0 as f64;
|
|
|
+ yval = yval + y0 as f64;
|
|
|
+ (xval,yval)
|
|
|
+}
|
|
|
+
|
|
|
+//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)
|
|
|
+pub fn gauss_proj_inv_cal(X: f64, Y: f64) -> (f64,f64){
|
|
|
+
|
|
|
+ let i_pi: f64 = 0.0174532925199433; ////3.1415926535898/180.0;
|
|
|
+ let a: f64 = 6378245.0;
|
|
|
+ let f: f64 = 1.0 / 298.3; //54年北京坐标系参数
|
|
|
+ ////a=6378140.0; f=1/298.257; //80年西安坐标系参数
|
|
|
+ let zone_wide: i32 = 6; ////6度带宽
|
|
|
+ let proj_no: i32 = (X / 1000000.0) as i32; //查找带号
|
|
|
+ let mut longitude0: f64 = (proj_no as f64 - 1.0) * zone_wide as f64 + zone_wide as f64/ 2.0;
|
|
|
+ longitude0 = longitude0 * i_pi; //中央经线
|
|
|
+ let x0: f64 = proj_no as f64 * 1000000.0 + 500000.0;
|
|
|
+ let y0: f64 = 0.0;
|
|
|
+ let xval: f64 = X - x0;
|
|
|
+ let yval: f64 = Y - y0; //带内大地坐标
|
|
|
+ let e2: f64 = 2.0 * f - f * f;
|
|
|
+ let e1: f64 = (1.0 - (1.0 - e2).sqrt()) / (1.0 + (1.0 - e2).sqrt());
|
|
|
+ let ee: f64 = e2 / (1.0 - e2);
|
|
|
+ let m: f64 = yval;
|
|
|
+ let u: f64 = m / (a*(1.0 - e2 / 4.0 - 3.0 * e2*e2 / 64.0 - 5.0 * e2*e2*e2 / 256.0));
|
|
|
+ let fai: f64 = u + (3.0 * e1 / 2.0 - 27.0 * e1*e1*e1 / 32.0)*(2.0 * u).sin() + (21.0 * e1*e1 / 16.0 - 55.0 * e1*e1*e1*e1 / 32.0)*(
|
|
|
+ 4.0 * u).sin()
|
|
|
+ + (151.0 * e1*e1*e1 / 96.0)*(6.0 * u).sin() + (1097.0 * e1*e1*e1*e1 / 512.0)*(8.0 * u).sin();
|
|
|
+ let c: f64 = ee * fai.cos()*fai.cos();
|
|
|
+ let t: f64 = fai.tan()*fai.tan();
|
|
|
+ let nn: f64 = a / (1.0 - e2 * fai.sin()*fai.sin()).sqrt();
|
|
|
+ let r: f64 = a * (1.0 - e2) / ((1.0 - e2 * fai.sin()*fai.sin())*(1.0 - e2 * fai.sin()*fai.sin())*(1.0 - e2 *
|
|
|
+ fai.sin()*fai.sin())).sqrt();
|
|
|
+ let d: f64 = xval / nn;
|
|
|
+ //计算经度(Longitude) 纬度(Latitude)
|
|
|
+ let longitude1: f64 = longitude0 + (d - (1.0 + 2.0 * t + c)*d*d*d / 6.0 + (5.0 - 2.0 * c + 28.0 * t - 3.0 * c*c + 8.0 * ee + 24.0 * t*t)*d
|
|
|
+ *d*d*d*d / 120.0) / fai.cos();
|
|
|
+ let latitude1: f64 = fai - (nn*fai.tan() / r)*(d*d / 2.0 - (5.0 + 3.0 * t + 10.0 * c - 4.0 * c*c - 9.0 * ee)*d*d*d*d / 24.0
|
|
|
+ + (61.0 + 90.0 * t + 298.0 * c + 45.0 * t*t - 256.0 * ee - 3.0 * c*c)*d*d*d*d*d*d / 720.0);
|
|
|
+ //转换为度 DD
|
|
|
+ let longitude: f64 = longitude1 / i_pi;
|
|
|
+ let latitude: f64 = latitude1 / i_pi;
|
|
|
+ (longitude,latitude)
|
|
|
+}
|