/**
* 空间大地直角坐标->大地坐标
*/
public GeoPoint XYZ_BLH(int ellipse, Point3D point)
{
GeoPoint geoPoint = new GeoPoint();
if(geoPoint == null || point == null)
{
System.out.println("对象为空！");
}
double a = 0,b = 0,e1 = 0,e2 = 0; //a为长半轴，b为短半轴，e1为第一偏心率，e2为第一偏心率的平方
double deta = 0.0000001;
switch(ellipse) //选择椭球体
{
case 0: //WGS84椭球体
{
a = 6378137.0; //长半轴
b = 6356752.3142;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006694384999588;
break;
}
case 1: //北京54椭球
{
a = 6378245.0; //长半轴
b = 6356863.0187730473;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006693421622966;
break;
}
case 2: //西安80椭球
{
a = 6378140.0; //长半轴
b = 6356755.2881575287;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006694384999588;
break;
}
default:break;
}
//double W = Math.sqrt(1.0-e1*e1*Math.pow(Math.sin(point.getX()), 2));
//double N = a/W;
double X = point.getX();
double Y = point.getY();
double Z = point.getZ();
double m = Math.sqrt(Math.pow(X, 2)+Math.pow(Y, 2));
double L = geoPoint.getL(); //大地经度
double B = geoPoint.getB(); //大地纬度
double H = geoPoint.getH(); //大地高
L = Math.atan(Y/X);
if (L < 0)
{
L += Math.PI;
}
double e2_ = e2/(1-e2);
double c = a*Math.sqrt(1+e2_);
double ce2 = c*e2;
double k = 1 + e2_;
double front = Z/m;
double temp = front;
int count = 0; //迭代次数
do
{
front = temp;
m = Math.sqrt(Math.pow(X, 2)+Math.pow(Y, 2));
temp = Z/m + ce2*front/(m*Math.sqrt(k+ Math.pow(front, 2)));
}
while(Math.abs(front-temp)<deta&&count<100000); //是否在误差范围之内
B = Math.atan(temp);//求纬度
if (B<0)
{
B += Math.PI;
}
double W = Math.sqrt(1-e1*e1*Math.sin(B)*Math.sin(B));
double N = a/W;
//N = (a*m - c*c)/(2*b*Z);
System.out.println("N = " + N);
H = m/Math.cos(B)-N;//求高
//H = Z/Math.sin(B)-N*(1.0-e1*e1);
//H = X/(Math.cos(B)*Math.cos(L))-N;
//转换为角度值
L = Math.toDegrees(L);
B = Math.toDegrees(B);
//H = Math.toDegrees(H);
geoPoint.setL(L);
geoPoint.setB(B);
geoPoint.setH(H);
return geoPoint;
}
/**
* 大地坐标->空间大地直角坐标
* @param ellipse 椭球体
* @param geoPoint 转换前的坐标
* @return 返回空间直角坐标
*/
public Point3D BLH_XYZ(int ellipse, GeoPoint geoPoint)
{
Point3D point1 = new Point3D(1,1,1);
if(geoPoint == null || point1 == null)
{
System.out.println("对象为空！");
}
double a = 0,b,e1 = 0,e2 = 0; //a为长半轴，b为短半轴，e1为第一偏心率，e2为第一偏心率的平方
switch(ellipse) //选择椭球体
{
case 0: //WGS84椭球体
{
a = 6378137.0; //长半轴
b = 6356752.3142;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006694384999588;
break;
}
case 1: //北京54椭球
{
a = 6378245.0; //长半轴
b = 6356863.0187730473;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006693421622966;
break;
}
case 2: //西安80椭球
{
a = 6378140.0; //长半轴
b = 6356755.2881575287;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = 0.006694384999588;
break;
}
default:break;
}
double B = geoPoint.getB()*Math.PI/180; //转换为弧度制
double L = geoPoint.getL()*Math.PI/180;
double H = geoPoint.getH();
//计算卯酉圈曲率半径
double N = a/Math.sqrt(1.0-e1*e1*Math.pow(Math.sin(B), 2));
//计算空间直角坐标
double X = (N + H)*Math.cos(B)*Math.cos(L);
double Y = (N + H)*Math.cos(B)*Math.sin(L);
double Z = (N*(1-e1*e1)+ H)*Math.sin(B);
point1.X = X;
point1.Y = Y;
point1.Z = Z;
return point1;
}
/**
* 高斯-克吕格坐标正算（从大地坐标到平面直角坐标）
* @param ellipse 椭球体
* @param zoneWide 带宽(3度或6度)
* @param geoPoint 大地坐标
* @param point 直角坐标
*/
public void BL_xy(int ellipse,int zoneWide,GeoPoint geoPoint,Point point)
{
if(geoPoint == null || point == null)
{
System.out.println("对象为空！");
}
double a = 0,b,e1 = 0,e2 = 0; //a为长半轴，b为短半轴，e1为第一偏心率，e2为第一偏心率的平方
switch(ellipse) //选择椭球体
{
case 0: //WGS84椭球体
{
a = 6378137.0; //长半轴
b = 6356752.3142;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006694384999588;
break;
}
case 1: //北京54椭球
{
a = 6378245.0; //长半轴
b = 6356863.0187730473;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006693421622966;
break;
}
case 2: //西安80椭球
{
a = 6378140.0; //长半轴
b = 6356755.2881575287;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006694384999588;
break;
}
default:break;
}
double B = geoPoint.getB();
double L = geoPoint.getL();
double centerL = 0.0;//投影带中央经线的度数
int n = 0; //投影带的带号（6°带是1-60，3°带是1-120）
double l = 0; //该地所在经度与中央经度的差值
if (6 == zoneWide) //如果是6度分带
{
n = (int)L/6;
if (L%6 > 0)
{
n = n + 1;
}
centerL = 6*n-3;
l = L - centerL;
}
if (3 == zoneWide)//如果是3度分带
{
n = (int)(L-1.5)/3;
if ((L-1.5)%3>0)
{
n = n +1;
}
centerL = 3*n;
l = L - centerL;
}
//转化为弧度
B = B*Math.PI/180;
L = L*Math.PI/180;
l = l*Math.PI/180;
double X;//从赤道起算的子午线弧长
//计算子午线弧长的系数
double A0 = 1.0+3.0/4.0*Math.pow(e1, 2)+45.0/64.0*Math.pow(e1, 4)
+350.0/512.0*Math.pow(e1, 6)+11025.0/16384.0*Math.pow(e1, 8);
double A2 = -3.0/4*Math.pow(e1, 2)+60.0/64*Math.pow(e1, 4)+
525.0/512*Math.pow(e1, 6)+17640.0/16384.0*Math.pow(e1, 8)/2.0;
double A4 = 15.0/64.0*Math.pow(e1, 4)+210.0/512.0*Math.pow(e1, 6)+
8820.0/16384.0*Math.pow(e1, 8)/4.0;
double A6 = -35.0/512.0*Math.pow(e1, 6)+2520.0/16384.0*Math.pow(e1, 8)/6.0;
double A8 = 315.0/16384.0*Math.pow(e1, 8)/8.0;
//计算子午线弧长X
X = a*(1.0-Math.pow(e1, 2))*A0*B + A2*Math.sin(2*B) + A4*Math.sin(4*B)
+ A6*Math.sin(6*B) + A8*Math.sin(8*B);
double t = Math.tan(B);
double anke = e2*Math.cos(B);
double N = a/Math.sqrt(1.0-Math.pow(e1, 2)*Math.pow(Math.sin(B), 2));
//坐标计算
double x = point.getX();
double y = point.getY();
x = X + 1.0/2.0*N*t*Math.pow(Math.cos(B), 2)*Math.pow(l, 2)+
1.0/24.0*N*t*(5.0-Math.pow(t, 2)+9*Math.pow(anke, 2)+4*Math.pow(anke, 4))*Math.pow(Math.cos(B), 4)*Math.pow(l, 4)
+ 1.0/720*N*t*(61.0-58*Math.pow(t, 2)+Math.pow(t, 4)+270.0*Math.pow(anke, 2)-330.0*Math.pow(anke, 2)*Math.pow(t, 2))*Math.pow(Math.cos(B), 6);
y = N*Math.cos(B)*l + 1.0/6*N*(1.0-Math.pow(t, 2)+Math.pow(anke, 2))*Math.pow(Math.cos(B), 3)*Math.pow(l, 3)+
1.0/120*N*(5.0-18*t*t+Math.pow(t, 4)+14*Math.pow(anke, 2)-58*anke*anke*t*t)*
Math.pow(Math.cos(B), 3)*Math.pow(l, 5);
y += 500000.0;//加上500KM
point.setX(x);
point.setY(y);
}
/**
* 高斯平面直角坐标转换为大地坐标
* @param ellipse 椭球体
* @param zoneWide 带宽
* @param point 平面坐标
* @param geoPoint 大地坐标
*/
public void xy_BL(int ellipse,double centerL,Point point,GeoPoint geoPoint)
{
if(geoPoint == null || point == null)
{
System.out.println("对象为空！");
}
double a = 0,b,e1 = 0,e2 = 0; //a为长半轴，b为短半轴，e1为第一偏心率，e2为第二偏心率
switch(ellipse) //选择椭球体
{
case 0: //WGS84椭球体
{
a = 6378137.0; //长半轴
b = 6356752.3142;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006694384999588;
break;
}
case 1: //北京54椭球
{
a = 6378245.0; //长半轴
b = 6356863.0187730473;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006693421622966;
break;
}
case 2: //西安80椭球
{
a = 6378140.0; //长半轴
b = 6356755.2881575287;//短半轴
e1 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/a;
e2 = Math.sqrt(Math.pow(a, 2)-Math.pow(b, 2))/b;
//e2 = 0.006694384999588;
break;
}
default:break;
}
double x = point.getX();
double y = point.getY()-500000;
double Bf = 0.0;//底点纬度
double B0 = 1.0+3.0/4.0*e1*e1 + 45.0/64.0*Math.pow(e1, 4)+350.0/312.0*Math.pow(e1, 6)+
11025.0/16384.0*Math.pow(e1, 8) + 43659.0/65536.0*Math.pow(e1, 10);
Bf = x/(a*(1.0-Math.pow(e1, 2))*B0);
double tf = Math.tan(Bf);
double Mf = a*(1.0-e1*e1)/Math.sqrt(Math.pow(1.0-Math.pow(e1, 2), 3));
double Nf = (a/Math.sqrt(1.0-Math.pow(e1, 2)))/Math.sqrt(1.0+Math.pow(e2*Math.cos(Bf),2));
double anke = e2*Math.cos(Bf);
double B = geoPoint.getB();
double L = geoPoint.getL();
//开始坐标计算
B = Bf - tf/(2*Mf*Nf*Math.cos(Bf))*Math.pow(y, 2) + tf/(24*Mf*Math.pow(Nf, 3))*
(5.0+3.0*tf*tf+anke*anke-9.0*anke*anke*tf*tf)*Math.pow(y, 4)-
1.0/(720.0*Math.pow(Nf, 5)*Math.cos(Bf))*(61.0+90.0*tf*tf+45*Math.pow(tf,4))*Math.pow(y, 6);
L = y/(Nf*Math.cos(Bf)) - (1.0+2*tf*tf+anke*anke)*Math.pow(y, 3)/(6.0*Math.pow(Nf, 3)*Math.cos(Bf))
+ (5.0+28*tf*tf+24*Math.pow(tf, 4)+6*anke*anke+8*anke*anke*tf*tf)*Math.pow(y, 5)/(120.0*Math.pow(Nf, 5)*Math.cos(Bf));
B = B*180/Math.PI;
L = L*180/Math.PI;
L = L + centerL;//中央经线+经度差
geoPoint.setB(B);
geoPoint.setL(L);
}这些代码可能有什么不合理的地方，还有就是由空间直角坐标转换为大地坐标的时候，大地高计算有很大偏差，一致找不出来是什么原因
