## GIS坐标转换及其编程实现

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