/** * 空间大地直角坐标->大地坐标 */ 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); } 这些代码可能有什么不合理的地方,还有就是由空间直角坐标转换为大地坐标的时候,大地高计算有很大偏差,一致找不出来是什么原因 |