package com.towery.smoke.data;
import java.text.DecimalFormat;
import java.util.HashMap;
import java.util.Map;
public class GpsToBj {
public static void main(String[] args) {
Map<String, String> gpStoBJ = GPStoBJ("39.9084152661928","116.502275001177");
String gpsx = gpStoBJ.get("bjx");
String gpsy = gpStoBJ.get("bjy");
System.out.println("gpsx="+gpsx+",gpsy="+gpsy);
}
/**
* GPS TO BJ
* @param pointx
* @param pointy
* @return
*/
public static Map<String,String> GPStoBJ(String pointx,String pointy){
Map<String,String> map = new HashMap<String, String>();
try {
DecimalFormat df = new DecimalFormat("#.##");
double L0 = 117;
double h = 150;
double L = Double.parseDouble(pointy);//经度
double B = Double.parseDouble(pointx);;//纬度
double[] bj = WGS842BJ(L, B, L0, h);
String pointxsde = df.format(bj[1]);
String pointysde = df.format(bj[0]);
map.put("bjx", pointxsde);
map.put("bjy", pointysde);
} catch (Exception e) {
e.printStackTrace();
}
return map;
}
/**
* WGS84与BJ54相互转换
* gps to bj
* @param L
* @param B
* @param L0
* @param h
* @return
*/
public static double[] WGS842BJ(double L, double B, double L0, double h){
double[] result = new double[2];
//第一步将WGS84的LB转换为WGS84的空间直角坐标KX84、KY84、KZ84
double Ta = 6378137; //6378245 //--椭球体长半轴-54=6378245;WGS84=6378137
double Tb = 6356752.3142; //6356863.0188 //--椭球体短半轴-54=6356863.0188;WGS84=6356752.3142
B = deg(B);
L0 = deg(L0);
L = deg(L);
double E = Math.pow((1 - Math.pow(Tb / Ta, 2)), 0.5); //-- 第一偏心率
double N = Ta / ( Math.pow(1 - E * E * Math.sin(B) * Math.sin(B), 0.5)); // '--卯酉圈曲率半径
double KX84 = (N + h) * Math.cos(B) * Math.cos(L);
double KY84 = (N + h) * Math.cos(B) * Math.sin(L);
double KZ84 = (N * (1 - E * E) + h) * Math.sin(B);
int iType = 7;
double X0 = 0.0; //平移X0
double Y0 = 0.0; //平移Y0
double Z0 = 0.0; //平移Z0
double tx = 0.0; //X轴旋转角度TX
double ty = 0.0; //Y轴旋转角度TY
double tz = 0.0; //Z轴旋转角度TZ
double k = 0.0; //缩放尺度K
if (iType == 7){
//第二步将WGS84空间直角坐标通过七参数转换为54空间直角坐标KX54、KY54、KZ54
X0 = 123.072; //平移X0
Y0 = 279.5004; //平移Y0
Z0 = 27.0206; //平移Z0
tx = -0.000959959166667; //X轴旋转角度TX
tx = deg(tx); //调用子过程将度换算到弧度
ty = 0.000614433888889; //Y轴旋转角度TY
ty = deg(ty); //调用子过程将度换算到弧度
tz = -0.001131175277778; //Z轴旋转角度TZ
tz = deg(tz); //调用子过程将度换算到弧度
k = -0.000014693; //缩放尺度K
} else {
//第二组三参数
X0 = -36.398; //平移X0
Y0 = 73.535; //平移Y0
Z0 = 71.968; //平移Z0
tx = 0; //X轴旋转角度TX
tx = deg(tx); //调用子过程将度换算到弧度
ty = 0; //Y轴旋转角度TY
ty = deg(ty); //调用子过程将度换算到弧度
tz = 0; //Z轴旋转角度TZ
tz = deg(tz); //调用子过程将度换算到弧度
k = 0; //缩放尺度K
}
double KX54 = X0 + (1 + k) * (1 * KX84 + tz * KY84 - ty * KZ84);
double KY54 = Y0 + (1 + k) * (-tz * KX84 + 1 * KY84 + tx * KZ84);
double KZ54 = Z0 + (1 + k) * (ty * KX84 - tx * KY84 + 1 * KZ84);
//第三步将54空间直角坐标反算到54大地坐标B54、L54
Ta = 6378245; // '--椭球体长半轴-54=6378245;WGS84=6378137
Tb = 6356863.0188; // '--椭球体短半轴-54=6356863.0188;WGS84=6356752.3142
E = Math.pow(1 - Math.pow(Tb / Ta, 2), 0.5);// '-- 第一偏心率
double L54 = Math.atan(KY54 / KX54);
double tt = KZ54 * E * E;
double sinB = (KZ54 + tt) / Math.pow(KX54 * KX54 + KY54 * KY54 + Math.pow(KZ54 + tt, 2), 0.5);
int myNum = 1;
double deltsinB = 1;
while (deltsinB > 0.0000001 || myNum <= 10){
myNum = myNum + 1;
double sinB0 = sinB;
N = Ta / (Math.pow(1 - E * E * sinB * sinB,0.5));// '--卯酉圈曲率半径
tt = N * E * E * sinB;
sinB = (KZ54 + tt) / Math.pow(KX54 * KX54 + KY54 * KY54 + Math.pow(KZ54 + tt,2),0.5);
deltsinB = Math.abs(sinB - sinB0);
}
double B54 = Math.asin(sinB);
N = Ta / (Math.pow(1 - E * E * sinB * sinB, 0.5));
L54 = L54 * 180 / 3.1415926536;
B54 = B54 * 180 / 3.1415926536;
if (L54 < 0){
L54 = 180 + L54;
}
if (B54 < 0){
B54 = 180 + B54;
}
//第四步将54大地坐标B54、L54进行高斯投影正算得到54平面坐标X54、Y54
//调用子过程将度分秒换算到弧度
L54 = deg(L54);
B54 = deg(B54);
Ta = 6378245;// ' 6378137'--椭球体长半轴-54=6378245;WGS84=6378137
Tb = 6356863.0188;// '6356752.3142 '--椭球体短半轴-54=6356863.0188;WGS84=6356752.3142
E = Math.pow(1 - Math.pow(Tb / Ta, 2),0.5);// '-- 第一偏心率
double e1 = Math.pow(Math.pow(Ta / Tb, 2) - 1, 0.5);// ' -- 第二偏心率
N = (Ta * Ta / Tb) / (1 + e1 * e1 * Math.cos(B54) * Math.cos(B54));// - -卯酉圈曲率半径
double t = Math.tan(B54) * Math.tan(B54);
double C = e1 * e1 * Math.cos(B54) * Math.cos(B54);
double AA = (L54 - L0) * Math.cos(B54);
double M = Ta * ((1 - Math.pow(E, 2) / 4 - 3 * Math.pow(E,4) / 64 - 5 * Math.pow(E,6) / 256) * B54 - (3 * Math.pow(E,2) / 8 + 3 * Math.pow(E,4) / 32 + 45 * Math.pow(E,6) / 1024) * Math.sin(2 * B54) + (15 * Math.pow(E,4) / 256 + 45 * Math.pow(E,6) / 1024) * Math.sin(4 * B54) - (35 * Math.pow(E,6) / 3072) * Math.sin(6 * B54));
N = Ta / Math.pow(1 - E * E * Math.sin(B54) * Math.sin(B54), 0.5);// '= (Ta * Ta / Tb) / (1 + e1 * e1 * Cos(B54) * Cos(B54)) ^ 0.5
double FE = 500000;
double X54 = M + N * Math.tan(B54) * (AA * AA / 2 + (5 - t + 9 * C + 4 * C * C) * Math.pow(AA, 4) / 24) + (61 - 58 * t + t * t + 270 * C - 330 * t * C) * Math.pow(AA, 6) / 720;
double Y54 = FE + N * (AA + (1 - t + C) * Math.pow(AA, 3) / 6 + (5 - 18 * t + t * t + 14 * C - 58 * t * C) * Math.pow(AA, 5) / 120);
//第五步将54平面坐标X54、Y54换算到北京地方坐标Xdf、Ydf
double dx = (X54 - 4414702.66) / 1000;
double dy = (Y54 - 500000) / 1000;
double spdx = 0;
double spdy = 0;
if (iType == 7){
spdx = 7.345806056;
spdy = 15.48276601;
}
double Xdf = 300202.06 + 1000.01161 * dx + 7.269 * dy - 9.98 * 0.000001 * (dx * dx - dy * dy) + 1.37 * 0.001 * dx * dy + spdx;
double Ydf = 555594.14 - 7.269 * dx + 1000.01161 * dy - 6.84 * 0.0001 * (dx * dx - dy * dy) - 2 * 0.000001 * dx * dy + spdy;
result[0] = Xdf;
result[1] = Ydf;
//result[2] = X54;
//result[3] = Y54;
return result;
}
public static double deg(double a){
double result = a;
//输入的角度是以度为单位的,XXX.XXXXXXX度
double k = Math.signum(a);
result = Math.abs(a);
result = result * Math.PI / 180;
result = result * k;
return result;
}
}