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)<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);  } 這些代碼可能有什麼不合理的地方，還有就是由空間直角座標轉換為大地座標的時候，大地高計算有很大偏差，一致找不出來是什麼原因